Short-Term Synchronous and Asynchronous Ambient Noise Tomography in Urban Areas: Application to Karst Investigation

Ya Liu; Jianghai Xia; Bo Guan; Chaoqiang Xi; Ling Ning; Hao Zhang

doi:10.1016/j.eng.2025.02.001

Engineering ›› 2025, Vol. 48 ›› Issue (5) :309 -325. DOI: 10.1016/j.eng.2025.02.001

Research Geotechnical Engineering—Article

research-article

Short-Term Synchronous and Asynchronous Ambient Noise Tomography in Urban Areas: Application to Karst Investigation

Author information +

History +

PDF (6337KB)

Abstract

Dense-array ambient noise tomography is a powerful tool for achieving high-resolution subsurface imaging, significantly impacting geohazard prevention and control. Conventional dense-array studies, however, require simultaneous observations of numerous stations for extensive coverage. To conduct a comprehensive karst feature investigation with limited stations, we designed a new synchronous–asynchronous observation system that facilitates dense array observations. We conducted two rounds of asynchronous observations, each lasting approximately 24 h, in combination with synchronous backbone stations. We achieved wide-ranging coverage of the study area utilizing 197 nodal receivers, with an average station spacing of 7 m. The beamforming results revealed distinct variations in the noise source distributions between day and night. We estimated the source strength in the stationary phase zone and used a weighting scheme for stacking the cross-correlation functions ( $C^{1}$ functions) to suppress the influence of nonuniform noise source distributions. The weights were derived from the similarity coefficients between multicomponent $C^{1}$ functions related to Rayleigh waves. We employed the cross-correlation of $C^{1}$ functions ( $C^{2}$ methods) to obtain the empirical Green’s functions between asynchronous stations. To eliminate artifacts in $C^{2}$ functions from higher-mode surface waves in $C^{1}$ functions, we filtered the $C^{1}$ functions on the basis of different particle motions linked to multimode Rayleigh waves. The dispersion measurements of Rayleigh waves obtained from both the $C^{1}$ and $C^{2}$ functions were utilized in surface wave tomography. The inverted three-dimensional (3D) shear-wave (S-wave) velocity model reveals two significant low-velocity zones at depths ranging from 40 to 60 m, which align well with the karst caves found in the drilling data. The method of short-term synchronous–asynchronous ambient noise tomography shows promise as a cost-effective and efficient approach for urban geohazard investigations.

Graphical abstract

Keywords

Seismic interferometry / Surface wave tomography / Asynchronous ambient noise / Geohazards / Seismic ambient noise

Cite this article

Download citation ▾

Ya Liu, Jianghai Xia, Bo Guan, Chaoqiang Xi, Ling Ning, Hao Zhang. Short-Term Synchronous and Asynchronous Ambient Noise Tomography in Urban Areas: Application to Karst Investigation. Engineering, 2025, 48(5): 309-325 DOI:10.1016/j.eng.2025.02.001

登录浏览全文

4963

注册一个新账户忘记密码

1. Introduction

In the early years of 21st century, seismic interferometry theory revealed that performing cross-correlation calculations of ambient noise (referred to as the

C^{1}

methods) can yield the empirical Green’s function (EGF) of subsurface structures [1], [2], [3]. The EGF approximates the Green’s function, which describes the signal received at one point from an impulse source at another point, representing the impulse response of the Earth’s subsurface. Surface wave energy dominates ambient noise, allowing researchers to retrieve surface waves from ambient noise data [4], [5], [6]. The study of ambient noise surface waves spans multiple scales and application scenarios, including detailed crustal and lithosphere characterization [7], [8], [9], [10], [11], [12], [13], global-scale tomography [14], near-surface structure imaging [15], [16], [17], [18], volcanic structure characterization [19], geothermal reservoir imaging [20], [21], [22], and fault structure investigations [23], [24], [25]. In near-surface applications, researchers have reported that high-quality surface wave imaging results can be obtained via relatively short-duration ambient noise recordings (from a few hours to a few minutes) of traffic-induced ground vibrations [18], [26], [27], [28], [29], [30], [31], greatly reducing the conventional duration of ambient noise observations. Ambient noise imaging offers several advantages, including its environmentally friendly nature, cost-effectiveness, and robustness against noise interference. These characteristics make it a highly promising method for applications in densely populated urban areas characterized by complex human activities.

Karst hazards are common geological disasters that often cause surface subsidence, ground collapse, groundwater pollution, and other damage [32]. These karst hazards can significantly impact the safety and stability of human habitation, causing severe economic losses and adverse social consequences [33]. Understanding the development and distribution of karst features is crucial for maintaining social safety and stability, as well as promoting sustainable economic development. Various geophysical methods are used to characterize abnormal regions, including acoustic arrival tomography [34], [35], [36], refraction wave methods [37], [38], [39], resistivity methods [40], [41], electromagnetic methods [42], [43], seismic waveform inversion [44], and multichannel analysis of surface waves [45], [46]. However, the complex urban environment (including restricted land accessibility and noise interference) poses significant challenges to geophysical imaging methods. This study uses ambient noise imaging to characterize underground karst features.

The combination of dense array observations and ambient noise surface wave tomography can provide high-resolution 3D shear-wave (S-wave) velocity structures of the subsurface medium [47], [48], [49]. However, a dense array requires many stations, which can be costly and sometimes challenging to achieve. In recent years, seismic interferometry theory has proposed the use of additional cross-correlations or convolutions on

C^{1}

functions to obtain EGFs [50], [51], [52]. The cross-correlation of the coda of the

C^{1}

functions is referred to as the

C^{3}

methods [50], [53], [54]. The cross-correlation or convolution of the direct surface waves on the

C^{1}

functions is referred to as the

C^{2}

methods [55], [56], [57]. The secondary cross-correlations or convolutions of

C^{1}

functions make it possible to utilize asynchronous observations [58].

In this study, we design a new synchronous–asynchronous ambient noise observation system (Fig. 1) in an abandoned parking lot to achieve short-term (two rounds of observations, approximately 24 h per round) dense array observations (average station spacing of 7 m) with a limited number of stations (197 nodal receivers) for high-resolution detection and imaging of underground karst features. We first analyze the complex spatiotemporal characteristics of urban ambient noise and propose a new

C^{1}

functions weighted stacking method to suppress the impact of nonuniform noise source distributions on the accuracy of

C^{1}

functions. Then, we demonstrate that the high-mode surface waves in

C^{1}

functions can be suppressed by exploiting the differences in the particle motions of fundamental- and higher-mode Rayleigh waves, thereby improving the accuracy of asynchronous

C^{2}

functions. Combining the Rayleigh wave dispersion measurement data from both the

C^{1}

and

C^{2}

functions, we obtain a high-resolution 3D S-wave velocity model of the subsurface via an ambient noise surface wave tomography method. The low-velocity zones identified below 40 m correspond well with the karst caves indicated by the drilling data. Finally, we examine the model resolution and the improvements achieved through the implementation of asynchronous observation.

2. Data

2.1. Design of synchronous and asynchronous seismic arrays

