# Bulletin of the Seismological Society of America

## Abstract

In surface‐wave dispersion measurement, we use cross correlation to compare the recorded signals at different stations to measure their phase difference. However, the measurement resolution is poor at low frequencies and the poor resolution was often incorrectly interpreted as large measurement errors. Here, we propose a new nonlinear signal comparison (NLSC) approach to achieve a uniform high‐resolution measurement across a wide frequency band. Furthermore, we can control the overall resolution in NLSC by an adjustable parameter. The traditional cross correlation is a special case of NLSC. We use both synthetic and recorded seismic data to demonstrate the effectiveness of NLSC by extracting surface Rayleigh‐wave phase‐velocity dispersion curves for different datasets, including synthetic multimode Rayleigh‐wave data, synthetic global Rayleigh waves of a Martian seismic model, and a real ambient noise correlation data processed from the USArray.

## Introduction

Signal comparison is a basic operation in seismic data analyses and processing, such as seismic waveform inversion, seismic imaging and migration, and surface‐wave analysis. For example, in seismic waveform inversion (e.g., Tarantola, 1984) we need to compare the synthetic seismic traces with the observed seismic traces to update the model. In seismic migration and imaging, we need to compare the forward propagated and the backward propagated wavefield to form an image in the subsurface (e.g., Claerbout, 1985). In surface‐wave (or body‐wave) analyses, we often need to compare surface waves recorded by two stations and determine their relative phase shift due to propagation to obtain the dispersion curves to invert for the subsurface structure. However, in most cases, signal comparison is mathematically formulated as cross correlation. The cross correlation naturally resides in the framework of Fourier analysis and the least‐squares minimization (e.g., the Fourier coefficient is computed by cross correlating the signal with the sines and cosines), and linear stacking of signals. We refer to traditional cross correlation and linear stacking as linear signal comparison (LSC).

However, using cross correlation to resolve a small time shift between two signals of similar wavelets depends on the frequency of the signals, with poor resolving power at low frequencies. This is widely observed for measuring low‐frequency surface‐wave dispersion (e.g., Foster *et al.*, 2014) and leads to difficulty in identifying the dispersion curve at low frequencies. In fact, cross correlation is just one technique to compare signals and is by no means the only method. The purpose of this article is to propose a new method, based on a nonlinear signal comparison (NLSC), to compare two signals to overcome those drawbacks inherent in the cross correlation. We choose to illustrate our NLSC in the context of surface‐wave dispersion curve extraction.

The surface wave is an interface phenomenon and its propagation is sensitive to velocity structures close to the interface. Surface‐wave dispersion inversion is a powerful tool extensively used in both global seismology (Nolet, 1975; Shapiro *et al.*, 2005; Herrmann, 2013; Foster *et al.*, 2014) and exploration seismology (Mcmechan and Yedlin, 1981; Park *et al.*, 1999; Xia *et al.*, 2007) to interrogate subsurface structures. In borehole acoustic logging, the type of surface wave is called the Stoneley wave, propagating along the fluid‐filled wellbore and its dispersion is critical to infer formation elasticity properties around the borehole (e.g., Tang and Cheng, 2004).

In the inversion of surface‐wave data, the first step is to estimate the frequency‐dependent phase‐velocity dispersion from the seismic data. Depending on the configuration of sources and receivers, many dispersion measurement techniques have been proposed. The two‐station method (e.g., Knopoff *et al.*, 1966; Yao *et al.*, 2006) has been widely used in seismology to extract dispersion curves; however, this method suffers the low‐resolution problem mentioned earlier and Yao *et al.* (2006) illustrated this nicely in their figure 14. For global seismology, if there is just one seismometer, the one‐station method was used to measure the globally averaged dispersion for surface‐wave tomography (e.g., Sato, 1957; Ekstrom *et al.*, 1997; Ritzwoller and Levshin, 1998). More recently, Zheng *et al.* (2015) proposed to use the same single‐station method but with global Rayleigh‐wave (*R* waves) dispersion to infer Martian lithospheric seismic velocity structure and internal temperature profile. To achieve this goal, it is very critical to obtain high‐resolution dispersion measurements for future Martian seismic data. When the seismometer distribution is dense in space, multichannel surface‐wave analysis using slant stacking can be used to obtain dispersion curve (Mcmechan and Ottolini, 1980; Stoffa *et al.*, 1980; Mcmechan and Yedlin, 1981; Nazarian and Stokoe, 1984; Xia *et al.*, 2007). The multioffset phase analysis can also be used (Strobbia and Foti, 2006). All these methods produce nonuniform resolution in the dispersion measurement and the resolution is the poorest at low frequencies. In fact, this dispersion measurement issue was previously recognized and several attempts had been made to sharpen the resolution using a high‐resolution Radon transform (Luo *et al.*, 2008), mode matching (Ross and Lee, 2012), or the power operation of dispersion map (Mcfadden *et al.*, 1986; Shen *et al.*, 2015). However, these attempts to increase the measurement resolution are not significantly beyond the scope of the linear signal analysis and can only achieve limited success. Furthermore, the resultant resolution is still relatively poor at low frequencies and nonuniform.

