A phase unwrapping approach in measuring surface wave phase velocities from ambient noise

SUMMARY

In the past two decades or so, ambient noise tomography (ANT) has emerged as a powerful tool for investigating high-resolution crustal and upper-mantle structures. A crucial step in the ANT involves extracting phase velocities from cross-correlation functions (CCFs). However, obtaining precise phase velocities can be a formidable challenge, particularly when significant lateral velocity variations exist in shallow subsurface imaging that relies on short-period surface waves from ambient noise. To address this challenge, we propose an unwrapping correction method that enables the accurate extraction of short-period dispersion curves. Our method relies on the examination of the continuity of phase velocities extracted from CCFs between a common station and other neighbouring stations along a linear array. We demonstrate the effectiveness of our approach by applying our method to both synthetic and field data. Both applications suggest our unwrapping correction method can identify and correct unwrapping errors in phase velocity measurements, ensuring the extraction of accurate and reliable dispersion curves at short periods from ambient noise, which is essential for subsequent inversion for subsurface structures.

1. INTRODUCTION

Ambient noise tomography (ANT) has become a powerful tool in imaging Earth's interior structures. ANT has gained widespread popularities in imaging the crustal and upper mantle structures in both regional (e.g. Ai et al. 2019a,b; Gu et al. 2019, 2022; Jiang & Denolle 2022; Shapiro et al. 2004; Lin et al. 2008; Yang et al. 2020; Wang et al. 2021; Zhang et al. 2022) and continental (e.g. Yang et al. 2007; Bensen et al. 2008; Shen & Ritzwoller 2016; Zhao et al. 2020; Han et al. 2021) scales. In this method, a surface wave dispersion curve is measured from a cross-correlation function (CCF) obtained by cross-correlating ambient seismic data recorded at a pair of stations. Subsequently, interstation surface wave dispersion curves from various station pairs within a seismic array are employed in tomography to constrain subsurface structures.

In most regional- and continental-scale ANT studies, surface wave dispersion curves extracted from CCFs typically span periods from a few seconds to several tens of seconds. These dispersion curves serve as valuable data for studying crustal and upper-mantle velocity structures (e.g. Lin et al. 2008; Luo et al. 2013; Shen & Ritzwoller 2016; Zhao et al. 2017; Li et al. 2018; Ai et al. 2019a). For these intermediate- and long-period ranges, global surface wave phase velocity models like GDM52 (Ekström 2011) with the shortest period of ∼25 s can be utilized as reference dispersion curves in unwrapping phases during measuring phase velocities from CCFs. Such a global reference model assists in the search for accurate short-period dispersion curves starting from the long-period ones as the dispersion curves of surface waves varies smoothly with periods.

However, when utilizing stations with noise data recorded only at short-period stations with cut-off periods of a few seconds, resulting CCFs from these stations lack intermediate- and long-period signals. In such a case, the previously mentioned intermediate-/long-period regional phase velocity reference model is not suitable for phase unwrapping. On the other hand, local-scale short-period phase velocity reference models are usually unavailable, making it challenging to accurately unwrap phase velocity dispersion curves.

To overcome this challenge, common approaches utilize array techniques to extract average phase velocity dispersion curves over a study region or subregions. For this purpose, various methods have been proposed, including the multichannel analysis of surface waves method (Park et al. 1999; Xia et al. 2000), the phase shift method (PSM, Park et al. 1998), the image transformation technique (EGFAnalysisTimeFreq) proposed by Yao et al. (2006, 2011), the high-resolution linear Radon transformation method developed by Luo et al. (2008) and the Frequency–Bessel transform (F-J) method (Wang et al. 2019; Hu et al. 2020). While these methods effectively address issues related to inaccurate phase velocity dispersion curve unwrapping, it is worth noting that these methods often rely on the assumption of laterally homogeneous velocity structures. Therefore, the phase velocity dispersion curves obtained through these methods typically represent an averaged or a highly smoothed model for a study region, lacking information about fine-scale velocity features.

