Application of Hilbert-Huang Transform in Generating Spectrum-Compatible Earthquake Time Histories

Abstract

Spectrum-compatible earthquake time histories have been widely used for seismic analysis and design. In this paper, a data processing method, Hilbert-Huang transform, is applied to generate earthquake time histories compatible with the target seismic design spectra based on multiple actual earthquake records. Each actual earthquake record is decomposed into several components of time-dependent amplitude and frequency by Hilbert-Huang transform. The spectrum-compatible earthquake time history is obtained by solving an optimization problem to minimize the relative difference between the response spectrum of the generated time history and the target seismic design spectra. Since the basis for generating spectrum-compatible earthquake time histories is derived from actual earthquake records by employing the Hilbert-Huang transform, the nonstationary characteristics and the natural properties of the seed earthquake records are well preserved in the generated earthquake time histories.

1. Introduction

The Hilbert-Huang transform (HHT) [1] provides a new method of analyzing nonstationary and nonlinear time series data and has been applied to a wide range of disciplines. Current major applications of the HHT include image processing, meteorology and atmospherics, ocean engineering, health monitoring, system identification, and earthquake engineering. In this paper, the HHT will be applied in generating spectrum-compatible earthquake time histories.

Seismic Response History Analysis (SRHA) is a major seismic response analysis method for critical structures, such as nuclear power plants [2–4]. The SRHA procedure is concerned with the determination of structural response as a function of time when the structural system is subjected to a given ground acceleration [5]. Hence, the representative input earthquake time histories are required for this analysis method.

A plot of the peak value of a response quantity, which can be displacement, velocity, or acceleration, as a function of the natural vibration period 𝑇_𝑛 of a linear single degree-of-freedom (SDOF) system, is called the earthquake response spectrum for that quantity. The concept of earthquake response spectra is illustrated in Figure 1 [5]. In the following, the response quantity is taken as acceleration. Based on the earthquake response spectra, the seismic design spectra, such as the conventional design spectra [4, 6] and the Uniform Hazard Spectra (UHS) [7, 8], can be derived by statistical analysis based on a large number of ground motion database [9] or the Probabilistic Seismic Hazard Analysis (PSHA) [10]. The seismic design spectra have been adapted as design ground motions in seismic design codes [4, 6, 11] for decades.

Figure 1 

Concept of earthquake response spectra.

An earthquake time history, whose response spectrum closely matches or envelopes a seismic design spectrum over a range of periods of engineering concern, is usually called a spectrum-compatible earthquake time history. The use of spectrum-compatible time histories instead of actual earthquake records is attractive for multiple reasons: (1) they are able to produce structural responses that present relatively lower dispersion; (2) there is only a small number of actual earthquake records available for many regions in the world [12].

Many research works have been done in generating spectrum-compatible artificial earthquake time histories. One of the most classical attempts is to use a Fourier series representation for the time histories to be generated based on random vibration theory. The sinusoidal motions are then summed, and a subsequent iterative process refines the solution to match the target response spectrum. In these works, various earthquake envelope/shape functions were used to characterize the nonstationary properties of the time histories approximately [13–20].

However, actual ground motions are complicated and contain a wealth of information about the nature of the earthquakes; they are influenced by, and consequently reflect, characteristics of the seismic source, the rupture process, the source-site travel path, and local site conditions. Although it is convenient to characterize them using a small number of parameters for generating artificial ground motions, such characterizations can never be complete [21]. With the increase of available strong ground motion records, using and scaling actual earthquake records has become one of the most popular research topics in this field.

Tsai [22] selected an existing actual earthquake record whose spectrum matches closely with the target design spectrum. The record was then passed successively through certain frequency-suppressing filters to reduce the spectrum ordinates wherever necessary. Similarly, sinusoidal motions were superposed over the selected record to increase the spectrum as required. However, several researchers found it convenient to work in the frequency domain rather than with the time histories in the time domain. The technique is to scale the Fourier amplitudes or phases of the actual record such that the resulting response spectrum is compatible with the target design spectrum [23–25].