The study area is an abandoned parking lot in Hangzhou, China, on top of which an office building is planned (Fig. 1(a)). The borehole data reveal that many karst caves exist underground (Fig. 1(b)), which threatens the safety of the planned construction. Ambient noise surface wave tomography is used to obtain the spatial distribution of subsurface anomalies and help address karst threats. Notably, a total of 36 lithological boreholes have been drilled within the parking area (Fig. 1(b)), three of which also feature S-wave velocity logging results. These pieces of information serve as auxiliary support to validate our inversion and model interpretations.

To achieve both extensive and intensive coverage over the area via a limited number of seismometers (197 nodal geophones in this study), we designed a new seismic observation system including synchronous and asynchronous seismic arrays (Fig. 1(a)). We conducted two rounds of observations within the confines of the parking lot, with each round lasting approximately 24 h. Each observation round consisted of 6 survey lines, with a spacing of approximately 10 m between the lines and approximately 5 m between the stations along each line. Upon completion of the first observation round, we promptly repositioned the stations to commence the second observation round. Simultaneously, outside the parking lot, we deployed 29 additional stations to record the ambient noise continuously throughout the entire observation period. These additional stations play a crucial role in the subsequent analysis, serving as bridges to reconstruct the surface waves between the two asynchronously observed station sets.

2.2. Ambient noise characteristics

We utilized a total of 197 three-component nodal geophones (with a natural frequency of 5 Hz) and a sampling rate of 500 Hz. In the two rounds of observations, we selected one node for each (outlined by black circles in Fig. 1(a)) to conduct a time–frequency analysis. As illustrated in Fig. 2, since the observations were carried out on two consecutive workdays, their time–frequency characteristics closely resembled each other. The frequency band of ambient noise was quite extensive, with higher energy during the daytime and reduced energy at nighttime, particularly for frequencies above 20 Hz, where the energy decreased significantly during nighttime. In the nighttime time–frequency spectra, numerous vertical white spectral lines were observed, and we tentatively attributed them to sporadic traffic events at night. Additionally, we observed a notable energy reduction during the midday period and speculate that it may be associated with human activities, particularly the reduction in traffic activity during lunchtime and siesta. The time–frequency analysis results indicate that the recorded ambient noise primarily originates from human activities and exhibits noticeable diurnal variations. This observation also reflects the complex time-varying characteristics of anthropogenic seismic sources.

In addition to the temporal and frequency characteristics, we also examined the spatial variations in the seismic sources. The distribution of ambient noise sources affects the retrieval of seismic waves [3], [59], [60], [61], [62]. We employed beamforming analysis to investigate the seismic source distribution using the first-round data of the internal array in the parking lot. By utilizing array data, beamforming can be used to scan the azimuth and velocity of seismic waves passing through the array. It delays and sums the array signals, with stronger summed energy indicating the true azimuth of the seismic waves. The ambient noise sources at different frequencies may originate from different directions. We analyzed the distribution of noise sources at frequencies of 3 and 10 Hz. Furthermore, we derived the temporal variations in the seismic source distribution on the basis of successive short time windows (Fig. 3).

The distributions of the 3 Hz noise sources throughout the entire observation period appear to closely approximate a random and uniform pattern (Fig. 3(a)). This phenomenon could be attributed to the low attenuation of low-frequency seismic waves, suggesting a possible origin from distant sources. In contrast, the energy of the 10 Hz noise source exhibits significant directional variations and notable diurnal variations (Fig. 3(b)). During the daytime, the energy is primarily concentrated between 35° and 120°, as well as between 200° and 300°, corresponding to the locations of nearby intersections in the vicinity of the site. On the basis of our field observations, there are no significant anthropogenic activities other than traffic in these directions. The site is surrounded by office buildings, and the observed energy concentrations are attributed primarily to traffic-induced vibrations at these intersections. At night, the energy becomes more uniformly distributed, without distinct directional patterns. This indicates that the 10 Hz noise sources predominantly originated from nearby road traffic. Notably, we normalize the beamforming analysis results for each time segment to better observe the characteristics of the source distribution; thus, the energy of the seismic sources between different time windows is not directly comparable.

On the basis of the spatial–temporal–frequency analysis of the recorded ambient noise, the ambient noise in urban environments has a highly complex nature. This complexity may challenge the underlying assumption of seismic interferometry, particularly the assumption of a randomly and uniformly distributed source distribution, which is crucial for reliable seismic wave retrieval. Therefore, it is imperative to exercise extra caution in data processing. In the next section, we introduce the methods employed for retrieving surface waves from synchronous and asynchronous records and for conducting surface wave dispersion measurements.

3. Methods

3.1. Calculation of C¹

The theory of stationary-phase approximation revealed that noise sources in stationary-phase zones (SPZs) play a crucial role in the retrieval of seismic waves from ambient noise. In SPZs, the noise sources undergo constructive interference, whereas noise sources outside these zones experience destructive interference. Enhancing the noise sources within the SPZs can significantly improve the accuracy of the retrieved surface waves [63], [64]. In the case of our data, the aforementioned analysis revealed an uneven distribution of noise sources, suggesting that continuing with traditional processing methods and relying on long-term stacking to obtain reliable signals may not be appropriate. Therefore, we shift our focus to preserving only the time segments where the SPZ sources are present, aiming to improve the accuracy of the retrieved surface waves.

The SPZ of surface waves encompasses the direction of the line connecting two stations as well as the surrounding area. In contrast to previous approaches that used array methods to analyze source strengths along the in-line direction, we analyze the in-line strengths by leveraging the similarity of cross-component

C^{1}

functions.

In the ideal far-field plane wave and uniform noise source distribution scenario, the cross-component

C^{1}

function of Rayleigh waves exhibits the following relationship [65]:

(1)

\begin{matrix} ϕ (Z Z) = ϕ (R R) \end{matrix}

(2)

\begin{matrix} |ϕ (Z R) - ϕ (R Z)| = π \end{matrix}

where

ϕ (:)

represents the phase term in the frequency domain. Z and R are the vertical and radial components, respectively. ZR represents the

C^{1}

function between the Z-component of one receiver and the R-component of another receiver. This convention can also be extended to ZZ, RR, and RZ, where ZZ and RR denote the C¹ functions between the same components of two receivers, and RZ represents the C¹ function between the R-component of one receiver and the Z-component of another receiver. One particular point that needs to be emphasized is that all the

C^{1}

functions in Eqs. (1), (2), as well as those in other instances mentioned in the subsequent text, are symmetric components (the average of the causal and time-reversed acausal parts) of the

C^{1}

functions to reduce the effect of the potential nonuniform source distribution [66]. Rayleigh waves are elliptically polarized, with a phase difference between the vertical and radial components. The phase difference between the vertical and radial components of Rayleigh waves can be either

\frac{π}{2}

- \frac{π}{2}

, depending on the type of particle motion, either retrograde or prograde [67]. Retrograde and prograde particle motions describe the movement of particles in seismic waves. Prograde motion indicates that the particle movement direction is consistent with the propagation direction of the surface wave, whereas retrograde motion indicates that it is opposite to the propagation direction. Eqs. (1), (2), however, hold for different types of particle motion when the noise sources are uniformly distributed or distributed along the in-line direction of two receivers. Therefore, we propose considering the similarity between ZZ and RR, as well as ZR and −RZ, as indicators to assess the strength of in-line seismic sources:

(3)

\begin{matrix} W = \max (P (Z Z, R R), P (Z R, - R Z)), if W < δ, W = 0 \end{matrix}

