Characterization of Seismic Responses in Mexico City Using Hilbert-Huang Transform

Abstract

In this investigation, we present the Hilbert-Huang transform (HHT) as an alternative technique, which has advantage over other methods for extracting useful data of seismic ground response. The HHT, integrated with the empirical mode decomposition (EMD) and the Hilbert transformation (HT), enables engineers to analyze data from nonlinear and nonstationary processes. The product of the transformation is a detailed description of time-varying frequency diagrams. The recordings of accelerations of soft-soil deposits in Mexico City are studied under this technique. Results of the analysis of accelerograms indicate that this adaptive decomposition permits the extraction motion characteristics, which cannot be effectively unraveled by other conventional data processing techniques. The findings and conclusions derived from studies such as the one presented here contribute to a better understanding of seismic response patterns.

1. Introduction

The assessment of ground motions due to potential earthquakes is perhaps the most challenging problem in geotechnical earthquake engineering. Analysis of available behavioral data collected during seismic phenomena is the natural first step to approach this demand. During any seismic event, there are variations of the intensity of the ground motion with time (acceleration, velocity or displacement). After the arrival of the first seismic wave, the intensity accumulates quickly till a maximum value is sustained (for some seconds) and later diminishes gradually until it gets vanished. Additionally exists a variation of the frequency content with time with a tendency to move to lower frequencies as time progresses, what is known as “frequency dependent dispersive effect” [1].

Recognized that linear transformed domain and response spectral analysis be ill with shortcomings when applied to ground motions [2, 3, 4, 5], this study proposes the use of the Hilbert-Huang transform (HHT) [6, 7] as an alternative for finding useful features into seismic recordings, granted its well-supported validity for dealing with nonstationary random processes [8, 9, 10, 11]. The subject under analysis is a well-known soft-clay deposit in Mexico City, the SCT site [12, 13, 14]. The seismic response of this location is examined in order to achieve detailed information about behaviors that could not be discovered with restrictive tools. The well-defined descriptions of the time-frequency-amplitude content obtained using the HHT better describe the nature of the ground motions. The HHT applied to seismic phenomena results in helpful data for safer designs and improved practical earthquake engineering.

2. Hilbert-Huang transform

2.1. Overview of Hilbert-Huang transform