However, Fourier-based techniques do not account for the instantaneous variations in the frequency contents of an earthquake time history that arise due to the arrivals of different types of seismic waves at different time instants and the phenomenon of dispersion in these waves. In order to simulate the nonstationary characteristics of earthquake time histories, a wavelet-based procedure has been used. It can decompose an actual earthquake record into a desired number of time history components with nonoverlapping frequency contents; each of the time history components is then suitably scaled for matching of the response spectrum of the modified earthquake record with a target design spectrum [26]. The wavelet-based technique, however, is still not good enough for characterizing the instantaneous variations in the frequency contents of an earthquake time history, since the bases of it are predefined functions and it is unable to locate the frequency distribution of the signal accurately [1].

An alternative method in preserving the nonstationary characteristics of earthquake time histories is to adjust actual earthquake time histories by adding wavelet functions in time domain to match the target design spectra. This method is based on an assumption that the time at which the spectral response of a time history occurs is not perturbed by a small adjustment of the time history [27–31]. The method, however, cannot generate spectrum-compatible time history based on multiple actual earthquake records. Also, the effectiveness of the spectral matching may highly depend on the initial time history.

In this paper, a new generation method [32], which can overcome some deficiencies of previous methods, is expanded for both the uniform hazard spectra and the conventional design spectra. The generated time histories preserve the nonstationary characteristics of the multiple seed actual earthquake records.

In Section 2, the Hilbert-Huang transform (HHT) and its characteristics in processing earthquake time history are first described and analyzed, respectively. In Section 3, HHT is used to generate spectrum-compatible time histories based on multiple actual earthquake records and optimization techniques. The proposed generation method is then demonstrated by numerical examples in Section 4. Some conclusions and observations are discussed in Section 5.

2. Hilbert-Huang Transform (HHT)

To obtain spectrum-compatible earthquake time histories for seismic response history analysis, the Hilbert-Huang transform (HHT), which decomposes an actual earthquake time history into a small number of components, is applied.

The Hilbert-Huang transform (HHT), developed by Huang et al. [1], can represent nonstationary and nonlinear data, such as earthquake records, by decomposing the data into several components and transforming the data from time-domain to frequency-domain.

In the past decade, some researchers have used the HHT in processing and analyzing recorded earthquake ground motions [34–37]. They concluded that the Hilbert-Huang spectral analysis gave the most detailed information in a time-frequency-energy presentation compared to traditional data processing techniques. The HHT is also used in simulating a large number of ground motions by randomly shifting the instantaneous frequencies of the actual earthquake records [38, 39].

Compared to Fourier transform and wavelet transform, the HHT can meet the necessary conditions for the basis to represent a nonstationary and nonlinear time series: complete, local, and adaptive. The first condition guarantees the degree of precision of the expansion. The requirement for locality is the most crucial for nonstationarity, which means all events have to be identified by the time of their occurrence. Consequently, it is required that both the amplitude (or energy) and the frequency be functions of time. The requirement for adaptivity is also crucial for both nonstationary and nonlinear data. It is impossible to expect a predetermined basis to fit all the phenomena in the data. An easy way to generate the necessary adaptive basis is to derive the basis from the data themselves. This is the substantial advantage of the HHT compared to other transform techniques.

The HHT consists of the empirical mode decomposition (EMD) and the Hilbert Spectral Analysis (HSA). The HHT uses the EMD method to decompose a time series into so-called intrinsic mode functions and uses the HSA method to obtain instantaneous frequency data.

2.1. Empirical Mode Decomposition (EMD)

The Empirical Mode Decomposition (EMD) step is based on the assumption that any earthquake record consists of different, simple, and intrinsic modes of oscillation, derived from the earthquake record [1]. Each of these oscillatory modes, called an intrinsic mode function, is defined by the following conditions: (1)over the entire earthquake time history, the number of extrema and the number of zero-crossings must be 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.