where

P

(…, …) represents the calculation of the Pearson correlation coefficient, W denotes the maximum value among two sets of correlation coefficients (ZZ and RR, as well as ZR and RZ), and

δ

represents the threshold value. If the value of W is below the threshold, then W is set to 0. The computation of the Pearson correlation coefficient involves the removal of the mean from the data and normalization by the standard deviation (SD), making it a suitable measure for assessing signal similarity. The weight is data adaptive, and only a threshold value needs to be set. The closer W is to 1, the stronger the seismic source is along the in-line direction. Conversely, the closer W is to 0, the more likely the seismic sources are to be located at off-line positions. We can obtain W for each ambient noise segment and use it as a weight for the subsequent stacking of

C^{1}

functions. This allows us to select data segments with strong in-line seismic sources for

C^{1}

function stacking, aiming to suppress the influence of nonuniform noise source distributions on the retrieval of surface waves and improve the reliability and quality of the signals.

To validate the effectiveness of the weighted stacking method, we present a simulated data case with a nonuniform distribution of noise sources. We employed a two-layer model as the subsurface medium (Table 1) and simulate one hour of three-component (Z, N, and E) ambient noise data [68], where N and E represent the north–south and east–west components, respectively. A total of 1400 seismic sources were distributed within a circular region, with an inner radius of 500 m and an outer radius of 1500 m (Fig. 4). Specifically, we deliberately introduced additional noise sources in the azimuthal ranges of 120° to 130° and 320° to 330° to create a nonuniform distribution of noise sources. A linear array comprising 25 receivers is positioned in the central area, with a receiver spacing of 5 m between each station. The intensities and orientations of all the point sources were randomly generated, resulting in the presence of both Rayleigh and Love waves in the ambient noise. Importantly, we only simulate fundamental-mode surface waves to ensure that any additional signals observed in the waveforms and dispersion energy are not caused by higher-mode surface waves but rather by the nonuniform distribution of noise sources.

We performed several standard preprocessing steps on the data [66], [69], including data segmentation (3 s long with 50% overlap), mean removal, detrending, and spectral whitening. We subsequently used the leftmost station as the reference station and computed

C^{1}

functions between this reference station and the other stations. As we were dealing with three-component data, we obtained cross-component

C^{1}

functions such as ZN, ZE, and NE, where ZN denotes the cross-correlation between the N and Z components, ZE represents the cross-correlation between the E and Z components, and NE corresponds to the cross-correlation between the N and E components. We applied component rotation to obtain

C^{1}

functions containing Rayleigh waves (vertical and radial components), namely, ZZ, RR, ZR, and RZ. We calculated the weights W for the data segments and performed weighted stacking on the ZZ

C^{1}

functions. Fig. 4(a) illustrates the distribution of sources and receivers, providing spatial context for the data simulation. Fig. 4(b) presents the comparison results between the virtual shot gathers (VSGs) obtained through linear stacking and weighted stacking. The VSG is formed by sorting the cross-correlation pairs on the basis of stations’ distances from the reference station. Compared with the results obtained via the new stacking approach, the VSG of linear stacking results in spurious early arrivals. We performed multichannel dispersion imaging on the VSGs [70]. The results demonstrate that the original stacking method results in a significant number of artifacts in the dispersion spectrum (Fig. 4(c)). In contrast, the dispersion trend of the new stacking method agrees well with the theoretical dispersion curve and is noise free (Fig. 4(d)). The influence of a nonuniform noise source distribution on the retrieval of surface waves has been eliminated through the application of our new stacking method.

For the actual field data, the new method was utilized to obtain the

C^{1}

functions. The data were first downsampled to 100 Hz to increase computational efficiency and then bandpass-filtered in the range of 1–25 Hz. The data preprocessing subsequently followed the procedures used in the synthetic tests, except for setting the segmentation window to 30 s. We calculated the stacking weights W using a threshold

δ

of 0.4. The

C^{1}

functions of different components were all weighted and stacked by the obtained weights. Fig. 5 shows the VSG with the first trace of the first-round data as the reference trace. The conventional linear stacking method results in numerous artifacts, particularly early arrivals before direct surface waves. However, after applying the weighted stacking method, the waveforms become clearer, and the early arrivals caused by the uneven noise source distribution are noticeably suppressed.

We statistically analyzed the distribution of the similarity between different components of the

C^{1}

functions after stacking. The definition of similarity is the same as in Eq. (3), except that no thresholding is applied. Figs. 6(a) and (b) indicate that the conventional linear stacking method failed to suppress the influence of nonuniform noise sources, resulting in low similarity of the

C^{1}

functions. The use of the new stacking method, however, significantly improved the similarity of the stacked

C^{1}

functions. We also calculated the signal-to-noise ratios (SNRs) of the

C^{1}

functions (Figs. 6(c) and (d)). The SNR is calculated as the ratio of the maximum amplitude within the signal window to the root mean square of the data outside the signal window. After applying the weighted stacking method, the number of

C^{1}

functions with an SNR greater than six considerably increases, and this improvement is consistent across multiple components.

3.2. Calculation of C²

The method for obtaining reliable

C^{1}

has been introduced above, which lays a solid foundation for computing

C^{2}

. Before computing

C^{2}

, it is essential to consider the influence of higher-mode surface waves. The multimode surface waves in

C^{1}

could generate fundamental and higher-mode cross-talk artifacts in

C^{1}

. For the VSG in Fig. 5(b), we computed its dispersion spectrum (Fig. 7(a)) and found highly prominent higher modes. Assuming that the fundamental mode and the higher modes involve retrograde and prograde elliptical particle motions, respectively, the higher modes can be suppressed by stacking the ZZ, RR, ZR, and RZ components with appropriate phase corrections, as extensively described in previous studies [65], [67], [71]. A positive

\frac{π}{2}

phase shift is applied to the ZR cross-correlation, and a negative

\frac{π}{2}

phase shift is applied to the RZ cross-correlation. We employ the phase-weighted stacking method [72] to fully leverage their phase coherency, suppressing higher-mode surface waves by stacking the ZZ, RR, and phase-corrected ZR and RZ cross-correlations (Fig. 7(b)). To mitigate the influence of higher modes on

C^{2}

, we utilize the filtered

C^{1}

to calculate

C^{2}

A schematic diagram for computing

C^{2}

is shown in Fig. 8. We treated the external stations located outside the parking area as intermediary stations and calculated these

C^{1}

functions with the first- and second-round stations within the parking area to derive

C^{1}

. We subsequently performed cross-correlation between these two sets of

C^{1}

to obtain

C^{2}

, which were selectively stacked to form the final

C^{2}

estimate, representing the surface waves between the two asynchronous stations. The selection is based on the angle

θ

(Fig. 8) formed between the lines connecting the internal stations and the external station. As indicated in Fig. 8, the station distance for

C^{2}

is denoted as

d_{1}

, while the distances between two

C^{1}

station pairs are denoted as

d_{2}

and

d_{3}

. The actual propagation distance for the signal obtained by cross-correlating

C^{1}

(d_{3} - d_{2})

, which we approximate as

d_{1}

. The larger the

θ

is, the greater the approximation error. In this study, we set the

θ

to 30°, meaning that only