Therefore, there is a demand for the development of a short-period phase velocity extraction or unwrapping technique capable of deriving short-period dispersion curves from CCFs in the absence of a short-period phase velocity reference model. This study introduces a short-period phase velocity unwrapping approach that does not require an accurate phase velocity reference model. Our method relies on an analysis of phase velocity continuities recorded at neighbouring stations along a linear array. We apply this technique to both synthetic and field data to illustrate its effectiveness in obtaining accurate short-period dispersion curves.

2. METHODS

2.1 Phase velocity measurements in ANT

One of the crucial steps in ANT is measuring phase velocities from CCFs, which is based on phase measurement of surface wave signals in a CCF. Phases (φ) of surface waves are converted into phase velocity (C) by unwrapping the cycles of propagation, as described by the following equation (e.g. Snieder 2004; Lin et al. 2008; Tsai 2009):

where r represents the station separation, ω denotes the angular frequency, φ denotes the phase of the surface wave signal in a CCF, the phase shift of π/4 arises from the interference of homogeneous sources distributed across a 2-D surface domain (e.g. Yao et al. 2006; Lin et al. 2008). In scenarios involving a 1-D linear domain, where sources are positioned exclusively along the lines extending outward from a station pair, this π/4 phase shift must be removed (e.g. Lin et al. 2008). This adjustment is applied in our 1-D simulation described in Section 3. And 2πN stands for the phase ambiguity that must be determined based on the cycle of propagation with reference to a phase velocity model.

As mentioned in ‘Introduction’ section, in cases involving broadband surface waves from ANT, N can be readily determined by comparing measured dispersion curves from CCFs with reference phase velocities at intermediate and long periods (> 25 s) according to GDM52 (Ekström 2011), or other established models. However, the challenge arises during short-period phase velocity measurements, especially at periods shorter than 5 s, where determining the appropriate N value is not so straightforward. In this study, our primary focus revolves around determining the suitable value of N, which plays a critical role in phase velocity measurements.

2.2 Unwrapping of surface wave propagation

To demonstrate our unwrapping method, we first utilize phase velocities over various station pairs along a linear array where the phase velocity is kept constant at 1 km s⁻¹. The station distribution is depicted in Fig. 1(a).

As illustrated in Fig. 1, taking station A as an example, we examine the phase velocities between station A and other stations by plotting phase velocity values at locations of individual stations along the linear array. In a scenario of accurate unwrapping with no cycle skipping, the phase velocities between station A and other stations are represented by the continuous dots shown in Fig. 1(b). The plus and minus signs positioned above or below the continuous dots denote phase velocities resulting from $\pm$1 cycle of wrong unwrapping.

By analysing the continuity of phase velocities between neighbouring stations relative to the common station A, we can determine if any inappropriate unwrapping occurs among the stations. In a case of wrong unwrapping, such as the examples shown in Fig. 1(c), where the phase velocities between stations A and B/C are wrongly unwrapped, phase velocity values plotted at stations B and C relative to other stations become discontinuous with apparent jumps appearing at the locations of stations B and C. To correct the phase velocities, we can simply unwrap the phase velocities by adding one to or subtracting one from the N value in eq. (1) and then recalculate the phase velocities, which results in the correct phase velocities for stations B and C as indicated by the arrows in Fig. 1(c).

The above illustration indicates that in cases where there is an issue of wrong unwrapping in phase velocity measurements, phase velocities plotted at a station with wrong unwrapping along a linear array exhibit significant phase velocity variation. By observing the discontinuities in the phase velocity values of neighbouring stations, we can assess the extent of wrong unwrapping in phase velocity measurements. This approach enables us to discern the overall trend of the phase velocity. If such issues are identified, correction measurements can be taken by adjusting the N value and subsequently recalculating the phase velocities.

Above, we demonstrate phase unwrapping correction with constant phase velocity along a linear array. However, real data applications typically involve significant variations in lateral seismic velocity. Based on this, we explore the application of phase unwrapping correction when there are differences in seismic velocity along a linear array.

Based on Fig. 1(b), we observe that, for example, taking station A as a reference, phase correction requires adjusting N and recalculating the phase velocities of other stations along the array according to their distances from station A. The closer a station is to station A, the greater the phase velocity difference caused by unwrapping errors, as depicted in Fig. 1(b). Therefore, when a reference phase velocity dispersion curve is available, the dispersion curve between station A and other station which is near station A can be accurately extracted without the need to perform our unwrapping correction.