An intrinsic mode function represents a simple oscillatory mode similar to a component in the Fourier-based simple harmonic function, but more general. A systematic way to decompose an earthquake record 𝑋(𝑡), designated as the sifting process of the EMD, is described as follows.

First, identify all the local extrema of the earthquake record. Connect all the local maxima by a cubic spline to produce the upper envelope of the data. Repeat the procedure for the local minima to produce the lower envelope of the data. The upper and lower envelopes should encompass all the data between them. The mean of these two envelopes is designated as 𝑚₁(𝑡), and ℎ₁(𝑡)=𝑋(𝑡)−𝑚₁(𝑡).

In the subsequent process, ℎ₁(𝑡) is treated as the data, then ℎ₁₁(𝑡)=ℎ₁(𝑡)−𝑚₁₁(𝑡), where 𝑚₁₁(𝑡) is the mean of the upper and lower envelopes of ℎ₁(𝑡).

Repeat the previous process until all the conditions in the definition of an intrinsic mode function is achieved. After repeated sifting, ℎ_1𝑘(𝑡) is given byℎ_1𝑘(𝑡)=ℎ_1(𝑘−1)(𝑡)−𝑚_1𝑘(𝑡),(1) where 𝑚_1𝑘(𝑡) is the mean of the upper and lower envelopes of ℎ_1(𝑘−1)(𝑡). ℎ_1𝑘(𝑡) is designated as the first intrinsic mode function 𝑐₁(𝑡) from the earthquake record, that is, 𝑐₁(𝑡)=ℎ_1𝑘(𝑡).

It is noted that the standard deviation SD, which is computed from two consecutive sifting results, is used as the criteria for terminating the sifting process for each intrinsic mode function. SD is defined asSD=^𝑇_𝑡=0ℎ_1(𝑘−1)(𝑡)−ℎ_1𝑘(𝑡)²ℎ^21(𝑘−1)(𝑡),(2) where 𝑇 is the total time length of the time history. When SD≤0.2, the sifting process for each intrinsic mode function is terminated [1].

Typically, 𝑐₁(𝑡) contains the shortest-period component of the record. 𝑐₁(𝑡) is then removed from the original data 𝑋(𝑡) to obtain the residue 𝑟₁(𝑡)=𝑋(𝑡)−𝑐₁(𝑡). The residue 𝑟₁(𝑡), which contains longer-period components, is treated as a new data and subjected to the same sifting process as described above. This procedure can be repeated to obtain all the subsequent 𝑟_𝑗 functions as follows:𝑟_𝑗−1(𝑡)−𝑐_𝑗(𝑡)=𝑟_𝑗(𝑡),𝑗=2,3,…,𝑛.(3) The sifting process can be terminated by either of the following predetermined criteria:(1)either the component 𝑐_𝑛(𝑡) or the residue 𝑟_𝑛(𝑡) becomes so small that it is less than a predetermined value, (2)the residue 𝑟_𝑛(𝑡) becomes a monotonic function.

The original earthquake record is the sum of the intrinsic mode functions 𝑐_𝑗(𝑡) plus the final residue 𝑟_𝑛(𝑡)𝑋(𝑡)=^𝑛_𝑗=1𝑐_𝑗(𝑡)+𝑟_𝑛(𝑡).(4) Hence, the earthquake record is decomposed into 𝑛 intrinsic mode functions 𝑐₁(𝑡), 𝑐₂(𝑡),…,𝑐_𝑛(𝑡) and a residue 𝑟_𝑛(𝑡) that can be either the mean trend or a constant.

Figure 3 shows all eight intrinsic mode functions of the El-Centro time history, recorded in 1940 Imperial Valley Earthquake as shown in Figure 2. It is seen that each component emphasizes a different oscillation mode with different amplitude and frequency content. The first intrinsic mode function has the highest-frequency content, and the frequency content decreases with the increase in intrinsic mode function until the 8th intrinsic mode function.