C^{1}

station pairs within

\pm

15° of the in-line direction are used to calculate

C^{2}

. On the basis of the station distances in this study, we calculate that for station pairs at the maximum angular edge, the approximation error is within 5%. Additionally, the

C^{2}

signals obtained within the angle range are stacked, which further reduces the approximation error. This selective stacking enhances the reliability and quality of

C^{2}

[55].

To further illustrate the impact of higher-mode surface waves on

C^{2}

, we compared the VSGs of

C^{2}

(Fig. 9). All asynchronous

C^{2}

pairs were sorted according to interstation distances and averaged at 3 m intervals. All waveforms were bandpass filtered between 2 and 15 Hz to ensure a fairer comparison. We observed that when higher modes in

C^{1}

are not filtered, the corresponding

C^{2}

functions exhibit numerous nonphysical waveforms (Fig. 9(a)). Nevertheless, once higher-mode surface wave filtering is applied to

C^{1}

, the resulting

C^{2}

shows only clear fundamental-mode surface waves (Fig. 9(b)), which ensures the reliability of the subsequent dispersion measurements.

3.3. Dispersion measurement

After the previous steps, we obtained the symmetric components of the

C^{1}

functions between the internal stations of the parking lot, including the

C^{1}

and

C^{2}

functions. Next, we applied the two-station dispersion analysis method [7], [73] to obtain the dispersion energy spectra for any cross-correlation pairs. This method involves performing time–frequency analysis on the

C^{1}

function between two stations and converting the time axis to a velocity axis on the basis of the distance between the stations, thereby obtaining the phase velocity of surface waves between the stations. To address phase ambiguity, which results in multiple peak energy clusters in the dispersion spectrum, during the picking of two-station phase-velocity dispersion curves, we utilized the multichannel dispersion analysis method to obtain a reference dispersion curve and selected dispersion trends near the reference curve in the two-station dispersion spectra. Only

C^{1}

functions with an SNR greater than six are selected for the dispersion curve picking.

To investigate the potential errors in phase velocity measurements obtained via the

C^{2}

method, we computed the

C^{2}

functions between station pairs in the first round of observation and automatically extracted the dispersion curves. We compared the phase velocities (

V_{C^{1}} - V_{C^{2}}

) extracted from the

C^{1}

and

C^{2}

functions along the same paths. At all frequencies where data were available, we calculated the relative velocity deviation (

\frac{V_{C^{2}} - V_{C^{1}}}{V_{C^{1}}} \times 100 %

). Our histogram (Fig. 10(a)) shows that 80% of the velocity differences fall within

\pm 6 %

. We also observed a greater proportion of positive velocity deviations, which is related to the approximation errors in the station paths inherent in the

C^{2}

method, as previously mentioned. The sources of these velocity differences are numerous, including signal quality, imperfect noise source distribution, and other factors. Fortunately, the error results indicate that there is no significant systematic bias in the velocity differences, which partially ensures the reliability of our results. Despite these errors, the

C^{2}

method allows for asynchronous observations, leading to denser station coverage. Fig. 10(b) shows the dispersion curves obtained from two rounds of observations via the

C^{1}

functions and the dispersion curves obtained from asynchronous station pairs via the

C^{2}

functions. The average dispersion curve of

C^{2}

is very close to that of

C^{1}

, indicating the reliability of our

C^{2}

results. We further removed data deviating from the mean velocity by more than

\pm 4

SDs to ensure that our data remained within a reasonable range. The final dispersion curve data used for inversion are shown in Fig. 11. To satisfy the plane-wave approximation theory, we retained only data with wavelengths not exceeding the station spacing [74]. The number of

C^{1}

dispersion measurements sharply decrease in the low-frequency band (Fig. 11(a)), primarily due to limited station spacing, making it challenging to obtain long-wavelength signals. The asynchronous array covers a larger geographical area, with an increased number of station pairs and greater station intervals, effectively boosting long-wavelength dispersion measurements. The phase velocity undergoes significant disturbances at various frequencies, and considerable variation in phase velocity is observed between low and high frequencies (Fig. 11(b)), indicating both the complexity and heterogeneity of the subsurface medium within the study area.

3.4. Surface wave tomography

After obtaining the surface wave dispersion curves, we proceed with surface wave tomography to derive the subsurface 3D S-wave velocity structure. With surface wave traveltime data in hand, we can select any reliable surface wave tomography algorithm and compare different imaging results; however, such comparisons are beyond the scope of this study. In this study, we employed the surface wave dispersion direct inversion method [75] to conduct surface wave tomography. This method uses the fast marching algorithm for ray tracing, accounting for the possibility of surface waves traveling along non-great circle paths, making it more suitable for imaging complex near-surface media. The direct inversion method involves solving the forward problem to compute theoretical traveltimes on the basis of an initial model and then iteratively adjusting the model to minimize the difference between observed and computed traveltimes. It establishes a linear relationship between 3D S-wave velocity perturbations and surface wave traveltime residuals by utilizing surface wave ray paths and one-dimensional (1D) surface wave depth sensitivity kernels. As a result, we can directly obtain the subsurface 3D S-wave velocity structure from surface wave dispersion measurements, eliminating the intermediate step of acquiring two-dimensional (2D) phase/group velocity maps via conventional approaches.

3.5. Workflow summary

The previous sections introduced the key steps involved in synchronous and asynchronous ambient noise surface wave tomography. Fig. 12 presents a complete workflow diagram. For traditional ambient noise surface wave tomography, the main steps include ambient noise data preprocessing,

C^{1}

function calculation, dispersion measurement, and surface wave tomography. In this study, we considered the impact of uneven noise source distributions on the

C^{1}

function during short-term observations and added a step for enhancing the noise sources within the SPZs. We also noted the presence of higher-mode surface waves in the

C^{1}

function, which affects both the

C^{1}

dispersion analysis and the

C^{2}

calculation. Therefore, we incorporated a step for suppressing higher-mode surface waves. To improve the accuracy of

C^{2}

retrieved surface waves, we designed an SPZ

C^{2}

stacking step, excluding

C^{2}

functions outside the SPZ from interfering with surface wave retrieval.

Notably, the workflow we designed serves as a framework, and the specific techniques used for various steps are not explicitly defined. This means that different methods can be applied to each step. For example, the surface wave tomography algorithm used in this study can be replaced with other surface wave tomography algorithms. Similarly, the suppression of higher-mode surface waves in this study is achieved through polarization filtering, but alternative methods, such as suppression based on high-resolution time–frequency analysis [76], can also be explored. The main innovation of this study lies in the design of several additional data processing steps, which are crucial for reliably obtaining

C^{1}

and

C^{2}

functions in short-term observations.

4. Results

4.1. Initial model and depth sensitivity

Model parameterization is conducted by setting equidistant grids in the longitudinal and latitudinal directions. However, the parking lot is not parallel to the east–west direction, and traditional parameterization methods can result in many grids being located outside the parking lot. The absence of surface wave paths across these grids can lead to reduced stability of the inversion system alongside increased computational costs. To address this issue, we utilize coordinate rotation to approximately rotate the survey lines in an east–west orientation. The true coordinates are restored once the inversion process is complete.