After determining the phase velocities of stations close to station A, we then check the phase velocities of stations farther away in order of their distances from station A. Taking stations A, B and C as examples, here we correct the phase velocity between stations A and C based on the phase velocity between stations A and B, as shown in Fig. S1 (Supporting Information). Assuming that the phase velocity between stations A and B is accurate, we can judge whether the phase unwrapping between stations A and C is accurate. In any period T, the difference between the traveltime between stations A and C and the traveltime predicted using the phase velocity between stations A and B should be less than T/2, because when the traveltime difference exceeds T/2 and the phase velocity would be close to another unwrapping erroneous phase velocity point. Therefore, the premise for the accurate correction is:

where $t_{\mathrm{A}\mathrm{C}}$ is the traveltime between stations A and C, $d_{\mathrm{A}\mathrm{C}}$ is the distance between stations A and C, $V_{\mathrm{A}\mathrm{B}}$ is the average phase velocity between stations A and B and $T$ is the period. Formula (2) can be rewritten as:

In which the interstation distance between stations A and B shown in $d_{\mathrm{A}\mathrm{B}}$ must meet the requirement of

Since we use adjacent station pairs to perform corrections one by one, $d_{\mathrm{B}\mathrm{C}}$ can be considered as the distance between adjacent stations $\mathrm{d}x$⁠. If we use $V_{\mathrm{A}\mathrm{B}}$ to approximate the wavelength, then the formula can be approximated as:

Thus, if the phase velocity difference between $V_{\mathrm{A}\mathrm{B}}$ and $V_{\mathrm{B}\mathrm{C}}$ is not so large, we can perform unwrapping correction successfully. Fig. S1 (Supporting Information) shows the relationship between the minimum interstation distance and the velocity difference between adjacent stations and indicates the conditions under which $V_{\mathrm{B}\mathrm{C}}$ can be directly unwrapped by referencing $V_{\mathrm{A}\mathrm{B}}$⁠. As revealed by Fig. S1 (Supporting Information), when the difference between $V_{\mathrm{A}\mathrm{B}}$ and $V_{\mathrm{B}\mathrm{C}}$ is not so large (e.g. the ratio of $V_{\mathrm{B}\mathrm{C}}$ and $V_{\mathrm{A}\mathrm{B}}$ is between 0.8 and 1.3), the unwrapping correction can still proceed smoothly even when $\mathrm{d}x$ is several wavelengths.

3. APPLICATION TO SYNTHETIC DATA

3.1 Simulation of CCFs

To illustrate the effectiveness of our approach in correcting phase velocity measurements resulting from improper unwrapping, we first apply our approach to synthetic data from a linear seismic array.

Numerous studies of ANT have been conducted to image regions characterized by a combination of low-velocity basins and high-velocity terrains, where the cycle skipping issue of surface waves significantly affects the measurements of phase velocities. Within this context, we conduct a simulation for a model that exhibits substantial lateral velocity variations. This model is depicted in Fig. 2(a), featuring low velocities on the left-hand side and high velocities on the right-hand side. In this model, Vs ranges from 0.5 km s⁻¹ on the upper left to 3.5 km s⁻¹ on the lower right. This model emulates a terrain with a combination of a basin and a mountain, akin to a case study presented in Section 4 that employs seismic noise data collected from a linear array traversing Taiwan Island from west to east.

In this simulation model, a total of 51 receivers are emplaced along a linear seismic array with a station interval of 2 km. To synthesize seismic data, we employ a staggered-grid finite-difference method developed by Wang et al. (2012, 2015). We simulate seismic waveforms generated by 51 sources colocated with the 51 stations. For each source, the source function is the first derivative of a Gaussian function with a peak frequency of 0.2 Hz, and a vertical point source is emplaced at the location of each receiver on the surface of the 2-D model. The synthetic seismograms from one source located at the far left side are depicted in Fig. 2(b). These simulated waveform records display a distinct dispersion phenomenon of surface waves with the signals also containing body waves, as portrayed in Fig. 2(b).