Figure 2 

Earthquake time history recorded in Imperial Valley Earthquake, 1940.

Figure 3 

Eight intrinsic mode functions of El-Centro time history through empirical mode decomposition.

The final residue 𝑟₈(𝑡) of the El-Centro time history through EMD procedure is shown in Figure 4. This final residue is smaller than a predetermined value, from which no more intrinsic mode function can be extracted. Figure 5 shows the difference between the original El-Centro time history and the reconstructed one given by (4). This small error indicates the completeness of the expansion, and the reconstructed data is numerically identical to the original one.

Figure 4 

Final residue of El-Centro time history through empirical mode decomposition.

Figure 5 

Difference between original El-Centro time history and reconstructed El-Centro time history.

2.2. The Hilbert Spectral Analysis (HSA)

For given data 𝐶(𝑡), the Hilbert transform 𝑌(𝑡) is define as𝑌(𝑡)=1𝜋ℙ^∞−∞𝐶𝑡^𝑡−𝑡^d𝑡^,(5) where ℙ denotes the Cauchy Principal Value. With this definition, 𝐶(𝑡) and 𝑌(𝑡) can be combined to form the analytic signal 𝑍(𝑡), given by𝑍(𝑡)=𝐶(𝑡)+i𝑌(𝑡)=𝑎(𝑡)e^i𝜃(𝑡),(6) where time-dependent amplitude 𝑎(𝑡) and phase 𝜃(𝑡) are𝑎(𝑡)=𝐶²(𝑡)𝑌²(𝑡),𝜃(𝑡)=tan⁻¹𝑌(𝑡)𝐶(𝑡).(7)

From the polar coordinate expression of (6), the instantaneous frequency can be defined as𝜔(𝑡)=d𝜃(𝑡)d𝑡.(8) Applying the Hilbert transform to the 𝑛 intrinsic mode functions of 𝑋(𝑡) in (4), the earthquake record 𝑋(𝑡) can be written as𝑋(𝑡)=ℜ^𝑛_𝑗=1𝑎_𝑗(𝑡)e^{i∫𝜔𝑗(𝑡)d𝑡},(9) where ℜ stands for the real part, 𝑎_𝑗(𝑡) and 𝜔_𝑗(𝑡) are the time-dependent amplitude and instantaneous frequency associated with the 𝑗th intrinsic mode function, respectively. The residue 𝑟_𝑛(𝑡) is not included because of its monotonic property.

Eight acceleration response spectra of intrinsic mode functions of the El-Centro time history are shown in Figure 6. For ease of comparison, the response spectrum of each intrinsic mode function is normalized by its Zero-Period Acceleration (ZPA), which is also known as the Peak Ground Acceleration (PGA) of the earthquake accelerogram. Corresponding to the intrinsic mode functions changing from high frequency to low frequency, the characteristic periods of the response spectra, on which the peak values of the spectral accelerations occur, move from short period to long period as shown in Figure 6.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 6 

Eight Normalized response spectra of intrinsic mode functions from El-Centro time history.

Based on (9), the amplitude and the instantaneous frequency can be represented as functions of time in a three-dimensional plot, in which the amplitude can be contoured or color mapped on the frequency-time plane. This frequency-time distribution of the amplitude 𝐻(𝜔,𝑡) is called the Hilbert amplitude spectrum or simply Hilbert spectrum. The Hilbert energy spectrum can be produced by substituting the amplitude with its squared values. The color mapped Hilbert energy spectrum of the El-Centro earthquake time history is shown in Figure 7. In the Hilbert energy spectrum, the energy increases with the color changing from warm color to cool color. The MATLAB codes for performing the HHT are programmed by authors.

Figure 7 

The Hilbert energy spectrum of El-Centro time history.