Among the 36 wells at the site, three contain S-wave velocity logging results (Fig. 13(a)). We utilized these results as the basis for constructing an initial layered model. The model comprises layers of increasing thickness, ranging from shallow to deep. Specifically, it consists of 15 layers with a thickness of 3 m, two layers with a thickness of 5 m, and three layers with a thickness of 10 m. The horizontal spacing of the inversion grid is set to be comparable to the station spacing, representing the physical resolution limit. Specifically, the horizontal spacing is 0.00005° (approximately 5.5 m) in both the longitudinal and latitudinal directions.

During tomography, the density and compressional wave (P-wave) velocity are both determined from the S-wave velocity via empirical relationships [77]. The same empirical relationships were used to derive the density and P-wave velocity of the initial layered model, and its depth sensitivity kernel for the fundamental mode Rayleigh-wave phase velocity is shown in Fig. 13(b). This clearly shows that surface waves become more sensitive to deeper layers with decreasing frequency. The frequency range covered by our dispersion data is 2 to 15 Hz and exhibits excellent sensitivity to media shallower than 60 m, whereas its sensitivity to deeper media is relatively weaker. Therefore, we restrict the interpretation of the inversion results to depths shallower than 60 m.

4.2. Checkerboard tests

To evaluate how well our path coverage can retrieve the subsurface structure, we set up two 3D checkerboard models. The first model has a horizontal anomaly size of approximately

14 m \times 17

m (Fig. 14), whereas the second model has a horizontal anomaly size of approximately

23 m \times 28

m (Fig. 15). The velocity of the checkerboard model is determined by applying a perturbation of

\pm 10 %

to the initial model used in the real-data inversion. We calculated travel time data on the basis of the theoretical model, which had the same station pairs and frequency coverage as the actual data. We added 1% Gaussian random noise to the simulated travel time data to make them more realistic. We used the same weight factors and number of iterations as in the actual inversion. The recovery results show that the

14 m \times 17

m anomaly can be clearly resolved at a depth of 30 m, and the maximum resolvable depth of the

23 m \times 28

m anomaly is approximately 55 m, which is analogous to the results obtained from sensitivity kernel analysis.

Importantly, the performance of the checkerboard test can be influenced by various factors, including noise levels, model parameterizations, and anomaly distributions. Hence, blindly relying on the test results is not recommended. A comprehensive analysis and interpretation of the actual inversion results should be performed by integrating other relevant information. The results obtained from the checkerboard test can offer a preliminary assessment of the resolution, mitigating the potential for overinterpretation of the inversion results.

4.3. S-wave velocity model

The objective function of tomography includes data residuals and regularization terms for model roughness, so it is necessary to determine an appropriate weight factor to balance their contributions to the objective function. We tested multiple weight factors and compared the data fit and model smoothness of the inversion results. Finally, a weight factor of one was determined. The average value of the surface wave traveltime residuals decreased from −0.0165 s to a final value of −0.0006 s after 10 iterations and rarely decreased further. The distribution of residuals exhibited a greater concentration at approximately 0 s (Fig. 16). The results demonstrate that the inverted model effectively captures and fits the dispersion measurement data and that the inversion has successfully converged.

The drilling results indicate the presence of karst caves at depths of 40 m or greater. Fig. 17 depicts the horizontal slices of the inverted 3D model at depths of 45, 50, and 55 m. Extensive low-velocity anomalies, which are consistent with the primary locations identified by the drilling results, are evident in the southwest and northeast regions of the site. Conversely, the southeastern and northwestern regions show high-velocity anomalies, which are indicative of shallower bedrock surfaces. On the basis of the results of sensitivity analysis and resolution testing, our model encounters challenges in precisely identifying anomalies that are smaller than 10 m in both the horizontal and vertical dimensions, especially at depths greater than 40 m. Consequently, it is unrealistic to expect our velocity model to completely uncover all karst cave locations. By superimposing the locations of the caves identified by drilling data onto the velocity slices, the majority of cave positions are associated with low-velocity anomalies, with the exception of a few small caves that do not show any velocity anomalies. A construction team previously attempted grouting treatment of the voids in the abandoned parking lot to construct a new office building. However, the treatment failed to meet building safety standards. The construction plan in the area has been at a standstill since its initiation in 2019, and as of 2023, construction has not yet commenced. The injection of cement slurry underground has had a minimal effect on improving the elastic strength of the underground medium. We speculate that this could be due to the presence of two widely distributed low-velocity zones in the velocity model, which may correspond to cavities and/or highly porous gravel. These observations suggest that our inversion results effectively capture the comprehensive distribution of underground karst features.

Fig. 18 depicts the two velocity profiles AA′ (Fig. 18(a)) and BB′ (Fig. 18(b)), whose locations in the 3D model are marked in Fig. 17. The lithological drilling results that align with the areas intersected by the profiles are overlaid. The analysis of the drilling data suggests that the underground medium beneath the parking lot can be classified into two distinct layers: sedimentary deposits overlaying bedrock. The S-wave velocity gradually increases from 100 to 400 m·s⁻¹ within the depth range of 0 to 50 m. Beyond this depth range, a significant increase in velocity is observed, which is attributed to the presence of bedrock. The locations of the velocity interface depressions closely align with the positions of the karst caves identified by the drilling data. The reasons for the absence of a few karst caves in the velocity profiles will be addressed in the subsequent discussion.

5. Discussion

5.1. Limiting factors of the model resolution

A significant aspect of inversion resolution pertains to the extent to which velocity anomalies of various sizes can be recovered. To investigate the impact of station density on resolution, an analysis was conducted to compare the inversion results when half of the stations were uniformly eliminated from the dataset. The results revealed a noteworthy decline in the recovery effectiveness of the checkerboard model (Fig. S1 in Appendix A). The identification of subsurface low-velocity anomalies was less discernible when relying on a sparse array (Fig. S2 in Appendix A). On the other hand, the complex structure of real-world media can also make it difficult for us to achieve ideal resolution. For example, in Fig. 18(b), well #18 indicates that the bedrock depth is approximately 40 m, whereas our velocity model shows that the steep velocity interface at that location is at 50 m. The vertical resolution is not strongly influenced by station density, whereas it is predominantly determined by the sensitivity of surface wave dispersion inversion. Surface wave waveform inversion has the potential to improve resolution [77], [78]. Additionally, incorporating higher-mode surface wave data can also increase the vertical resolution of the inversion [49], [79], [80], [81], [82]. This paper does not discuss the use of new inversion methods and new types of data, such as multimode Rayleigh and Love waves, as well as the spectral ratio of the horizontal to vertical components of ambient noise, to improve tomography resolution. Future research will explore these avenues further.

Another question arises when the inversion resolution is considered: How accurate is it to interpret the anomaly distribution using the velocity model? We have 36 lithological boreholes on site, three of which also include S-wave velocity logging results. We seek to examine the relationship between

V_{s}

and lithology by analyzing data from three boreholes. As depicted in Fig. 19, with the exception of gravel and limestone, which exhibit distinctive high-velocity characteristics, and karst caves, which exhibit low-velocity characteristics, there are no other notable velocity differences among various lithologies. The velocity of different soils fluctuates, and these discrepancies in velocity are also nonconstant. Considering the complexity of the inversion resolution mentioned above, we limit the interpretation of the results in this study to the spatial distribution of low-velocity anomalies and the depth variations in the bedrock, to avoid overinterpretation.