Rather than using the LSC, we propose a new idea using nonlinear signal analysis to attack this resolution problem. Here, we define the measurement resolution as the ability to resolve (i.e., resolving power) the phase shift between two identical but slightly time‐shifted signals. The resolving power depends on both the data quality and the method used. Obviously, whether the method is sensitive enough to tell a small time shift between the two signals is critical for the high‐resolution measurement of surface‐wave dispersion. Unfortunately the cross correlation is a very insensitive method if the time shift is much less than the period of the signal, and this limits its usefulness at low frequencies. If the surface‐wave data quality is high and contains no noise, in principle, we should be able to measure the relative time shift for surface waves recorded by two stations to an accuracy on the order of the time sampling interval. In practice, depending on the noise level, the resolution should be reduced and it should be greater than the sampling interval.

Before we give detailed exposition of our method, we have to point out that for the surface‐wave dispersion we may expect that the measured low‐frequency dispersion could be more reliable because the low‐frequency waves are less affected by scattering. However, the broad lobes (low resolution) at low frequencies produced by the cross correlation are often interpreted as large errors and low confidence, and this interpretation is clearly problematic.

## Nonlinear Signal Comparison

### Traditional Linear Signal Comparison

We shall consider a surface wave excited by a source that propagates to the receiver 1 and then to the receiver *x*, in the same azimuth. Receiver 1 records a seismogram *d*_{1}(*t*), and receiver *x* records *d*_{x}(*t*). The distance between these two receivers is given (=*x*) and if we can measure the relative time delay between the two signals (at some frequency), we can obtain the phase velocity *V*_{ph}, defined as the distance divided by the measured time delay. Often both the time delay and the resultant phase velocity are frequency *ω* dependent. In the traditional method, we search for a range of *V*_{ph} and calculate the resultant time delay. Then, we shift the trace *d*_{x}(*t*) accordingly and compute their cross correlation: (1)in which *S*_{LSC}(*ω*,*V*_{ph}) is the correlation for the frequency *ω* and the phase velocity *V*_{ph}(*ω*); *d*_{1}(*t*;*ω*) and *d*_{x}(*t*;*ω*) are the seismic waveforms at the frequency *ω* in the time *t* domain (can be thought of as filtered data); is the time‐shifted trace according to the distance *x* between receiver *x* and the reference trace at receiver 1 (*x*=0) and the scanning phase velocity *V*_{ph}; *T* is the length of the time window of interest; *σ*_{1} and *σ*_{2} are the variance of the signals defined as (2)

We expect that at the true time shift, the two signals should achieve the maximum cross correlation for *S*_{LSC}(*ω*,*V*_{ph}). To see more clearly how *S*_{LSC} changes with the frequency *ω* and the time delay *τ*, we consider two cosine signals and calculate *S*: (3)In general, at the low frequencies, *S*_{LSC}(*ω*,*τ*) is insensitive to *τ* (Fig. 1a). It means that to use *S*_{LSC}(*ω*,*τ*) as a measure to resolve the small time shift at low frequencies is not particularly useful. On the other hand, the resolving power of *S*_{LSC}(*ω*,*τ*) is nonuniform across different frequencies.

### Nonlinear Signal Comparison