3. Generation of Spectrum-Compatible Time Histories Based on Hilbert-Huang Transform

A conventional design spectrum (CDS) is the statistical result of processing actual earthquake records obtained from worldwide ground motion database [9]. A Uniform Hazard Spectrum (UHS) is derived from a Probabilistic Seismic Hazard Analysis (PSHA) [10, 33]. The low-frequency portion of the UHS represents large far-field earthquakes, whereas the high-frequency portion of the UHS represents small near-field earthquakes [7]. Since both CDS and UHS represent the design ground motions by combining the seismic hazard contributions from multiple earthquakes, the selection of seed actual earthquake records can be based on the identification of dominant magnitude-distance pairs that impact the site of interest.

The first step in generating spectrum-compatible time history is to select several actual earthquake records. Small near-field earthquakes and large far-field earthquakes are used to characterize the high-frequency portion and the low-frequency portion of the target design spectrum, respectively. The selection of representative actual earthquake records refers to the results of Seismic Hazard Deaggregation (SHD). SHD can provide the relative hazard contributions of the earthquake sources in terms of distance and magnitude. The selected actual earthquake records, induced from earthquakes with such distances and magnitudes, have relatively the most seismic hazard contributions to the UHS of the site of interest at specified vibration periods. An introduction of the seismic hazard deaggregation can be found in McGuire [10].

The second step is to decompose each of the selected actual earthquake records into several intrinsic mode functions via EMD. The 𝑛 generated intrinsic mode functions from the selected actual earthquake records are treated as the bases to represent a nonstationary earthquake time history. Compared to Fourier transform and wavelet transform, the bases to represent a time series are derived from the actual earthquake records themselves by HHT. Hence, these intrinsic mode functions preserve the information of the nature of the selected earthquakes, which is consistent with the seismic hazard circumstances at the site of interest.

The third step is to transfer each intrinsic mode function from time-domain to frequency-domain. The time-dependent amplitude and the instantaneous frequency of each intrinsic mode function are then generated through HSA. The HHT amplitudes and the instantaneous frequencies of the actual earthquake records are thus scaled to obtain the generated earthquake time historyTH(𝐱,𝑡)=ℜ^𝑛_𝑗=1𝑥_𝑗𝑎_𝑗(𝑡)e^{i𝑥𝑛+𝑗∫𝜔𝑗(𝑡)d𝑡},(10) where 𝑥₁,𝑥₂,…,𝑥_𝑛 are the amplitude scaling parameters and 𝑥_𝑛+1,𝑥_𝑛+2,…,𝑥_2𝑛 are the frequency scaling parameters. 𝐱={𝑥₁,𝑥₂,…,𝑥_2𝑛}^𝑇 is the column vector of the scaling parameters.

The purpose of generating spectrum-compatible earthquake time history is to make the time history rich in all frequencies. By scaling the amplitudes of the intrinsic mode functions, the frequency contents over corresponding frequency ranges can be changed to enrich the time history in all frequencies as required. Since each intrinsic mode function contains a range of frequencies, adjacent intrinsic mode functions of each actual earthquake record may overlap or separate in the frequency domain. Slightly adjusting the instantaneous frequencies of the intrinsic mode functions can ensure that there are no significant gaps in the frequency contents of the intrinsic mode functions.

The fourth step is to construct an optimization model:minimize∶𝑉=^𝑚_𝑖=1𝑓𝑓_𝑖(𝐱),_𝑖(𝐱)=||𝑆^𝑀TH(𝐱,𝑡),𝑇_𝑖−𝑠𝑇_𝑖||subjectedtoconstraints∶𝑥_𝑗𝑠𝑇_𝑖,>0,𝑗=1,2,…,2𝑛,(11) where 𝑆_𝑀[TH(𝐱,𝑡),𝑇_𝑖] is the spectral acceleration of the generated time history TH(𝐱,𝑡) at a specific period 𝑇_𝑖, 𝑠(𝑇_𝑖) is the spectral acceleration of the target design spectrum at 𝑇_𝑖, and 𝑉 is called the objective function.