5.2. Advantages of adding asynchronous observations

This study presents a new synchronous–asynchronous observation system that enables dense array observations of the research area with a limited number of stations. The ray path density of surface waves at frequencies of 2.0 and 5.4 Hz in this study is depicted in Fig. 20. To simulate traditional synchronous observations, half of the stations were purposefully reduced. The addition of asynchronous observations results in a noticeable increase in the ray path density of surface waves, leading to improved coverage for the study area. We successfully achieved a significant increase in the ray path density with only two observation rounds in this study. Theoretically, the ray path density can be further increased by adding more rounds of observations. The new observation system significantly reduces the number of stations necessary for dense-array observations, thereby resulting in a substantial cost reduction. Such a reduction is almost unattainable through traditional synchronous observations. This synchronous–asynchronous observation system has the potential to serve as an alternative option for future studies in dense-array ambient noise imaging, particularly when the number of stations is restricted.

6. Conclusions

This study demonstrates the application of short-term synchronous–asynchronous ambient noise tomography in urban karst investigations. By utilizing a synchronous–asynchronous ambient noise observation system, we are able to achieve extensive spatial coverage for the study area with a limited number of stations. The similarity found among different components of

C^{1}

functions is used as the stacking weight, and we successfully mitigate the influence of a nonuniform source distribution. The

C^{2}

method effectively retrieves the EGFs between asynchronous stations when the high modes in

C^{1}

functions are filtered by a multicomponent Rayleigh wave stacking technique. Through ambient noise surface wave tomography, we obtain a 3D S-wave velocity structure beneath the study area up to a depth of 70 m. The results reveal significant low-velocity anomalies below 40 m in the northeastern and southwestern regions, which could be attributed to karst caves and/or fragmented gravel with high porosity, as supported by well information. The incorporation of

C^{1}

and

C^{2}

data results in a significant enhancement in spatial coverage and resolution, surpassing the utilization of only

C^{1}

functions. We conclude that the newly proposed synchronous–asynchronous ambient noise tomography enables rapid and high-resolution seismic imaging, even in complex ambient noise environments with limited seismometer resources. This method has the potential to provide valuable support for urban geological hazard prevention and control.

CRediT authorship contribution statement