The HHT was proposed by Huang and coworkers [6] and consists of two parts: (1) empirical mode decomposition (EMD) and (2) Hilbert spectral analysis. Signals to be analyzed are decomposed into a finite number of intrinsic mode functions or IMFs in what is called the EMD process, the crucial part of the transformation. An IMF is described as a function satisfying the following conditions: (1) the number of extrema and zero-crossings must either equal or differ at most by one; (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero [6]. An IMF permits directly Hilbert transforms.

2.2. Hilbert spectral analysis

The Hilbert transform allows the computation of instantaneous frequencies and amplitudes and describes the signal more locally. Eq. (1) displays the Hilbert transform ŷ which can be written for any function x(t) of Lp class (the Lp spaces are function spaces described by a natural generalization of the p-norm for finite-dimensional vector spaces) [5]. The PV denotes Cauchy’s principle value integral.

$$H\left[x\left(t\right)\right]\equiv \widehat{y}\left(t\right)=\frac{1}{\pi }\mathit{PV}{\int }_{-\infty }^{\infty }\frac{x\left(\tau \right)}{t-\tau }\mathit{d\tau }$$E1

Huang et al. [6] and Attoh-Okine et al. [7] determined that an analytic function can be formed with the Hilbert transform pair as shown in Eq. (2).

$$z\left(t\right)=x\left(t\right)+i\widehat{y}\left(t\right)=A\left(t\right){\mathrm{e}}^{\mathit{i\theta }\left(t\right)}$$E2

where

$$A\left(t\right)={\left(x^2+{\widehat{y}}^2\right)}^{1/2},\theta \left(t\right)={\text{tan}}^{-1}\left(\frac{\widehat{y}}{x}\right)\text{and}i=\sqrt{-1}$$E3

A(t) and θ(t) are the instantaneous amplitudes and phase functions, respectively [15]. The instantaneous frequency can then be written as the time derivative of the phase, as shown in Eq. (4).

$$\omega =\frac{\mathit{d\theta }\left(t\right)}{\mathit{dt}}$$E4

Note that the analytic function z(t) is the mathematical approximation to the original signal x(t). Due to the fact that amplitude and frequency functions are expressed as functions of time, the Hilbert spectrum, which displays the relative amplitude or energy (square of amplitude) contributions for a certain frequency at a specific time, can be constructed as H(ω,t). Figure 1 shows the complete block diagram of the HHT. Researchers seeking for deeper knowledge in this subject should refer to the works of [7, 16, 17].

Figure 1.

Hilbert-Huang transform operating diagram (modified from [18]).

2.3. Empirical mode decomposition

The EMD is the first stage in the HHT method. The algorithm starts to decompose almost any signal into a finite set of functions. The resulting functions are called intrinsic mode functions (IMFs) whose Hilbert transformation gives physical instantaneous frequency values. The algorithm employs an iterative sifting process which successively subtracts the local mean from a signal, as follows:

1. Settle on the local extrema (maxima, minima) of the signal.

2. Link the maxima with an interpolation function, producing an upper envelope of the signal.

3. Link the minima with an interpolation function, producing a lower envelope of the signal.

4. Calculate the local mean as half the difference between the upper and lower envelopes.

5. Locate the local mean from the signal.

6. Iterate on the residual.

These six steps are repeated until the produced signal converges to the definition of an IMF, which is a signal with a zero mean and whose number of extrema and zero-crossings differs by at most one, considered as monocomponent functions with no riding waves [5]. Then, the IMF is subtracted from the original signal, and the sifting process is repeated on the remainder. This is iterated until the final residue is a monotonic function. The ultimate extracted IMF is the lowest frequency component of the signal, better known as the trend. The characterization of an IMF was created to guarantee that the signal gives physical frequency values when using the Hilbert transform.

Once a signal has been fully decomposed, the signal D(t) can be written as the finite sum of the IMFs and a final residue, as shown in Eq. (5).

$$D\left(t\right)=R_n\left(t\right)+{\displaystyle \sum _{j=1}^n}{\mathit{IMF}}_j\left(t\right)$$E5

Using Eqs. (5) and (7), the analytic function can be formed as shown in Eq. (2).

$$D\left(t\right)-R_n\left(t\right)=\text{Re}\left[{\displaystyle \sum _{j=1}^n}{A}_j\left(t\right)e^{i\int {\omega }_j\left(t\right)\mathit{dt}}\right]$$E6

Also, for reference, Eq. (6) shows the Fourier decomposition of a signal, x(t).

$$D\left(t\right)=\text{Re}\left[{\displaystyle \sum _{j=1}^n}{A}_j{e}^{i{\omega }_{j}t}\right]$$E7

The EMD decomposition can be considered a generalized Fourier decomposition, because it describes a signal in terms of amplitude and basis functions whose amplitudes and frequencies may fluctuate with time [5, 16]. Detailed examples of EMD can be found in [16].

3. Hilbert-Huang transformation of seismic response

3.1. Geotechnical database

To illustrate the HHT for studying seismic response of clayey soils, data recorded in the SCT site, a location in Mexico City, are analyzed. SCT is located within a densely urbanized zone, where earthquake-related damages have concentrated recurrently in the past. The data set selected for this investigation consisted of a set of 22 acceleration time series recorded during earthquakes that happened from 1985 to 2011. The characteristics of these events are enlisted in Table 1. The signals are of 150 s length (on average) at a measuring rate of 0.01 s. Magnitudes go from M3 to M8.1, and epicenters are from the Pacific subduction zone (in front of the Guerrero and Michoacán coasts) and normal intraplate earthquakes (located in central Mexico).

Table 1.

Characteristics of seismic events.

The formation is a sequence of soft clay strata interspersed with layers of harder clayey silts with sands (Figure 2a). The evolution of pore pressures, measured over a 12-year period (1990–2002) [19], illustrates a gradual depletion of pressures from 0.014 to 0.002 kPa/year at different depths. The water content, volumetric weight and undrained strength profiles at two different dates (1952 and 1986) are shown in Figure 2b. These profiles show that (1) thicknesses of the relevant clay strata have decreased as a consequence of regional subsidence; (2) water content reductions are especially significant below 15 m approximately which is consistent with the fact that pumping-induced consolidation propagates upwards from the base of the clayey soils to the surface; and (3) changes in soil density (volumetric weight) confirm that soil strata have densified as the clay masses consolidate.

Figure 2.

Stratigraphic column in SCT site, the evolution of water content W, volumetric weight ɣ, undrained strength q_c and shear wave velocity Vs profiles. (a) Description (b) W, % (c) ɣ, kN/m³ (d) q_c, kPa (e) Vs, (m/s).

Concerning the dynamic properties, the evolution on shear wave propagation velocities is presented through the field evidence given in Figure 2d (suspension logging tests performed in 1986 and 2001 [20]). The expected stiffening of the clay strata over the 15 years that lapsed between dates is more evident for greater depths.

3.2. Decomposition of accelerograms

According to the procedure described above (EMD algorithm), the selected accelerograms were decomposed into their IMF components. The results are like those shown, by way of example, in Figure 3 (IMFs obtained from the September 30, 1999 event). For this accelerogram, the algorithm yielded nine components and a residue.

Figure 3.

IMFs obtained from September 30, 1999, event.

To discriminate between the IMFs, they will be studied under two geotechnical-seismological parameters that are crucial for anti-seismic designs: fundamental frequency f_n and peak ground acceleration (PGA).

In Mexico City, several past earthquakes have resulted in devastating failures of built environment, which in most of the cases were attributed to strong amplification of seismic waves. For analyzing the phenomena of resonance, it is necessary to determine f_n of a soil deposit. f_n is dependent on its thickness, low strain stiffness and density. Assessment of these characteristics has been the subject of several research studies, yet inconsistencies and uncertainties associated with their estimation have not been resolved conclusively [21, 22]. Hence, many researchers tend to corroborate their calculations by using data recorded on surface during earthquake events.

The best tactics to establish f_n are to compute it using an abstract method and then compare this value with that shown by recorded data. Theoretically, fundamental frequency for SCT (clay deposit ≈ 40 m thick) calculated using average of shear wave velocity is f_n = 0.5 Hz. Using microtremors [23, 24] and multi-degree-of-freedom system (MDOF) methods [25], the computed values occur in the neighborhood of f_n ~0.47 Hz.

Attempting to relate the intrinsic oscillations with this important deposit parameter, the predominant frequency of each IMF in data set is analyzed. What follows this predominant frequency will be called f_IMF. In Figure 4, boxplots of the f_IMF corresponding to each mode are depicted. The first (IMF1) and last modes (IMF6 to IMF10) are not in the graph. The IMF1s have flat spectra where it is not possible to determine a predominant frequency. The last modes (from IMF6 to IMF10) have very short box sections with exceptionally low frequency values (near to zero). These modes, in terms of their fundamental frequency, are considered irrelevant, that is, they do not have effect on the analysis of the measured response.

Figure 4.

Boxplot of IMF, excluding the corrupting and irrelevant modes.

Despite the range of magnitudes, epicentral depths, generating mechanisms, directivity and dates, the f_IMF of IMF5 is closely related to the natural frequency of SCT site. The f_IMF of IMF5 practically does not change, and this is noticeable in the condensed box. The dispersion is minimal, and it should be noted that this graph includes minor events (M < 4) and the extreme earthquake of 1985 (M8.1). On the other hand, the median of the f_IMF of IMF3 is comparable to the second vibrate mode f₂ (using the MDOF method [23], f₂ is ~1.4 Hz). Its boxplot is taller than the IMF5 and has evident dispersion in lower and upper whiskers.

Mexico City’s seismic response has changed over the years. The numerous efforts made to analyze this phenomenon showed that the response may be due to land subsidence (consequence of water extraction) [19, 20]. Since the pumping and the capacity of the clay strata to expel water are quite different at each geotechnical zone of the city, the effect on seismic response is a complex problem.

To guarantee safe structural designs, the most significant aspect to be demonstrated is the effect of the subsidence on the spectral ordinates (with respect to the current values specified for design purposes) and the relative values of the structures and site periods. With specific geotechnical and geological information, Avilés et al. [26] presented an incremental procedure to evaluate the progressive evolution of f_n over time, in order to evaluate the effects of these changes on the earthquake ground response. The adjustments expected in selected values of f_n, from 1980 to 2020, can be appreciated in Figure 5. What can be anticipated, as the evolution of this phenomenon, is that after the initial reference year the parameter decreases monotonically as the times span increases [27].

Figure 5.

Evolution of f_n, continuous line results from [26]. Dashed line represents the behavioral trend of f_IMF-IMF5.

Observe what occurred when the values of f_n, from original data, are superimposed in the graph. The frequencies obtained from the recorded accelerograms hardly follow the trend set by [26]. Most events do not fit the specific curve for SCT (1/f_n = 2), and the expected drop patterns of f_n cannot be verified. However, the behavior shown by the IMF3 and IMF5 fits very well to the decrement curves obtained with the iterative model. The f_IMF of IMF5 fits the theoretical curve of 1/f_n = 2.0. This finding reinforces the hypothesis about stiffening of the clayey strata that are not possible to verify using the f_n from registered data. Actually, a deep analysis of the path shown by f_IMF-IMF5 could drive to two stages of the evolution: the first one with a more dramatic decrement of the 1/f_n since 1985 (or before) and toward the year 2000, and a second one where the frequencies change more slowly to the present day. The model’s tendency that drops the values of 1/f_n less dramatically over time is clearly illustrated with the aligned values of f_IMF-IMF3.

PGA, as the maximum ground acceleration that occurred during earthquake shaking at a location [28], is used as an intensity measurement, and the design-basis-earthquake ground motion is often defined in terms of PGA [29]. Even damage to buildings and infrastructure is more closely related to ground motion, of which PGA is a measure, rather than the magnitude of the earthquake itself [30].

The boxplots of the maximum instant amplitudes, associated with each IMF, are shown in Figure 6a. The most extreme points (marked with the numbers 2 and 15) are two of the most intense seismic events registered in Mexico City (M > 7.0), being the extreme value the devastating 1985 Michoacán earthquake (M8.1).

Figure 6.

Boxplot of MIAs according to the IMF, (a) the whole set, including 19/09/85 Michoacán earthquake and (b) events with M ≤ 7.5.

Isolating the M8.1 event to properly observe the characteristics of each box (Figure 6b), it is evident that last modes have very small boxes. Considering this minor amplitudes and their fundamental frequencies (very close to zero), we conclude that they are irrelevant modes. Examining the relatively small box of the IMF1 maximum amplitudes, it was found that the extraordinary earthquake of 1985 is found inside the rectangle. It seems that this first oscillation does not respond to the input energy as the other modes do, so it is inferred that a mechanism, other than the ground motion itself, is the one that controls this oscillation. Based on this observation and their flat Fourier spectra, the IMF1 is labeled as a mode that “corrupts” the registered motion of the soil deposit.

It can be seen that the energy is concentrated mainly on fourth and fifth oscillations. These modes most react to the extreme earthquake of 1985. The dispersion is minimal, and the atypical event (M8.1) is certainly at a distance, which does not happen with the other oscillations. It appears that the recognized amplification potential of this soft-clay site [20] is mainly evident for the central modes.

The “corrupted” (IMF1) and the irrelevant modes (IMF6 to IMF10) were removed from the registered accelerograms. The response spectra of accelerations from both sets of signals, the original and the “cleaned,” are calculated. Some selected responses are depicted in Figure 7. These are examples of responses that can be reproduced with the IMF3 or IMF5 with adequate approximation. The energy-frequency content of the original response spectra is represented with one intrinsic mode, the one with a MIA near to the PGA registered. It has been found that the most intense events from the Pacific subduction zone tend to be related mainly with the fifth mode, while for minor and intraplate events, the ground motion can be represented using only the third mode. Deep examinations are necessary to clarify this finding, but the advantages of EMD when analyzing complex responses are clear.

Figure 7.

Examples of spectra evaluated using original signals and selected IMFs, (a) events better represented with IMF5 and (b) events better represented with IMF3.

It should be emphasized that the “cleaned” signals represent simpler and easier time series, with desirable density for minimizing computational costs (when running complex ground-motion models).

3.3. Hilbert spectra of ground motions

The nonrecommended practice of using Fourier for studying nonstationary-nonlinear data is revealed as a limitative tool for these kinds of analyses since it cannot show the detailed information about the peculiarities of the energy-frequency distribution, aspect that is utilized, for example, to conclude about site effects.

From the Hilbert spectra, events in the database have direct important characteristics such as the arrival time and the duration of intense phase of the earthquakes can be detected. In addition to the well-defined zone of maximum accelerations, the huge variations of the response frequencies between the beginning (the maximum amplitude cycles) and the end of the event can be clearly detected and quantified

(see, for example, the HS events in Figure 8). Although most of the events found in the database have disseminated energy in the high-frequency range, this is minor compared to the highest levels recorded in the entire time series. The persistent energy resides along horizontal packets below 1 Hz. It is easy to detect the peak accelerations of each response (in addition to the concentric energy surrounding these maxima), the fundamental frequencies, and the duration of the arrival packets. The “corrupted” oscillations can be distinguished, with low frequencies and low constant amplitudes that are maintained from the beginning to the end of the record.

Figure 8.

The HS for selected intraplate and interplate events. Significant frequency deviation during the strong motion showed that the energy is concentrated below 1 Hz.

The difficulty of transforming signals from near-field motions (time series of short duration) to Fourier spectra is not a problem using HS (see, for example, the HS of the June 15, 1999 event). The frequency distribution clearly shows the energy content along the frequency axis (from 0.4 to 0.8 Hz), and it is quite easy to determine the resonant frequency of the site (at ≈0.55 Hz) and the associated energy density.

In Figure 9, three of the most intense earthquakes registered in Mexico City are shown. They are the well-known 1985 and the most recent devastating events, September 8, 2017, and September 19, 2017. While the responses to minor earthquakes have a much wider energy distribution, the content for these events is quite different. The spikes are more evident, and they can be located around the PGA. Although these Hilbert spectra have the decreasing trend in energy density as a function of frequency, as the other events, these spectra show a visible lack of energy in the low-frequency range (less than 0.2 Hz).

Figure 9.

The historical 1985 earthquake (M8.1) and the responses registered during September 8, 2017 (M8.2) and devastating September 19, 2017 (M7.1).

Consider the HS of September 19, 2017 earthquake, the one that hit the city and collapsed a large number of buildings, houses and underground facilities. It is well defined at the low-frequency range (from 0.4 to below 1 Hz) where maximum amplitudes were present, a serious situation for the high-rise structures. The peak of energy, located at a very narrow frequency range (0.5–0.8 Hz), could be responsible for the huge number of damaged buildings in a defined geotechnical area of the city. Similarly, to the spatial variation of the destruction at Mexico City during the 1985 earthquake, for this event some areas have been recognized, with special stratigraphies, which presented amplification potentials never seen before. As many investigators have concluded, any linear model could not predict this kind of ground responses.

4. Conclusions

The decomposed components, namely IMFs, contain observable, physical information about the ground motion. Owing to the fact that the intrinsic oscillations can be related to the phenomena without restrictive hypothesis, the conclusions obtained from the behavior of each mode are more natural than those achieved from theoretical positions faraway the system, which generates the signal.

Despite the fact that the soil properties in the SCT site have changed because of regional subsidence, the recorded accelerations do not seem to exhibit the expected stiffening effect. What is shown in this research is that the time series, which serve as a basis for showing the soils hardening, are corrupted by external oscillations that mask actual responses. It was recognized that IMF5 and IMF3 actually showed the decrease in rate of the frequencies postulated in theoretical models. By detecting the oscillations that do mainly correspond to the real ground motion, the EMD is presented as an alternative tool for finding directly the frequencies related with the highest energy concentrations.

The usefulness of these finding points to the need for further investigation in evaluating the competency of tools for analyzing the recordings and for obtaining conclusions, since they are the basis of the process of calibration of all the theoretical, geotechnical and seismological models.

The Hilbert spectral transformation permits the recognition of oscillation modes that would have been masked by Fourier spectral analysis.

Acknowledgments

The accelerogram records for this research were provided by the Unidad de Instrumentación Sísmica, Instituto de Ingeniería, UNAM.