Surface-Wave Dispersion Measurements Using Hilbert-Huang Transform

An efficient procedure using the Hilbert-Huang transform (HHT) method for the determination of the dispersion curves of seismic surface waves is investigated in this study. The HHT method is compared with a Fourier-based time-frequency analysis in the extraction of dispersion data from a synthetic seismic waveform generated by a point source in a layered medium. The HHT method provides high-resolution measurement of the spectral content as a function of arrival time, and the dispersion curve of group velocity is determined more accurately than that derived using con­ ventional Fourier-based frequency-time analysis. Using a simple polyno­ mial fitting technique, the dispersion curve of the phase velocity can also be derived by integrating the dispersion curve of the group velocity. Numeri­ cal tests show that the HHT method is capable of mapping the energy dis­ tribution in the frequency-time domain accurately, and the procedure re­ quires much less computer time than other methods for surface-wave dis­ persion analysis.


INTRODUCTION
Surface-wave dispersion measurement is essential to studies of the crust and upper mantle structures, earthquake source mechanisms, and inelastic properties of the earth.Several methods, essentially all Fourier-based, have been in common use to estimate the dispersion curves since the first analysis on this dispersion phenomenon was studied by a simple peak-and-trough method (Ewing and Press 1952).It measures the arrival time of peaks, zeros, and troughs for each oscillation corresponding to a narr ow-band wave packet with an average period given by the period of the particular cycle.Nonetheless, for most observations, the dispersive wave-form is composed of closely packed energy arrivals of adjacent periods, making it virtually impossible to measure the group velocity at a certain period.Sato (1955) was among the first to measure this dispersion property applying the Fourier transformation to the signal.His method was limited to highly-dispersed and high signal-to-noise ratio data.The moving win dow analysis (MWA, i.e., Landisman et al. 1969) and the multiple filter technique (MFT, i.e., Dziewonski et al. 1969) are two different but both Fourier-based methods that are often used to measure surface-wave dispersion.The filtering techniques applied to these methods are equivalent: MWA in the time domain, and MFT in the frequency domain, both methods ex tend the resolution of dispersion measurement to a broader period range.They are particularly useful in treating dispersive signals over continental paths where the dispersion curves are nearly flat over a broad period band, in which wave energies arrive at nearly the same time.
However, there is a difficult resolution problem inherent in all Fourier-based analysis: the product of frequency uncertainty and arrival time uncertainty is a constant (Feng and Teng 1983).Several filtering techniques to optimize the resolution of dispersion measurement have been proposed (e.g., Inston et al. 1971;Cara 1973;Nyman and Landisman 1977).Unfortunately, the goal of improving the resolution of dispersion measurement by these optimization tech niques is still limited by the above-mentioned uncertainty principle of the Fourier analysis.
There is, .inaddition, some basic difficulty in the signal stationarity and linearity assumptions on which the Fourier theory is based.For example, the span of any seismic data is time limited: clearly, seismic data are non-stationary.Moreover, the seismic waves coming through the attenuating earth and recording instruments, both of which may not be distortion-free.The linearity assumption inay be compromised.These problems are usually ignored for most Fou rier-based analysis.One has to assume the data to be piecewise stationary and nearly linear, which is not quite appropriate to the seismic data when the MWA or the MFT is applied.
Furthermore, the conflicting requirements occur unavoidably, when the frequency resolution requires longer time series but the time resolution requires a narrow window width.
In order to overcome the limitations presented by the traditional Fourier-based analysis, the wavelet analysis (WA), which itself is also Fourier-based and essentially an adjustable window spectral analysis, was introduced.Among others, and in seismology, Chakraborty 1996 used the wavelet transformation and related methods to perform the frequency-time analy sis of seismic reflection data.Pyrak-Nolte and Nolte (1995) applied the wavelet transforma tion to measure the group velocity of the long-period Rayleigh waves.Although these studies showed that the wavelet transformation gives better resolution in the analysis of non-station ary data, the problems of non-adaptive nature, which fits the data only when we have pre knowledge of the wavelet physically, and physically meaningless to the linear phenomena, are still not solved completely.Huang et al. (1996) introduced a new approach called the Hilbert Huang transform (HHT) technique to perform general analysis of highly transient time series from nonlinear, dispersive, and non-stationary systems.For extremely transient and some what nonlinear seismic data, this method provides a proper tool to investigate the characteris tics of ground motion signals in frequency and time domain simultaneously.In this paper, we demonstrate an efficient determination of the group velocity dispersion curves with the HHT technique.To show the improvement in the resolution and accuracy of the dispersion measurement, the results are compared with those obtained by the conventional MFT method.