After completing the waveform simulation, we collect all the vertical-component simulated waveform data from 51 sources recorded at each station, which serve as seismic noise data from the 51 noise sources. Then, following the noise data processing procedures of Bensen et al. (2007), we calculate CCFs for each station pair using the seismograms from each of the 51 noise sources and then stack the 51 CCFs between each station pair to obtain the final CCFs. As we only use vertical-component data in our calculation of CCFs, the resulting surface waves in CCFs are Rayleigh waves propagating between station pairs. We plot a record section of the stacked CCFs (Fig. 2c) between a reference station, located at the leftmost position of the seismic array, and other stations, which displays clear surface waves with a lower moveout velocity at the section of ∼0–30 km and a higher moveout velocity at the ∼30–100 km section, mimicking the Vs transition from the low-velocity left-hand side to the high-velocity right-hand side in Fig. 2(a).

3.2 Measurements of phase velocities

With the simulated CCFs, we adopt a frequency–time analysis method (FTAN, Levshin & Ritzwoller 2001) to measure phase velocities of surface waves from these CCFs. This FTAN method needs a reference phase velocity model to unwrap phases of surface waves in order to pick the right branch of phase velocity dispersion curves. For a local-scale region, accurate reference phase velocity models are typically not available. However, based on all the CCFs, we can readily obtain an average phase velocity dispersion curve using array-based methods as we discuss in ‘Introduction’ section. Here, we adopt the PSM of Park et al. (1998), which scans phase shifts of surface waves in CCFs along the profile to get an average phase velocity dispersion curve.

With the average phase velocity dispersion curve as the reference model, we measure dispersion curves of surface waves from all the CCFs using the FTAN method. In Fig. 3, we plot examples of the dispersion curves from CCFs for stations 01 and 20 with other stations. For comparison, we also calculate the theoretical dispersion curves from our simulation model using CPS package (Computer Programs in Seismology, Herrmann 2013), which are shown in Figs 3(b) and (d). In the phase velocity measurements, due to the strong lateral variation of Vs structures, the average phase velocity dispersion curve cannot properly represent the accurate interstation dispersion curves, resulting in a lot of wrongly picked dispersion curves, denoted as red or blue curves in Figs 3(a) and (c). Dispersion curves for other station pairs show similar patterns with a few dispersion curves picked with wrong branches.

To identify and correct these wrongly picked dispersion curves, we perform phase velocity unwrapping and correction procedures for each station and each period by selecting a series of separate periods along the dispersion curves. Here, we choose stations 01 and 20 for phase velocity measurement at 3 s period as examples to illustrate our procedures.

In the procedures, we first plot the variations of phase velocities at a period for surface waves between a common station and other stations along the linear array. The phase velocities are plotted as a function of station index ordered by their locations along the linear array from the left to the right as shown in Figs 4(a) and (c). If the phase velocities are all correctly measured, they would form a continuous line along the linear array. If phase velocities at some stations are wrongly measured due to improper unwrapping, the variations of phase velocities along the linear array exhibit strong discontinuities, as shown in Fig. 4(a) for station 01 and in Fig. 4(c) for station 20. These jumps in phase velocities can be readily identified with the guidance of the continuous phase velocity trends along the linear array.

With the identified phase velocity jumps, we then correct the wrongly picked measurements of phase velocity by adjusting the N value and subsequently recalculating the phase velocity to get the accurate dispersion curves. After the correction, phase velocities display as continuous dots as shown in Fig. 4(b) for station 01 and Fig. 4(d) for station 20, respectively. We repeat the above correction process for each station along the linear array, and finally calculate the phase velocity differences between measured dispersion curves and theoretical dispersion curves, as shown in Fig. S2 (Supporting Information), displaying that the dispersion curves after correction are more consistent with the theoretical dispersion curves.

3.3 Phase velocity tomography

To further evaluate the impact of unwrapping correction on tomography results, we conduct ANT using dispersion data before and after applying our unwrapping correction as well as utilizing predicted dispersion data from the simulation model.