To overcome the resolution limitation and to achieve a uniform resolution across a wide band of frequencies, we propose a new NLSC similarity measurement: (4)in which (5)in which the variance *σ*_{1} and *σ*_{x} are defined in equation (2). In the new measure (equation 4), *σ* is a nonnegative real‐valued parameter to control the overall resolution. We define a background value *S*_{π} (see the Appendix) for *S*_{NL} when the two signals have a phase difference of *π*: (6)in which *I*_{0} is the modified Bessel function of the zeroth order and (7)Finally, we define the normalized NLSC by scaling *S*_{NL} as (8)The range of *S*_{NLSC} is from 0 to 1. Equation (8) is the formula for our high‐resolution dispersion measurement. We will show that S_{NLSC} has a uniform resolution over a wide band of frequencies and this resolution is controlled by *σ*∈[0,+∞) (Fig. 1b). When *σ*→∞, the NLSC reduces to the traditional cross correlation (i.e., *S*_{NLSC}→*S*_{LSC}), but when *σ*→0, the smeared dispersion lobe shrinks to a geometric line with zero thickness, which is the dispersion curve.

## Synthetic and Field Data Examples

### Example 1: Rayleigh‐Wave Overtones for a Multilayered Model

In the first example, we consider a layered model with elastic layers (Fig. 2a). Rayleigh‐wave overtones are modeled to test the NLSC for multimode surface waves. We model only the fundamental, the first and second overtone using the method by Herrmann (2013) (Fig. 2b). The traditional dispersion analysis *S*_{LSC} (Fig. 2c) is compared with our *S*_{NLSC}(*ω*,*V*_{ph};*σ*) for different *σ* values (Fig. 2d–f). The three Rayleigh‐wave modes can be clearly seen in both traditional LSC as well as in our NLSC. With decreasing *σ*, *S*_{NLSC} resolution has been progressively sharpened and the dispersion curve can now be picked more readily. However, if *σ* is too small, frequency discretization features can be seen (Fig. 2f). NLSC can allow us to choose a continuous range of *σ* to select the dispersion curves.

### Example 2: Global Rayleigh‐Wave Dispersion for Mars Using a Single Station

The National Aeronautics and Space Administration (NASA) InSight mission will put a seismometer on Mars in the near future (Panning *et al.*, 2015). Zheng *et al.* (2015) showed that with just a single seismometer on Mars, we may be able to infer if there is a possible low‐velocity zone (LVZ) in the Martian lithosphere, if we can measure the Rayleigh‐wave group velocity dispersion accurately. In particular, we can measure the dispersion using global Rayleigh waves *R*_{1} and *R*_{3} (or *R*_{2} and *R*_{4}). The LVZ is related to a possible large thermal gradient in the lithosphere. Therefore, if we are able to detect such an LVZ using seismology, it has an important implication in the inference of Martian internal thermal structure and its planetary evolution. However, the key step is to achieve a high‐resolution dispersion measurement. Here, our goal is to demonstrate the ability of NLSC in extracting a high‐resolution dispersion curve using just the seismic recording from one station. First, we generate synthetic seismograms using the 1D Martian seismological model constructed by Zheng *et al.* (2015) with an LVZ. In their paper, Zheng *et al.* (2015) used Mineos (Masters *et al.*, 2011) and the direct solution method (DSM; Geller and Ohminato, 1994) to calculate synthetic seismograms. Here, we use the 3D spectral element method (Komatitsch and Tromp, 2002) to calculate the synthetic seismogram (Fig. 3). We also added noise to the synthetic seismic data (Fig. 3b,c). We use *R*_{1} and *R*_{3} to measure the dispersion. Because the wave propagation distance for *R*_{1} and *R*_{3} is the great circle distance of Mars, which is independent of the epicenter, these two phases can be used to measure the dispersion curve and detect the possible global LVZ as a function of frequency. We obtain the dispersion map using the data with noise (Fig. 3b,c) by the *S*_{LSC} (Fig. 3d) and by the *S*_{NLSC} (Fig. 3e,f). To test the stability of NLSC performance under the influence of noise, we produce 100 dispersion maps with 100 different random noise realizations and stack them together (Yin and Yao, 2016). NLSC shows a significant improvement of the measurement resolution (Fig. 3f), which can facilitate the phase‐velocity picking. Based on the picked phase‐velocity dispersion (Fig. 3f), we can compute the group velocity dispersion (Aki and Richards, 2002; p. 255; Fig. 3g). The maximum of the group velocity around 100 s indicates the existence of the LVZ. Our result here is consistent with the dispersion curve extracted by Zheng *et al.* (2015) using normal‐mode Mineos and the DSM codes.