Hilbert-Huang Transform
Unlike the traditional Fourier analysis, in which the frequency is defined for the sine or cosine function spanning the whole data length with constant amplitude, a new method called the Hilbert-Huang transform (Huang et al. 1996) was developed for analyzing nonlinear and non-stationary time-series.The key to the method is a process Huang et al. called Empirical Mode Decomposition (EMD) by which any complicated time series can be decomposed into a finite and often small number of Intrinsic Mode Functions (IMF) that admit well-behaved Hilbert transforms.By applying the Hilbert-Huang transformation to get the IMFs ( Xj (t) ), an arbitrary time series X(t) can be expressed as where n is the instantaneous frequency, dr , re -=t-r P is the Cauchy principal value of the integral, Since both the amplitude and frequency are functions of time, the Hilbert spectrum, H(m, t), which is the amplitude contoured on the frequency-time plane in a three-dimensional plot, is more appropriate to accommodate non-stationary data than the Fourier spectrum, because the Fourier expansion is made in a global sense and HHT expansion is made in a local sense.
Since the Hilbert transform gives the best-fit to a local sine or cosine waveform to the data, the frequency resolution for any point is uniformly defined by the stationary phase method or local derivative of the phase function.This advantage is especially effective in extracting the time-varying information of an arbitrary signal (Huang et al. 1996).

Procedure
To determine the dispersion curve of the group velocity of seismic surface waves, we directly apply the EMO decomposition to the seismogram to get the intrinsic mode function (IMF), from which the Hilbert spectrum is calculated.The Hilbert spectrum usually shows sharp-banded amplitude contours illustrating energy distribution of IMFs in the frequency time space; this represents the dispersion of group velocity for the waveform concerned.For each frequency, the time of the maximum amplitude corresponds to the group arrival at which the group velocity can be obtained; moreover, the velocity error can be estimated from the width of the amplitude contour band in the spectrum.According to Papoulis (1962), and where �r is the travel time of a wave group, and cp is the phase function.Then, the phase velocity can be obtained as where X is epicentral distance, ell S R is the initial phase of the source.In order to accomplish the integral, we extrapolated the dispersion curve of the group velocity to lower frequency (even zero Hz) using a polynomial fitting, and obtained the dispersion curve of the phase velocity through integration.Owing to the underestimation of phase velocity in the low frequency range where few data can be obtained from the spectrum, the phase derived from the integral of group delay time is over-estimated.Similar to the correction used in the single station method (Brune et al. 1960), we corrected the derived phase by subtracting a integer (N ) to constrain the phase velocity of low frequency to a reasonable value.In general, N is small ( <5) and approaches zero as the sharp velocity variation of dispersion curve occurs at lower frequency (e.g., fundamental mode shown in the test below).