We adopt the ray-theory-based ANT method developed by Zhao et al. (2020) to invert dispersion curves for phase velocities. Originally, this method was used to invert for 2-D phase velocity maps. In this method, the phase velocity at any geographic point is calculated by averaging the phase velocity coefficients from various neighbouring grid nodes using a 2-D Gaussian function. The characteristic length of the Gaussian weighting function acts as the smoothing parameter during inversion. Thus, the method first inverts for phase velocity coefficient at each grid node and then uses the 2-D Gaussian weighting function to calculate the phase velocity at any geographic point. This method can be readily modified to suit the inversion for 1-D phase velocity profiles along a linear array. In the case of a 1-D linear array, the grid consists of linear nodes, and the 2-D Gaussian smoothing function in a 2-D case is adapted into a 1-D Gaussian smoothing function, which is used to control the smoothness of phase velocity variations along the linear array. Similar to its application in a 2-D array, this method first inverts for the phase velocity coefficients along the linear nodes, and then calculates the phase velocity at any point along the linear array by averaging the surrounding phase velocity coefficients based on the 1-D Gaussian weighting function.

In the tomography, the size of grid node along the linear array is set to 2 km and regulation parameters such as smoothing and damping are all set to the same values. After getting phase velocity profiles along the linear array at individual periods using the inversion method described above, we then combine these phase velocity profiles at periods from 1 to 6 s to form phase velocity transects as shown in Figs 5(a)–(c). The differences between the phase velocities based on dispersion data before and after unwrapping correction, relative to the predicted phase velocities calculated from the simulation model are shown in Figs 5(d) and (e). In addition, we conduct tests to assess the impact of unwrapping correction on the phase velocity distribution when utilizing measured dispersion curves based on regionalized average reference models, as illustrated in Fig. S3 (Supporting Information). All the results show that the phase velocity transects based on dispersion data with the unwrapping correction are very close to the theoretical ones (Fig. 5e and Fig. S3e, Supporting Information). But the phase velocity transects based on dispersion data without the unwrapping correction are significantly different from the theoretical ones (Fig. 5d and Fig. S3d, Supporting Information). The above comparisons suggest our unwrapping correction procedures are effective in correcting wrongly picked dispersion curves due to wrong unwrapping and obtaining accurate dispersion curves from CCFs.

4. APPLICATION TO FIELD DATA IN TAIWAN

4.1 Noise data processing

To further demonstrate the effectiveness of our method in correcting unwrapping errors in phase velocity measurements, we choose noise data deployed in Central Taiwan as a case study. This region features a complex terrain with both basins and mountains, resulting in significant variations of surface wave phase velocities. We collect continuous noise data from a dense seismic array of the TAIGER (TAiwan Integrated GEodynamics Research) project. This array consists of 62 short-period seismometers with a cutoff period of 5 s, covering a ∼135 km profile with an average station spacing of ∼2 km (Fig. 6a). The noise data span a duration of three months from 2009 March to June.

The procedures for obtaining CCFs from noise data closely follow the method outlined by Bensen et al. (2007) as we employ for the simulation noise data. First, the continuous noise data are cut into a series of daily segments. The mean, trend, instrument responses are removed from the daily segments. Subsequently, a bandpass filter of 0.2–20 s period range is applied to all daily segments. Following this, the data undergo temporal normalization and spectral whitening processes to suppress effects from energetic sources, such as earthquakes. Lastly, daily CCFs are computed using these preprocessed data for each station pair, and they are then linearly stacked.

Fig. 6(b) displays examples of CCFs between station M01 and other stations. We can observe clear Rayleigh wave signals at positive time lags, suggesting noise sources come from the western coastline of Taiwan. Therefore, we employ the positive-component CCFs to measure interstation surface wave phase velocities. Higher-mode signals seem visible in the positive-component CCFs as outlined by the red lines in Fig. 6(b). In this study, to illustrate our unwrapping correction method, we only measure phase velocities from the fundamental-mode surface waves by isolating the fundamental-mode signals using the time windows of blue lines depicted in Fig. 6(b).

4.2 Phase velocity measurement