There are a great number of optimization techniques, such as the Quasi-Newton method and the Nelder-Mead method [40], which can be used to achieve the minimization of the objective function 𝑉 to make the response spectrum of the generated earthquake time history as close to the target design spectrum as possible. For the optimization problem (11) in the proposed generation method, the nonlinear Quasi-Newton algorithm is applied. The detailed description of the Quasi-Newton algorithm can be found in [40].

By solving the optimization problem (11) using a suitable optimization algorithm, a spectrum-compatible earthquake time history, whose response spectrum closely matches the target design spectrum, can be obtained. It contains the desired characteristics of ground motions from actual earthquakes that contribute relatively the most seismic hazard to the site of interest. The generated spectrum-compatible earthquake time history also preserve the nonstationary characteristics of the original earthquake records through the HHT. The schematic block diagram illustrating the proposed generation procedure is shown in Figure 8.

Figure 8 

Flow chart of procedure for generating spectrum-compatible earthquake time histories.

4. Numerical Examples

For illustration purpose, seven numerical examples are presented in this section.

4.1. Target Seismic Design Spectra

The Uniform Hazard Spectra (UHS) for Quebec City, La Malbaie, and eastern North America region (ENA), respectively, are used as the target UHS. The UHS at Quebec City and La Malbaie were derived for ordinary building structures by the Geological Survey of Canada (GSC) [8]. The ENA UHS was derived for nuclear power plants in eastern North America region [7].

The conventional design spectra (CDS) for nuclear energy facilities in Quebec City, La Malbaie, New York City, and Chicago, respectively, are used as the target CDS. The CDS at Quebec City and La Malbaie are from the Geological Survey of Canada (GSC) and Standard CAN3-N289.3-M81 [4]. The CDS at New York City and Chicago are from the U.S. Geological Survey (USGS) and Regulatory Guide 1.60 [6].

In general, nuclear structures (existing plants and new designs) are very stiff with fundamental frequencies in the range of 3–15 Hz [41]. It is adequate for the frequency range of the conventional design spectra for nuclear energy facilities to include from 0.1 Hz to 25 Hz. The selected uniform hazard spectra for nonnuclear building structures also stop at frequency 25 Hz, since these structures are usually more flexible than nuclear structures [42].

The new ENA UHS for nuclear power plants stops at high frequency of 100 Hz because its spectral shape has much larger amplitudes at frequencies between 20 Hz and 100 Hz compared to the conventional design spectra for nuclear energy facilities. It is widely recognized that the new ENA UHS is consistent with the seismic hazard circumstances in eastern North America region and this high-frequency ground motion may cause damaging response stresses in high-frequency sensitive equipments of nuclear power plants [41].

Information of the selected seismic design spectra is listed in Table 1. All seven seismic design spectra are for rock site conditions (average shear wave velocity to a depth of 30 m greater than 750 m/s) and for 5% critical damping.

Table 1 

Seismic design spectra.

4.2. Selection of Seed Earthquake Time Histories

The seed actual earthquake records for each site are selected based on the results of the Seismic Hazard Deaggregation (SHD) [7, 43, 44]. The results of the SHD and the selected earthquake ground motions for all sites are listed in Tables 2 and 3, respectively.

Table 2 

Results of seismic hazard deaggregation.

Table 3 

Selected actual earthquake records.

It should be mentioned that the SHD results are derived from the UHS of the sites of interest. For the conventional design spectra for Quebec City, La Malbaie, New York City, and Chicago, the seed earthquake time histories are still selected based on the SHD results. Hence, the UHS for these four sites were determined first by the probabilistic seismic hazard analysis [8, 45]. The SHD results for these sites are then derived based on their corresponding UHS by the seismic hazard deaggregation [43, 44].