TESTS OF THE HHT TECHNIQUE BY KNOWN NUMERICAL FUNCTIONS
We shall use three known time series from Huang et al. (1996) to illustrate the high resolution property of the HHT method.
First, let us consider the time series: which is the summation of three distinct cosine waves with frequencies of 0.1, 0.02, and 0. OOlHz, respectively (Fig. 1).The EMD method successfully decomposed the signal into three IMFs, which are identical with three assumed cosine waves, except for mimute differences on boundaries because of the endpoint effect.The residue is much smaller than IMFs and repre sents a monotonic trend decided by cutoff of the data window.The resulting Hilbert spectrum (Fig. 2), which shows three constant-frequency energy bands at 0.1, 0.02, and 0.001 Hz, is given at the lower middle panel.The vertical bars on the right of this panel give the relative amplitude scale, while the panel on its left gives the marginal spectrum.This example shows that the HHT method can analyze the stationary signal simply from the data without any as sumption of harmonic components (e.g,, sine or cosine waves in the Fourier transform).
In the second example, we consider an amplitude modulated wave given by (2 n t) x(t) = exp(-0.01t)cos W , for t from 0 to 500 sec.The time series is a non-stationary signal and given in the middle of Fig. 3.This amplitude-exponentially-decayed signal has a zero mean and is symmetric, so Data I Signal 0.12.------.---��----------------� 0.1 . . ... ... .... ..,.�=: :"1 there is no need for invoking EMD since the data is already an IMF.The time�frequency energy distributions are given for the Fourier-based analysis (top panel) and for the HHT method (bottom panel).The vertical bars on the right of these two panels give the relative amplitude scale, while the two panels on their left give the Fourier spectrum of the Fourier based frequency-time analysis and the marginal spectrum for the HHT.The Fourier-based analysis gives a smeared average frequency range over which the main wave energy resides due to the harmonic wave presumption.The time variation of the amplitude modulation is also insignificant in the Fourier-based result.On the•other hand, the HHT method gives a much sharper identification of energy distribution in the frequency-time domain.Moreover, the time variation of amplitude decay is clearly shown as fading of energy distribution in the Hilbert spectrum.The amplitude modulation has introduced a small instantaneous frequency modula tion around the mean of the carrier frequency, which can be detected using the HHT method.The range of variation, however, is insignificant compared with the frequency smearing phe nomena of the Fourier-based analysis.
The third example further demonstrates the resolution power of the Hilbert spectrum.Let us consider the case of a simple cosine wave with one frequency (l/16Hz) suddenly switching to another frequency (l/32Hz) at t = 500s; i.e., the time series 113:5 ' 0 " ' 0.35 •.,:.;.: .. :. , .. ... ,,,,, .,., ., .. The vertical bars on the right of these two panels give the relative amplitude scale, while the two panels on their left give the Fourier spectrum and the marginal spectrum.Again, the Fourier spectrum gives smeared average frequency range over which the main wave energy resides.Adding and canceling of the harmonic components using the in Fourier analysis produce artificial energy outside the main frequencies in which we expect all the energy to exist.Also, the Fourier based analysis has to compromise between time and frequency resolution, as a result, the frequency transition at t = 500s is not precise.In contrast, we can see the sharp frequency definitions and the time location of the frequency switch using the Hilbert spectrum.Small frequency fluctuations are observed on the endpoints and point of transition due to the Gibbs phenomena and the finite data length.Fortunately, the range of variation is insignificant.Com parison of the marginal Hilbert spectrum and the Fourier spectrum is shown in the left panels.
Because of the Fourier transform assumption of global frequency distribution (due to stationarity assumption), its spectrum introduces numerous spurious frequency contents that make up many side lobes of the spectrum.The marginal spectrum, on the other hand, gives sharp identifica tion at two main frequencies, as expected.
These three examples clearly illustrate the different resolution powers of Fourier-based analysis and HHT.We shall apply this new method to the extraction of seismic surface-wave dispersion data.