Following the same workflow employed to process synthetic data presented in Section 3, we first compute an average regional dispersion curve by employing PSM. Using the resulting average dispersion curve as the reference model, we obtain phase velocity dispersion curve from each CCF using the FTAN method. Figs 7(a) and (c) illustrate examples of the measured phase velocity dispersion curves between station M01/M20 and other stations.

Subsequently, we examine the phase velocity measurements at various individual periods by plotting the phase velocity values along the linear array centred at each individual station. Figs 8(a) and (c) showcase two examples of the phase velocity profiles at 3 s period from CCFs centred at stations M01 and M20. These figures reveal various incorrectly picked phase velocities as denoted as red or blue dots, corresponding to the wrong dispersion curves in red or blue lines shown in Figs 7(a) and (c). To correct the identified erroneous phase velocities, we apply the same unwrapping correction procedures as used for the synthetic data. The resulting phase velocity variations along the linear profile after the correction are shown in Figs 8(b) and (d), displaying continuous variation of phase velocities along the profile. We apply the correction method to all phase velocity measurements centred at each station. Examples of dispersion curves after the correction are plotted in Figs 7(b) and (d) for CCFs with the common station M01 and station M20, displaying tighter distribution of dispersion curves compared with the dispersion curves before the correction.

Then, we select those dispersion curves with the interstation distances longer than one wavelength (Luo et al. 2015) and the signal-to-noise ratios of surface waves in CCFs greater than 10 for the subsequent tomography as presented in the following section.

4.3 Phase velocity tomography

Like the case of synthetic data, we also perform ANT for profiles along the linear array in Central Taiwan using the two sets of dispersion curves before and after the unwrapping correction. We obtain phase velocity variations along the profile at periods from 1 to 6 s by employing the ANT method of Zhao et al. (2020).

Subsequently, we assemble the resulting phase velocities at periods between 1 and 6 s along the linear array to create phase velocity transects, as depicted in Figs 9(a) and (b). A clear distinction emerges when comparing two phase velocity transects. In mountainous regions in the eastern part of Taiwan, phase velocities obtained before unwrapping correction are 5–15 per cent higher than those after unwrapping correction. Near the basin-mountain junction area, phase velocities before unwrapping correction are higher than those after unwrapping correction at shorter periods, and lower at higher periods.

In contrast to synthetic data where the true model is known, assessing model accuracy in practical data applications is more challenging because we lack precise knowledge of the subsurface structures. Nevertheless, we can evaluate the resulting phase velocity model by examining the misfits in phase traveltimes of surface waves. To this end, we calculate misfits in phase traveltimes at several periods and plot the histograms of these misfits in Figs 10(a)–(d). The histograms reveal that phase velocity results based on uncorrected dispersion data exhibit significantly larger misfits with several separate misfit peaks in the histograms, indicating phase velocity errors caused by incorrect unwrapping. In contrast, misfits of phase traveltimes based on the corrected unwrapping data display Gaussian distributions with almost zero means and less than 0.8 s standard deviations. This comparison of misfits in phase traveltimes suggests that our unwrapping correction leads to more accurate phase velocity measurements, echoing the same conclusion reached with synthetic data.

5. DISCUSSIONS

5.1 Selections of reference phase velocity models

The main reason for incorrect unwrapping in phase velocity measurement is due to that the reference model is significantly different from the true one. If the reference model is sufficiently accurate, erroneous unwrapping can be avoided. Unfortunately, in practical applications, especially for a region with strong lateral variations of subsurface seismic velocities, it is challenging to have a reference model featuring the strong variations as we do not know the subsurface structures. As discussed before, it is quite straightforward to obtain an average phase velocity model of a study area, which can serve as the first-order reference model in picking the first-round dispersion curves. Then, our method can work as the second-order unwrapping correction to make sure the phase velocities measured from CCFs are accurate.