Compared to the target design spectra, the high-frequency portion of the response spectrum of an actual earthquake time history is highly nonsmooth. For Quebec City, La Malbaie, New York City, and Chicago, three earthquakes are selected for each of these four sites. To make the resulting response spectrum closely match the target design spectrum, two small near-field earthquakes are used to characterize the high-frequency plateau portion of the target design spectrum, and one large far-field earthquake is used to describe the low-frequency portion of the target design spectrum.

To characterize the wide range of high frequency of the ENA UHS, three small near-field earthquakes are selected. One large far-field earthquake is used for the low-frequency portion of the ENA UHS.

The actual earthquake accelerograms used in this paper are searched from strong motion databases of the Pacific Earthquake Engineering Research Center (PEER) and the Natural Resources Canada (NRC). Each selected earthquake has two orthogonal horizontal acceleration time histories. All the selected actual earthquake time histories are recorded by free-field instruments or instrument shelters at rock site (average shear wave velocity to a depth of 30 m greater than 750 m/sec). The magnitude-distance pairs of the selected earthquake records are approximately consistent with the results of the SHD and, as a result, generally contribute the most relative seismic hazard to their corresponding sites of interest.

4.3. Spectrum-Compatible Earthquake Time Histories

The response spectra of the generated earthquake time histories and the selected seed time histories for each site are shown in Figures 9, 10, 11, 12, 13, 14, and 27, respectively. The resulting response spectra closely match the target design spectra. The generated spectrum-compatible earthquake time histories and their corresponding seed time histories for each site are shown in Figures 15, 17, 19, 21, 23, 25, and 29, respectively. The generated spectrum-compatible earthquake time histories generally preserve the temporal characteristics of the actual seed time histories.

Figure 9 

Quebec City: target uniform hazard spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 10 

La Malbaie: target uniform hazard spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 11 

Quebec City: target N289.3 design spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 12 

La Malbaie: target N289.3 design spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 13 

New York City: target RG 1.60 design spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 14 

Chicago: target RG 1.60 design spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 15 

Seed and generated earthquake time histories for UHS at Quebec City.

In matching the high-frequency portion of the response spectrum of the generated time history to the target spectrum, the high-frequency intrinsic mode functions, that is, the first one or two components, are the main contributors. However, because these high-frequency intrinsic mode functions cover a relatively wide range of frequency, it is very difficult to just match the high-frequency portion of the generated and target spectra without significantly affecting the match in the mid- and low-frequency ranges. Since the optimization model (11) can be regarded as a comprise in selecting the scaling parameters to reach a minimal overall relative difference between the generated and the target spectra, it is not hard to see that the optimization model tends to sacrifices the match in the high-frequency portion.

As shown in Figure 27, the difference between the resulting spectrum and the target spectrum is relatively large over the frequency range from 60 Hz to 100 Hz. Further studies may be needed to improve the empirical mode decomposition and to decompose the signal into more components in the high-frequency range. The results will also be improved if more earthquake ground motions with higher-frequency contents are available.

The color mapped Hilbert energy spectra of the generated spectrum-compatible earthquake time histories for each numerical example are shown in Figures 16, 18, 20, 22, 24, 26, and 28, respectively. It can be seen that the generated time histories have sufficient energy over the entire frequency range. More energy can also be observed at the beginning of the time histories through this frequency-time distribution of energy, which is consistent with the nonstationary characteristics of the generated time histories.

Figure 16 

The Hilbert energy spectrum of generated earthquake time history for UHS at Quebec City.

Figure 17 

Seed and generated earthquake time histories for UHS at La Malbaie.

Figure 18 

The Hilbert energy spectrum of generated time history for UHS at La Malbaie.

Figure 19 

Seed and generated earthquake time histories for CDS at Quebec City.

Figure 20 

The Hilbert energy spectrum of generated time history for CDS at Quebec City.