USING HHT TECHNIQUE TO EXTRACT SURFACE-WA VE DISPERSION DATA
To show the capability of the HHT technique in obtaining time-varying characteristics from seismic waves, a synthetic seismogram of the vertical component calculated (Wang and Herrmann 1980) with a layer model (modified from IASPI91 radial heterogeneity earth model, Fig. 5) was used for this test.The source is placed at a depth of 10 km and recorded at the distance of 6000 km.In order to compare the results, the contours in time-frequency domain are generated and compared using both HHT and MW applied to the same vertical component of the synthetic seismogram.The waveform and the spectrum obtained from HHT and MW analysis are shown in Figs. 6 and 7.In Fig. 6, the upper panel shows the synthetic seismic signal and the marginal spectrum (Huang et al. 1996), which corresponds to Fourier amplitude spectrum Gust as shown in Fig. 7), and is plotted in the lower left panel.The 20 (Frequency Time) spectrum is shown in the lower middle panel.For noise free signal, the results from the HHT (Fig. 6) give excellent resolution while the MW results (Fig. 7) give a broad diffused band of energy arrivals.The frequencies at any time instant are pinpointed sharply and accu rately in the HHT results, with the amplitude well represented by the energy level of the en-  A simple transformation gives the group velocity dispersion curve (dots in Fig. 9).Moreover, the time variation of amplitude decay is clearly shown as a fading pattern in energy distribution in the HHT spectrum.The resolution power of the HHT spectrum allows an accurate estimate of group velocity dispersion for the fundamental mode and even for the first-higher mode easily.On the other hand, the result obtained using the MW (Fig. 6) shows a smeared average frequency range over which the main wave energy resides.This smearing is due to the harmonic wave assumption for a stationary time series, and represents the global property, which is non-existence.For an extremely noisy signal (Fig. 8, top trace), 30% ran dom noise is added.The HHT spectrum analysis can still provide the better resolution (Fig. 8 middle lower panel).This demonstrates that the HHT analysis is far superior to the Fourier based moving window (MW) method in the determination of the group velocity dispersion curve of noisy non-stationary surface waves, which is actually what observed in nature.
The dispersion of the group velocity can be easily determined from the HHT energy frequency-time distribution in Fig. 6.While the resolution in frequency is good in a time (or distance) range, we can directly find the arrival time of the maximum of energy at each fre-quency (Fig. 9), which corresponds to the group arrival for each frequency.To determine the dispersion of the phase velocity, we can fit these points with a simple second order polynomial followed by a numerical integration.Applying the procedure described above, the dispersion curve of the phase velocity is also determined (Fig. 9).

CONCLUSION
The HHT technique can yield energy-frequency-time distribution with high-resolution for a non-stationary surface-wave signal, from which the group velocity and phase velocity dispersion curves can be measured to much higher accuracy than obtainable by conventional Fourier-based methods.This improvement has direct bearing on the advanced analysis of crustal and upper mantle structure.Without the sharp identification of frequency, the group velocity measurements using the Fourier-based methods, such as the moving window approach are shown to be much poorer.Although this study shows some encouraging results using the   The HHT method can also give us a satisfactory result for extreme noise data.

183
HHT technique, some improvements such as the polynomial fitting of the group velocity curve and estimation of error of the dispersion curves need to be further investigated.Moreover, a systematic study with a multi-station complete seismic data set is required to examine the full application of the phase velocity dispersion measurement using this HHT method.

Fig. 1 .
Fig. 1.The EMD result showing three IMFs and a residue for an analytical sig nal with summation of three cosine waves with frequencies of 0.1, 0.02, and 0.001 Hz.

Fig. 2 .
Fig. 2. Illustration of the signal (top) and the Hilbert spectrum (bottom).The vertical bars on the right of this panel give the relative amplitude scale, while the panel on its left gives the marginal spectrum.

Fig. 3 .
Fig. 3.The F-T distribution and its corresponding spectrum using the MW method (top) and HHT (bottom) of an amplitude-modulated cosine wave.Note that the HHT method gives a much sharper identification of energy dis tribution in the Frequency-Time domain.

Fig. 4 .
Fig. 4. The F-T distribution and its corresponding spectrum using the MW method (top) and HHT (bottom) of a signal with abruptly-changed frequency.The HHT method again gives a much sharper identification of energy distribution in the Frequency�Time domain.

Fig. 8 .
Fig. 8.The F-T distribution and its corresponding spectrum using the HHT of a synthetic seismogram (used in Figs. 6 and 7) with 30% random noise.

Fig. 9 .
Fig. 9. HHT derived group velocity dispersion data (dots), the theoretical (solid line) and calculated (dash line) dispersion curves of group and phase velocity.