Ya Liu: Writing – review & editing, Writing – original draft, Software, Methodology, Formal analysis, Conceptualization. Jianghai Xia: Writing – review & editing, Supervision, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. Bo Guan: Writing – review & editing, Software, Investigation, Formal analysis. Chaoqiang Xi: Writing – review & editing, Investigation, Formal analysis. Ling Ning: Writing – review & editing, Investigation, Formal analysis. Hao Zhang: Writing – review & editing, Investigation, Formal analysis.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study is supported by the National Natural Science Foundation of China (41830103) and the Project of Nanjing Center of China Geological Survey (DD20190281). We appreciate Prof. Hongjian Fang for sharing the DSurfTomo tomography package (https://github.com/HongjianFang/DSurfTomo). We thank the members of the Subsurface Imaging and Sensing Laboratory (SISL), Institute of Geophysics and Geomatics, China University of Geosciences for their help in data collection. We would like to thank Managing Editor Qizhuang Xing and the three anonymous reviewers for their valuable feedback, which has significantly improved the quality of this manuscript.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2025.02.001.

References

Publishing order | Descend order by publishing year | Descend order by cited within

[1]	Weaver RL, Lobkis OI.Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies.Phys Rev Lett 2001; 87(13):134301.

[2]	Campillo M, Paul A.Long-range correlations in the diffuse seismic coda.Science 2003; 299(5606):547-549.

[3]	Snieder R.Extracting the Green’s function from the correlation of coda waves: a derivation based on stationary phase.Phys Rev E Stat Nonlin Soft Matter Phys 2004; 69(4):046610.

[4]	Sabra KG, Gerstoft P, Roux P, Kuperman WA, Fehler MC.Extracting time-domain Green’s function estimates from ambient seismic noise.Geophys Res Lett 2005; 32(3):L03310.

[5]	Shapiro NM, Campillo M, Stehly L, Ritzwoller MH.High-resolution surface-wave tomography from ambient seismic noise.Science 2005; 307(5715):1615-1618.

[6]	Wapenaar K, Draganov D, Snieder R, Campman X, Verdel A.Tutorial on seismic interferometry: part 1—basic principles and applications.Geophysics 2010; 75(5):75A195-75A209.

[7]	Yao H, van RD der Hilst, de MV Hoop.Surface-wave array tomography in southeast Xizang from ambient seismic noise and two-station analysis—I.Phase velocity maps. Geophys J Int 2006; 166(2):732-744.

[8]	Yang Y, Ritzwoller MH, Levshin AL, Shapiro NM.Ambient noise Rayleigh wave tomography across Europe.Geophys J Int 2007; 168(1):259-274.

[9]	Rawlinson N, Salmon M, Kennett BLN.Transportable seismic array tomography in southeast Australia: illuminating the transition from Proterozoic to Phanerozoic lithosphere.Lithos 2014; 189:65-76.

[10]	Bao X, Song X, Li J.High-resolution lithospheric structure beneath China’s mainland from ambient noise and earthquake surface-wave tomography.Earth Planet Sci Lett 2015; 417:132-141.

[11]	Luo S, Yao H, Li Q, Wang W, Wan K, Meng Y, et al.High-resolution 3D crustal S-wave velocity structure of the Middle-Lower Yangtze River metallogenic belt and implications for its deep geodynamic setting.Sci China Earth Sci 2019; 62(9):1361-1378.

[12]	Nimiya H, Ikeda T, Tsuji T.Three-dimensional S wave velocity structure of central Japan estimated by surface-wave tomography using ambient noise.J GeophysRes Solid Earth 2020; 125:e2019JB019043.

[13]	Gaite B, Iglesias A, Villaseñor A, Herraiz M, Pacheco JF.Crustal structure of Mexico and surrounding regions from seismic ambient noise tomography.Geophys J Int 2012; 188(3):1413-1424.

[14]	Haned A, Stutzmann E, Schimmel M, Kiselev S, Davaille A, Yelles-Chaouche A.Global tomography using seismic hum.Geophys J Int 2016; 204(2):1222-1236.

[15]	Li C, Yao H, Fang H, Huang X, Wan K, Zhang H, et al.3D near-surface shear-wave velocity structure from ambient-noise tomography and borehole data in the Hefei urban area, China.Seismol Res Lett 2016; 87(4):882-892.

[16]	Sun R, Li J, Yan Y, Liu H, Bai L, Chen Y.Three-dimensional urban subsurface space tomography with dense ambient noise seismic array.Surv Geophys 2024; 45(3):819-843.

[17]	Nakata N.Near-surface S-wave velocities estimated from traffic-induced Love waves using seismic interferometry with double beamforming.Interpretation 2016; 4(4):SQ23-SQ31.

[18]	Xia S, Zhang C, Cao J.Ambient noise tomography for coral islands.Engineering 2023;25 Jun:182–93.

[19]	Carvalho J, Silveira G, Dumont S, Ramalho R.3D-ambient noise surface wave tomography of Fogo volcano, Cape Verde.J Volcanol Geotherm Res 2022; 432:107702.

[20]	Lehujeur M, Vergne J, Schmittbuhl J, Zigone D, LeChenadec A.EstOF Team. Reservoir imaging using ambient noise correlation from a dense seismic network.J Geophys Res Solid Earth 2018; 123(8):6671-6686.

[21]	Cheng F, Xia J, Ajo-Franklin JB, Behm M, Zhou C, Dai T, et al.High-resolution ambient noise imaging of geothermal reservoir using 3C dense seismic nodal array and ultra-short observation.J Geophys Res Solid Earth 2021; 126(8):e2021JB021827.

[22]	Zhou C, Xia J, Pang J, Cheng F, Chen X, Xi C, et al.Near-surface geothermal reservoir imaging based on the customized dense seismic network.Surv Geophys 2021; 42(3):673-697.

[23]	Gu N, Wang K, Gao J, Ding N, Yao H, Zhang H.Shallow crustal structure of the Tanlu fault zone near Chao Lake in eastern China by direct surface wave tomography from local dense array ambient noise analysis.Pure Appl Geophys 2019; 176(3):1193-1206.

[24]	Guan B, Mi B, Zhang H, Liu Y, Xi C, Zhou C.Selection of noise sources and short-time passive surface wave imaging—a case study on fault investigation.J Appl Geophys 2021; 194:104437.

[25]	Ning L, Dai T, Liu Y, Xi C, Zhang H, Zhou C.Application of multichannel analysis of passive surface waves method for fault investigation.J Appl Geophys 2021; 192:104382.

[26]	Park CB, Miller RD.Roadside passive multichannel analysis of surface waves (MASW).J Eng Environ Geophys 2008; 13(1):1-11.

[27]	Nakata N, Snieder R, Tsuji T, Larner K, Matsuoka T.Shear wave imaging from traffic noise using seismic interferometry by cross-coherence.Geophysics 2011; 76(6):SA97-SA106.

[28]	Cheng F, Xia J, Xu Z, Hu Y, Mi B.Frequency–wavenumber (FK)-based data selection in high-frequency passive surface wave survey.Surv Geophys 2018; 39(4):661-682.

[29]	Liu Y, Xia J, Xi C, Dai T, Ning L.Improving the retrieval of high-frequency surface waves from ambient noise through multichannel-coherency-weighted stack.Geophys J Int 2021; 227(2):776-785.

[30]	Mi B, Xia J, Tian G, Shi Z, Xing H, Chang X, et al.Near-surface imaging from traffic-induced surface waves with dense linear arrays: an application in the urban area of Hangzhou, China.Geophysics 2022; 87(2):B145-B158.

[31]	Cheng F, Xia J, Shen C, Hu Y, Xu Z, Mi B.Imposing active sources during high-frequency passive surface-wave measurement.Engineering 2018; 4(5):685-693.

[32]	Waltham T, Bell FG, Culshaw MG, Knez M, Slabe T.Sinkholes and subsidence: karst and cavernous rocks in engineering and construction.Springer, Berlin (2005)

[33]	Weary D.The cost of karst subsidence and sinkhole collapse in the United States compared with other natural hazards.In: Proceedings of the Fourteenth Multidisciplinary Conference; 2015 Oct 5–6; Rochester, MN, USA. Tampa: University of South Florida Tampa Library; 2015. p. 433–46.

[34]	Zhang R, Ai T, Ren L, Li G.Failure characterization of three typical coal-bearing formation rocks using acoustic emission monitoring and X-ray computed tomography techniques.Rock Mech Rock Eng 2019; 52(6):1945-1958.

[35]	Dong L, Tong X, Ma J.Quantitative investigation of tomographic effects in abnormal regions of complex structures.Engineering 2021; 7(7):1011-1022.

[36]	Dong L, Pei Z, Xie X, Zhang Y, Yan X.Early identification of abnormal regions in rock-mass using traveltime tomography.Engineering 2023; 22(3):191-200.

[37]	Hiltunen DR, Cramer BJ.Application of seismic refraction tomography in karst terrane.J Geotech Geoenviron Eng 2008; 134(7):938-948.

[38]	Aati AHA, Shabaan SH.Detection of karstic limestone bedrock by shallow seismic refraction in an area west of Assiut, middle Egypt.Lead Edge 2013; 32(3):316-322.

[39]	Peng D, Cheng F, Liu J, Zong Y, Yu M, Hu G, et al.Joint tomography of multi-cross-hole and borehole-to-surface seismic data for karst detection.J Appl Geophys 2021; 184:104252.

[40]	Zhang J, Liu S, Chen Q, Wang B, Ren C.Application of cross-borehole integrated geophysical methods for the detailed investigation of karst in urban metro construction.JEEG 2019; 24(4):525-536.

[41]	Shangxin F, Yufei Z, Yujie W, Shanyong W, Ruilang C.A comprehensive approach to karst identification and groutability evaluation—a case study of the Dehou reservoir, SW China.Eng Geol 2020; 269:105529.

[42]	Gan F, Han K, Lan F, Chen Y, Zhang W.Multi-geophysical approaches to detect karst channels underground—a case study in Mengzi of Yunnan Province, China.J Appl Geophys 2017; 136:91-98.

[43]	Foudili D, Bouzid A, Berguig MC, Bougchiche SS, Abtout A, Guemache MA.Investigating karst collapse geohazards using magnetotellurics: a case study of M’rara basin, Algerian Sahara.J Appl Geophys 2019; 160:144-156.

[44]	Kiernan M, Jackson D, Montgomery J, Anderson JB, McDonald BW, Davis KC.Characterization of a karst site using electrical resistivity tomography and seismic full waveform inversion.JEEG 2021; 26(1):1-11.

[45]	Xia J, Chen C, Li P, Lewis M.Delineation of a collapse feature in a noisy environment using a multichannel surface wave technique.Geotechnique 2004; 54(1):17-27.

[46]	Parker EH Jr, Hawman RB.Multi-channel analysis of surface waves (MASW) in karst terrain, southwest Georgia: implications for detecting anomalous features and fracture zones.J Environ Eng Geophys 2012; 17(3):129-150.

[47]	Chang JP, de SAL Ridder, Biondi BL.High-frequency Rayleigh-wave tomography using traffic noise from Long Beach, California.Geophysics 2016; 81(2):B43-B53.

[48]	Wang Y, Lin FC, Schmandt B, Farrell J.Ambient noise tomography across Mount St. Helens using a dense seismic array.J Geophys Res Solid Earth 2017; 122(6):4492-4508.

[49]	Fu L, Pan L, Li Z, Dong S, Ma Q, Chen X.Improved high-resolution 3D vs model of Long Beach, CA: inversion of multimodal dispersion curves from ambient noise of a dense array.Geophys Res Lett 2022; 49:e2021GL097619.

[50]	Stehly L, Campillo M, Froment B, Weaver RL.Reconstructing Green’s function by correlation of the coda of the correlation (C³) of ambient seismic noise.J Geophys Res 2008; 113(B11):B11306.

[51]	Curtis A, Halliday D.Source-receiver wave field interferometry.Phys Rev E Stat Nonlin Soft Matter Phys 2010; 81(4):046601.

[52]	Froment B, Campillo M, Roux P.Reconstructing the Green’s function through iteration of correlations.C R Geosci 2011; 343(8–9):623-632.

[53]	Spica Z, Perton M, Cal Mò, Legrand D, Córdoba-Montiel F, Iglesias A.3-D shear wave velocity model of Mexico and south US: bridging seismic networks with ambient noise cross-correlations (C¹) and correlation of coda of correlations (C³).Geophys J Int 2016; 206(3):1795-1813.

[54]	Ansaripour M, Rezapour M, Saygin E.Shear wave velocity structure of Iranian plateau: using combination of ambient noise cross-correlations (C¹) and correlation of coda of correlations (C³).Geophys J Int 2019; 218(3):1919-1938.

[55]	Chen Y, Saygin E.Empirical Green’s function retrieval using ambient noise source-receiver interferometry.J Geophys Res Solid Earth 2020; 125(2):e2019JB018261.

[56]	Zhang S, Feng L, Ritzwoller MH.Three-station interferometry and tomography: coda versus direct waves.Geophys J Int 2020; 221(1):521-541.

[57]	Rao H, Luo Y, Zhao K, Yang Y.Extracting surface wave dispersion curves from asynchronous seismic stations: method and application.Geophys J Int 2021; 226(2):1148-1158.

[58]	Chen Y, Saygin E, Kennett B, Qashqai MT, Hauser J, Lumley D, et al.Next-generation seismic model of the Australian crust from synchronous and asynchronous ambient noise imaging.Nat Commun 2023; 14(1):1192.

[59]	Boschi L, Weemstra C.Stationary-phase integrals in the cross correlation of ambient noise.Rev Geophys 2015; 53(2):411-451.

[60]	Liu Y, Xia J, Cheng F, Xi C, Shen C, Zhou C.Pseudo-linear-array analysis of passive surface waves based on beamforming.Geophys J Int 2020; 221(1):640-650.

[61]	Cheng F, Xia J, Luo Y, Xu Z, Wang L, Shen C, et al.Multichannel analysis of passive surface waves based on crosscorrelations.Geophysics 2016; 81(5):EN57-EN66.

[62]	Zhang H, Mi B, Xi C, Liu Y, Guan B, Ning L.Extraction of surface-wave phase velocities from ambient noise in the presence of local noise sources based on matched-field processing.J Appl Geophys 2022; 204:104755.

[63]	Ning L, Xia J, Dai T, Liu Y, Zhang H, Xi C.High-frequency surface-wave imaging from traffic-induced noise by selecting in-line sources.Surv Geophys 2022; 43(6):1873-1899.

[64]	Liu Y, Xia J, Xi C, Zhang H, Guan B, Dai T, et al.Enhancing noise sources in stationary-phase zones for accurate phase-velocity estimation of high-frequency surface waves.Geophysics 2023; 88(1):L1-L9.

[65]	Nayak A, Thurber CH.Using multicomponent ambient seismic noise cross-correlations to identify higher mode Rayleigh waves and improve dispersion measurements.Geophys J Int 2020; 222(3):1590-1605.

[66]	Lin FC, Moschetti MP, Ritzwoller MH.Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps.Geophys J Int 2008; 173(1):281-298.

[67]	Gribler G, Mikesell TD.Methods to isolate retrograde and prograde Rayleigh-wave signals.Geophys J Int 2019; 219(2):975-994.

[68]	Herrmann RB.Computer programs in seismology: an evolving tool for instruction and research.Seismol Res Lett 2013; 84(6):1081-1088.

[69]	Bensen GD, Ritzwoller MH, Barmin MP, Levshin AL, Lin F, Moschetti MP, et al.Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements.Geophys J Int 2007; 169(3):1239-1260.

[70]	Park CB, Miller RD, Xia J.Imaging dispersion curves of surface waves on multi-channel record.In: SEG technical program expanded abstracts 1998. Houston: Society of Exploration Geophysicists; 1998. p. 1377–80.

[71]	Gribler G, Liberty LM, Mikesell TD, Michaels P.Isolating retrograde and prograde Rayleigh-wave modes using a polarity mute.Geophysics 2016; 81(5):V379-V385.

[72]	Schimmel M, Paulssen H.Noise reduction and detection of weak, coherent signals through phase-weighted stacks.Geophys J Int 1997; 130(2):497-505.

[73]	Yao H, Gou Pédard, Collins JA, McGuire JJ, van RD der Hilst.Structure of young East Pacific rise lithosphere from ambient noise correlation analysis of fundamental- and higher-mode Scholte–Rayleigh waves.C R Geosci 2011; 343(8–9):571-583.

[74]	Luo Y, Yang Y, Xu Y, Xu H, Zhao K, Wang K.On the limitations of interstation distances in ambient noise tomography.Geophys J Int 2015; 201(2):652-661.

[75]	Fang H, Yao H, Zhang H, Huang YC, van RD der Hilst.Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application.Geophys J Int 2015; 201(3):1251-1263.

[76]	Luo Y, Xia J, Miller RD, Xu Y, Liu J, Liu Q.Rayleigh-wave mode separation by high-resolution linear Radon transform.Geophys J Int 2009; 179(1):254-264.

[77]	Brocher TM.Empirical relations between elastic wavespeeds and density in the earth’s crust.Bull Seismol Soc Amer 2005; 95(6):2081-2092.

[78]	Li J, Feng Z, Schuster G.Wave-equation dispersion inversion.Geophys J Int 2017; 208(3):1567-1578.

[79]	Pan Y, Gao L, Bohlen T.Random-objective waveform inversion of 3D–9C shallow-seismic field data.J Geophys Res Solid Earth 2021; 126:e2021JB022036.

[80]	Xia J, Miller RD, Park CB, Tian G.Inversion of high frequency surface waves with fundamental and higher modes.J Appl Geophys 2003; 52(1):45-57.

[81]	Pan L, Chen X, Wang J, Yang Z, Zhang D.Sensitivity analysis of dispersion curves of Rayleigh waves with fundamental and higher modes.Geophys J Int 2019; 216(2):1276-1303.

[82]	Liu Q, Chen X, Gao L, Yu Z, Chen J.Direct image dissimilarity inversion of ambient noise multimodal dispersion spectrograms.Bull Seismol Soc Amer 2023; 113(5):1960-1981.

RIGHTS & PERMISSIONS

THE AUTHOR

PDF (6337KB)

Part of a collection:

Supplementary files

Appendix A. Supplementary data

4762

Accesses

Citation

Detail

Sections

Recommended

Browse

Authors & reviewers

About the journal

Abstract

Graphical abstract

Keywords

Cite this article

1. Introduction

2. Data

2.1. Design of synchronous and asynchronous seismic arrays

2.2. Ambient noise characteristics

3. Methods

3.1. Calculation of C1

3.2. Calculation of C2

3.3. Dispersion measurement

3.4. Surface wave tomography

3.5. Workflow summary

4. Results

4.1. Initial model and depth sensitivity

4.2. Checkerboard tests

4.3. S-wave velocity model

5. Discussion

5.1. Limiting factors of the model resolution

5.2. Advantages of adding asynchronous observations

6. Conclusions

CRediT authorship contribution statement

Declaration of competing interest

Acknowledgments

Appendix A. Supplementary data

References

RIGHTS & PERMISSIONS

3.1. Calculation of C¹

3.2. Calculation of C²