Figure 21 

Seed and generated earthquake time histories for CDS at La Malbaie.

Figure 22 

The Hilbert energy spectrum of generated time history for CDS at La Malbaie.

Figure 23 

Seed and generated earthquake time histories for CDS at New York City.

Figure 24 

The Hilbert energy spectrum of generated time history for CDS at New York City.

Figure 25 

Seed and generated earthquake time histories for CDS at Chicago.

Figure 26 

The Hilbert energy spectrum of generated time history for CDS at Chicago.

Figure 27 

Eastern North America: target uniform hazard spectrum, response spectrum of generated earthquake time history, and response spectra of seed actual earthquake records for 5% damping and rock site.

Figure 28 

The Hilbert energy spectrum of generated time history for ENA UHS.

Figure 29 

Seed and generated earthquake time histories for ENA UHS.

4.3.1. Selection of Intrinsic Mode Functions as Basis for Time History Generation

As shown in the numerical examples, the seed earthquake time histories are decomposed into approximate 50 intrinsic mode functions as the basis for generating a spectrum-compatible earthquake time history. However, there are overlaps between many response spectra of these intrinsic mode functions, that is, the intrinsic mode functions occupy similar frequency contents. These overlapping intrinsic mode functions may be redundant for the purpose of generating spectrum-compatible earthquake time histories.

The purpose of generating spectrum-compatible earthquake time history is to make the time history rich in all frequencies. Hence, a small number of intrinsic mode functions, whose narrow-band response spectra cover from low frequency to high frequency without significant gaps between each other, can be selected as the basis for time history generation. To include the crucial information of all seed earthquakes, at least one high-frequency intrinsic mode function is selected from each seed time history, induced by small near-field earthquake, as the high-frequency basis and at least one low-frequency intrinsic mode function is chosen from each seed time history, induced by large far-field earthquake, as the low-frequency basis.

Figure 30 shows 21 narrow-band normalized response spectra of the intrinsic mode functions, which are selected from a pool of 59 intrinsic mode functions for Quebec City following the criteria described in previous paragraph. As shown in Figure 31, a satisfactory resulting spectrum is achieved based on the selected 21 intrinsic mode functions. This procedure can significantly reduce the numbers of intrinsic mode functions with overlapping frequencies, hence the dimension of the optimization problem, and speed up the convergence of the computation.

Figure 30 

21 normalized response spectra of intrinsic mode functions for Quebec City.

Figure 31 

Comparison of response spectra based on full intrinsic mode functions and selected intrinsic mode functions.

5. Conclusion and Discussion

A new method, which can generate an earthquake time history compatible with the target seismic design spectrum (UHS or CDS) based on Hilbert-Huang transform of multiple actual earthquake records and optimization technique, is proposed. The generation procedure is suitable for any type of seismic design spectra. The generated earthquake time histories preserve the nonstationary characteristics of the seed actual earthquake records.

From the numerical examples, the following observations can be made.(1)The success of the proposed generation procedure in any practical situation will depend much on the availability and the selection of strong ground motion records. Although the results of seismic hazard deaggregation in terms of magnitude-distance pairs provide much useful information for earthquake record selection, strong ground motions, which are recorded around the sites of interest and reflect the most realistic seismic hazard of the sites, are still desired. This difficulty will be alleviated with the increase of the number of the available earthquake records over time.(2) A classical optimization algorithm, Quasi-Newton method, is used to minimize the objective function in this paper. Since the objective function of the proposed method may be highly nonlinear or non-smooth, it is challenging to find global minima by the traditional optimization algorithms. The initial values of the optimization parameters may also influence the results of the generation computation. Further studies may be needed to find the most suitable optimization algorithm for the proposed generation method.

Acknowledgments

The research for this paper was supported, in part, by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the University Network of Excellence in Nuclear Engineering (UNENE). The authors are grateful to the referees for the constructive comments which helped to improve the paper.