One might wonder if we can circumvent the unwrapping errors by dividing a linear profile into several subregions and deriving regional average phase velocity model for each subregion, which could serve as a more accurate model for phase velocity measurements. Using the synthetic data presented in Section 3, we explore this procedure by dividing the study area into three sub-regions from left to right: the basin area, the basin-mountain junction area and the mountainous area. We adopt the PSM to extract average dispersion curves for subregions. Then, based on the central point of each CCF path, we select the most suitable one among three reference models to calculate phase velocities for each CCF. However, this approach of subregionalization still results in a high incidence of incorrect phase velocity unwrapping. Figs S3(b) and (d) (Supporting Information) display the phase velocity transects and the phase velocity differences derived using the three reference models for dispersion data, demonstrating that the use of multiple average models fails to eliminate the issue of incorrect unwrapping.

In the above process, we compare our unwrapping results using a single average phase velocity model with the corresponding results using several regional average models. We discover that employing a single average model for the entire region can yield the same phase velocity transects after unwrapping correction as using several regional models. Therefore, in the practical application, using an average model for the entire study area as an initial reference model is adequate to obtain accurate phase velocity measurements after the unwrapping correction.

5.2 Applications to practical seismic arrays

In this study, we demonstrate the applicability of our unwrapping correction method to dense linear arrays, particularly in regions with complex subsurface structures.

In our method, we select a few individual periods along dispersion curves to construct phase velocity profiles along a linear array for correction purpose. A smaller interval when choosing periods for correction results in a greater number of correction periods, leading to more accurate dispersion data after correction. If there are still discontinuities in phase velocity distributions at the shorter or longer period ends, but the phase velocity distributions remain continuous and stable at most periods in between, we need to remove these discontinuities during the unwrapping correction process. Therefore, by employing multiperiod unwrapping correction, we can determine the actual period range of the dispersion curve associated with each station pair and obtain accurate dispersion curves. Additionally, in cases where the deviation of N exceeds 1 (eq. 1), we can repeat the correction process multiple times until we achieve a continuous phase velocity distribution between adjacent stations.

This unwrapping correction technique can also be applied to dense 2-D arrays, an increasingly popular choice for high-resolution subsurface imaging due to the cost-effectiveness of portable seismometers. When implementing our method on 2-D arrays, any seismic station may be designated as a ‘central station’. For each selected central station, we can generate a 2-D phase velocity variation surface based on the phase velocities between individual stations and the central station. Phase velocities unwrapped incorrectly can be identified by observing discontinuities in the 2-D phase velocity variation across individual stations, similar to the patterns observed in the 1-D linear array. Subsequently, a similar unwrapping correction can be applied. This process can be then repeated for all seismic stations designated as central stations, ensuring accurate unwrapping correction across the entire seismic array.

6. CONCLUSIONS

In this study, we introduce an unwrapping correction method designed to identify and rectify inaccuracies in short-period phase velocity measurements resulting from erroneous unwrapping. This method is based on the examination of phase velocity continuity along a linear array. By applying our method to both synthetic and field data, we demonstrate its effectiveness in correcting phase velocity measurements affected by unwrapping errors. Our method is proved to be particularly advantageous in extracting accurate phase velocities of surface waves for terrains with substantial variation of subsurface structures. Furthermore, we suggest that our method can also be applied to dense 2-D arrays by generating 2-D phase velocity variation surfaces for unwrapping correction, showcasing its versatility and potential applicability in various geological settings. Our unwrapping correction method holds promise for improving the quality of dispersion curves, facilitating more precise characterizations of shallow subsurface velocity structures based on short-period seismometers.

Acknowledgement

This work is supported by the National Science Foundation of China (NSFC, nos 42074055, 42374057, 42304067 and 42074069), Natural Science Foundation of Hubei Province of China (no. 2023AFD211), and supported by MOST Special Fund from the State Key Laboratory of Geological Processes and Mineral Resources and by the Fundamental Research Funds for the Central Universities, China University of Geosciences (MSFGPMR2022-4).

CONFLICT OF INTEREST

All authors declare that they have no conflicts of interest.

DATA AVAILABILITY

The seismic waveform data are downloaded from the IRIS Data Management Center at www.iris.edu (last accessed in 2023 June). All plots are made using Generic Mapping Tools version 6.3.0 (https://www.generic-mapping-tools.org/). The code and data used in this study will be shared upon reasonable request to the authors.

References