### Example 3: Ambient Noise Data

The ability to extract surface waves by cross correlation of the ambient noise between receivers is a significant development in the last decade in seismology. The obtained surface waves can provide valuable new information about the subsurface velocity structure without the natural earthquake sources (Campillo and Paul, 2003; Shapiro *et al.*, 2005). In this example, we cross correlate the noise (12 months of noise data) recorded by three stations of the USArray to extract the station pairwise Rayleigh waves (Fig. 4). The noise cross correlation yields two seismograms. Because the three stations are along a line, we use the two‐station method to extract the phase‐velocity dispersion. We calculate the dispersion map by computing *S*_{LSC} and *S*_{NLSC} with different *σ*s (Fig. 4) and compare the results. LSC produces a dispersion map with a broad lobe, especially at the low‐frequency end (Fig. 4a). This has been commonly observed in surface‐wave dispersion measure in numerous previous studies. Often the broad lobe is interpreted as error in the phase‐velocity measurement. This is incorrect, as noted earlier in the Introduction. With our NLSC method (Fig. 4c–f), with decreasing *σ* value, we are able to achieve a uniform high‐resolution dispersion map, which can greatly improve the dispersion curve identification.

## Discussions and Conclusions

In addition to the surface‐wave dispersion measurement, NLSC has a wide range of potential applications in which a signal comparison is needed. For example, in seismic imaging and migration, we need to cross correlate (i.e., compare) the downgoing wavefield from the source and the backpropagated wavefield from the receivers at the image point. The insensitivity of the cross correlation directly translates to low spatial resolution at low frequencies. It is possible to use NLSC to achieve uniform and high spatial resolution in seismic imaging.

In our article, we point out that the poor resolution in surface‐wave dispersion measurement at low frequencies is due to specific LSC techniques such as cross correlation and linear stacking. We propose a new method NLSC to produce high‐resolution dispersion measurement. Furthermore, our new method can achieve a uniform resolution across a wide band of frequencies, which is controlled by an adjustable frequency‐independent parameter, ranging from zero to infinity. When the parameter approaches infinity, the NLSC reduces to the traditional cross correlation. On the other hand, when the parameter approaches zero, an infinite sensitivity can be achieved. We demonstrated the effectiveness and performance of the NLSC using a number of synthetic and field data examples, in the context of global seismology and planetary seismology for future Martian seismological mission.

## Data and Resources

We used the SPECFEM3D_GLOBAL code to calculate wavefield in Mars. The code is publically accessible from the Computational Infrastructure for Geodynamics (CIG) repository https://geodynamics.org/cig/software/specfem3d_globe/ (last accessed February 2016). We used the Mineos package to calculate the wavefield in Mars based on normal mode theory. Mineos is accessible also from the CIG repository (developed by Masters, Barmine, and Kientz, v. 1.0.2; https://geodynamics.org/cig/software/mineos/#release, last accessed February 2014). The ambient noise data used in Figure 4 were downloaded from Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC).

## Appendix

*S*_{π} is calculated using equation (4) by shifting one of the signals by *T*/2, in which *T* is the period of the signal corresponding to the angular frequency *ω* (*ω*=2*π*/*T*): (A1)in which *I*_{0}(·) is the modified Bessel function of the zeroth order of the third kind. The newly defined (A2)is from 0 to 1 for a single frequency wave.

## Acknowledgments

We thank Associate Editor Cleat Zeiler for his insightful reviews and an anonymous review for useful comments. We thank Yao Yao for providing us the USArray ambient noise correlation used in Example 3. We thank Leon Thomsen for useful discussions. We applaud the effort by the Computational Infrastructure for Geodynamics and the authors of SPECFEM3D and Mineos for making their code available for public use. All the data used in the examples will be made available to anyone upon written request to the authors. This work is partially supported by National Science Foundation EAR‐1621878.

- Manuscript received 27 July 2016.