# Combination of Spectral Representation and Wavelet Packets for Generating Long-Period Ground Motions

Authors: Minghui Dai; Yingmin Li
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Far-field long-period ground motions (hereafter long-period ground motions) featuring low-frequency components are responsible for the resonant responses of high-rise buildings. In this context, it is beneficial to assess the dynamic performance of these buildings under long-period ground motions with the aid of time history analysis. This paper proposes a method for generating long-period motions by combining long-period components synthesized by spectral representation with high-frequency components simulated by wavelet packets. Later-arriving long-period surface waves (LALP surface waves), which are determined on the grounds of phase dispersion, represent the main long-period properties in sense of velocity spectrum at longer periods of interest. An analytical expression for power spectrum density is employed to capture the narrowband properties of LALP velocity surface waves. Meanwhile, modification of the Gaussian random process is performed in time and frequency domains to attain a modulated initial seed motion, which shows the variability of the targeted ground motion. A simulation of high-frequency components is accomplished by means of iteration, in which wavelet coefficients of the modulated seed motion are adjusted to match the targeted response spectrum and cumulative energy plot. Furthermore, comparisons between an ensemble of realizations and target motions demonstrate the feasibility of the proposed method to generate long-period simulations sharing similar properties to target motions.


---

## Body

## 1. Introduction

Far-field long-period ground motions, which consist of dominant low-frequency components, are likely to enable high-rise buildings to excite resonance responses lasting several minutes [1]. It has been reported that damage to either structural components or non-structural components under long-period ground motions gives rise to unsafe occupancy and considerable maintenance expenses [2–5]. Thus, it is necessary to evaluate the performance of such structures subjected to long-period ground motions by means of time history analysis. However, networks with the capability to record a high signal-to-noise ratio over a long-period range are insufficient. In addition, large-magnitude events in regions with active seismicity have low probabilities of occurrence, accounting for the low availability of long-period ground motions. In response to this challenge, artificial earthquake ground motions are essential to enrich electronic libraries of long-period ground motions.The generation of ground motions compatible with a prescribed spectrum is attractive to engineers, since it is considered an extension of stochastic simulations. In essence, this process modifies the amplitudes of Fourier components in the frequency domain to attain a realization which closely matches a given response spectrum [6, 7]. However, Fourier-based modification without involving time-varying frequency content fails to generate nonstationary of acceleration time series. To circumvent these disadvantages, wavelet analysis characterized by time-frequency discretization is introduced to generate artificial ground motions. For example, Spanos et al. [8, 9] employed harmonic wavelets to synthesize artificial ground motions compatible with the design response spectra mandated by either Chinese or European aseismic regulation through iteratively modifying wavelet coefficients. As referred from the wavelet-packet-based stochastic model suggested by Yamamoto and Baker [10], Huang and Wang [11] accomplished the generation of artificial accelerograms conditional on agreement with the response spectrum in frequency domain and cumulative energy in time domain. However, a limitation of these methods is that the extent of compatibility with a target spectrum decreases as the periods shift towards a long-period range.According to the seismological mechanism, long-period ground motions are characterized by LALP surface waves, which are converted from body waves at the margins of either basins or plains [12–15]. Because of these phenomena, the identification of long-period ground motions with broadband long-period characteristics is developed considering that energy from LALP components is predominantly relative to that of high frequencies [16]. However, it has been suggested that LALP surface waves in the form of velocity traces are equivalent to the summation of sinusoids with narrow-frequency components in special cases [17], as amplifications to the amplitudes of LALP surface waves occur at specific long periods. Thus, long-period ground motions share similar properties with pulse-like motions in view of long-period characteristics in frequency domain. Moustafa and Takewaki [18] adopted a function shaped-like Kronecker delta to describe frequency components of the pulse-like traces with one or few predominant frequencies on the basis of the stationary process [19–21]. In contrast, finite difference methods for long-period motions at the expense of computational resources and seismological information [2, 22, 23] are not applicable for engineering communities. In this regard, the spectral representation for long-period components of long-period motions can be useful to satisfy demands in engineering practices.The objective of this paper is to generate long-period ground motions through spectral representation for generating long-period components and wavelet packets for simulating high-frequency components. To achieve this aim, a time segment consisting of LALP components is truncated from acceleration time series on the basis of phase dispersion. Then, LALP surface velocity traces are simulated on the assumption of a stationary narrowband process. On the other hand, an initial seed motion is generated as a result of the modification of the Gaussian random process in time and frequency domains. Wavelet coefficients of the initial seed motion are adjusted by an iterative procedure to generate high-frequency components which are compatible with the target spectrum and cumulative energy plot. Finally, comparisons of simulations with target motions demonstrate that the proposed method has the potential to generate ground motions sharing similar long-period characteristics with target motions.


## 2. Long-Period Components

### 2.1. Determination of Long-Period Components

Similar to artificial generation for pulse-like ground motions [20, 24–26], determination of which fraction of a time history captures the gross long-period characteristics is a key issue in the separation of surface waves from body waves. This facilitates subsequent numerical modelling. Actually, there is no principle to conduct separation of long-period surface waves from its time history as accurately as possible, since the components which are dependent on frequencies always overlap in most cases. However, with reference to articles [14, 15], variations in the speed of wave trains in propagation media account for phase dispersion, which are referred to as “envelope delays” [13]. Specifically, the speeds of long-period surface waves increase with periods, as opposed to those of body waves, which decrease with periods.To illustrate phase dispersion in nature, we estimate the arrival times of phases which are dependent on periods by means of an array of narrowband Gaussian filters [27, 28]. Two acceleration time histories for CHB002EW and TKY023EW with epicentral distances 367 km and 377 km during the 2011 Tohuku, M9 earthquake, are considered representatives of long-period ground motions. Thus, we use them as examples to illustrate the proposed method in the following sections. As shown in Figure 1, phases with periods of 2, 5, 7, and 10 sec are individually extracted from the two records by using Gaussian filters. The variations of phases are shown in the form of arrival times corresponding to envelope peaks. It is seen that the arrival times of the long-period phases tend to arrive later with periods decreasing, while the 2 sec phase arrives at nearly the same time as that for PGA. This general trend coincides with the results from the analysis of long-period motions in the Po Plain of Italy [27].Figure 1
Phases dependent on periods from 2, 5, 7, and 10 sec are attained from the acceleration traces of (a) CHB002EW and (b) TKY023EW, by using a set of Gaussian filters.
(a)(b)The principle of determining long-period components is defined herein. It is conditioned that the truncated time series should match well with its original record in terms of long-period characteristics. Following this line, the steps are as follows.First, the maximum intervaltmc of zero-level crossings for acceleration time series is calculated such that the period p=2tmc of the earliest-arrival long-period component is determined based on phase dispersion illustrated above; second, the phase with the period of p is extracted from its acceleration time series by using the designed narrowband Gaussian filter; finally, a reference time point tl which is used to determine the LALP surface waves should be satisfied with the condition, tl∈tpga,tp. At this point, we could possibly search for an appropriate tl, after which the time segment consisting of sufficient surface waves enables the velocity response spectrum of the time segment to show an agreement with that of its original record around dominant long periods. Note that tp is the first wave pack of the filtered phase shown in Figure 2 and tpga is the arrival time for PGA.Figure 2
Illustration of determination of long-period components for (a) CHB002EW and (b) TKY023EW.
(a)(b)The detailed explanation of the condition fortl is that the tp is often treated as an arrival time for the filtered phase with a period of p [15, 27]; however, due to phase dispersion and scatting effects, it is expected that the truncated time segment contains sufficient surface waves energy at the dominant period as much as possible, thus resulting in the upper bound for tp; meanwhile, it is reasonable that later-arriving surface waves are converted from body waves [14, 15], thus resulting in the lower bound for tpga.As a result, the examples are shown in Figure2. The tmc for CHB002EW is 3.04 sec, and the corresponding p=2×3.04=6.08sec; a Gaussian filter with a centre period of p is employed to extract the phase shown in Figure 2(a); then, as shown in Figure 2(b), a time segment mainly consisting of long-period surface waves is truncated from its acceleration time series with reference to the tl, and this is conditioned by its velocity spectrum matching well with that of the original record (illustrated in Figure 3). Similarly, the process for TKY023EW is illustrated in Figure 2(b). The identified parameters of tl for these records are listed in Table 1.Figure 3
Comparison of truncated time series with corresponding original records in terms of velocity response spectrum. (a) CHB002EW and (b) TKY023EW.
(a)(b)Table 1
Parameters of determination and PSD for truncated long-period waves.
Recordtl (sec)αβpqCHB002EW150.58275.410.6718.1411.47TKY023EW138.86107.070.3512.019.96SIT003EW157.339999.034.628.963.59It is observed from Figure3 that the velocity response spectra of the truncated waves at long predominant periods show an agreement with those of the corresponding original records. In contrast, there are obvious deviations at short periods. In this regard, the truncated waves are capable of capturing the gross long-period characteristics of the original records.


### 2.2. Representation of Long-Period Components

As illustrated above, truncated waves with gross long-period properties are characterized by narrowband properties. Because of this, the simulation of long-period components consisting of LALP surface waves is generated on the basis of the stationary power spectrum density function (PSDF). Accordingly, an exponential form for modelling the PSDF of LALP surface waves is written as in equation (1), as is derived from the works regarding long-period waves:(1)Sω=αωpexp−β1ωq,where α is the scaling constant to modulate the amplitudes; ω denotes the angular frequency (rad/s); β, p, and q are coefficients related to the frequency domain.The derivative of the form with respect to angular frequencyω is given as follows:(2)dsωdω=α−pωp+1exp−βωq+αωpexp−βωq−β−qωq+1.It is easy to know that the maximum value forsω is attained when dsω/dω equals zero, and the predominant angular frequency is determined. Then, equation (2) can be rewritten as(3)αpωp+1exp−βωq=αpωpexp−βωqβqωq+1.Solving this equation, the predominant frequency is shown as follows:(4)fd=1q2πβqp.For equation (4), coefficients β, p, and q collectively control the predominant frequency, and the variability of the PSDF is available by modulating the relative proportion among three coefficients. As shown in Figure 4, the normalized PSDs with various parameters of β, p, and q demonstrate the ability of the function to capture narrow bandwidth where the greatest portion of energy is concentrated. At the same time, it is found that the bandwidth becomes wider when the predominant frequency shifts towards higher frequencies.Figure 4
Normalized power spectrum density with various parameters ofβ, ρ, and q.An actual PSD is estimated from LALP velocity surface waves using the modified periodogram of the Welch technique with Gaussian smoothing. Identification of the parameters is performed by fitting the analytical function equation (1) to the actual PSD in a nonlinear least-square sense. Thus, the fitting formula is presented as follows:(5)α^,β^,p^,q^=argminα,β,p,q∑i=1nPSDω−sω;α,β,p,q2.As a consequence, the parameters ofα, β, p, and q for the records CHB002EW and TKY023EW are tabulated in Table 1. Comparisons of the resulting PSDs with the actual PSDs are illustrated in Figures 5(a) and 5(b), respectively. In general, it is clear that the analytical functions perform well in modelling actual PSDs considering the narrowband properties.Figure 5
PSD for truncated LALP velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)For the temporal domain, we use Hilbert transform to obtain instantaneous amplitudes of the velocity time series. The envelopeelt for LALP velocity surface waves in time series is calculated as follows:(6)elt=vt2+Hvt2,in which Hvt is the Hilbert transform, vt is the velocity time series, and t is the time instant.As the instantaneous amplitudes are shown in Figure6, the multiple local maxima increase the difficulty in selecting an appropriate envelope model. To reduce deviations, the third level envelope of LALP velocity surface waves is attained by using the multimodal point method [29]. Then, the double-exponent function in the following equation is used to reflect the change in amplitudes with time:(7)e¯lt=c1expλ1t+c2expλ2t,where c1, c2, λ1, and λ2 are the fitted coefficients and e¯l denotes the fitted envelope.Figure 6
Instantaneous amplitudes for truncated velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)Correspondingly, the fitted envelopese¯l presented in Figure 6 are attained as a result of the fitting process, and fitted parameters are listed in Table 2.Table 2
Envelope parameters for truncated long-period velocity surface waves.
Recordc1λ1c2λ2CHB002EW5.83−0.00313.21−0.0030TKY023EW12.80−0.1112.90−0.0062SIT003EW13.88−0.01110.55−0.011It is observed that the fitted envelopese¯l are capable of capturing the trend of amplitudes with time evolution, and thus, the function can adapt to the different shapes by modulating the parameters.Considering the characteristics of PSD for LALP velocity time series, the time series is treated as a summation of cosine waves with narrowband frequencies in the sense of spectral representation, as defined in the following equation:(8)xt=2e¯l∑n=1Nsω;α,β,p,q2πΔf1/2cos2πfn+ϕn,where e¯l is the fitted envelope of the LALP velocity time series; sω;α,β,p,q denotes the power spectral density function with parameters of α,β,p,q; fn∈fl,fu is the frequency band; fl and fu are the lower and upper cutoff frequencies of the PSDF, respectively; the frequency step is Δf=fu−fl/N; fu is dependent on sampling rate of the time series; and fl is 0.01 Hz; each discrete frequency is fn=fl+nΔf; ∅n is a set of random phase angles distributed uniformly over the range of 0,2π; and N is the number of time steps.Note that the acceleration time history should be performed in baseline correction through a fourth-order Butterworth filter with bandwidth ranging from 0.1 to 1 Hz. This enables the removal of spurious drifts and the retention of the long-period properties. Consequently, it is seen from Figure7 that the realizations present stochastic properties in time domain, but predominant frequencies controlling the LALP surface waves cannot vary over time.Figure 7
Simulated long-period surface waves in the form of velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)Comparisons of the simulations with truncated long-period surface waves in terms of 5% damping spectra are shown in Figure8. It is found that the simulated spectra share a similar predominant period with the truncated waves. This is consistent with the results in Figure 5. Meanwhile, the truncated waves almost fluctuate between the median plus or minus one standard deviation ±σ. Because of this, the truncated waves are possibly treated as one sample of the realizations on the basis of stochastic properties.Figure 8
Response spectra for an ensemble of realizations of truncated long-period surface waves from two records of CHB002EW and TKY023EW (solid blue line). Of 100 simulated ground motions (thin gray lines), their median (thick red line) and median plus and minus one standard deviation±σ (dashed lines) are involved; acceleration spectra on the left side; velocity spectra on the right side. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 2.1. Determination of Long-Period Components

Similar to artificial generation for pulse-like ground motions [20, 24–26], determination of which fraction of a time history captures the gross long-period characteristics is a key issue in the separation of surface waves from body waves. This facilitates subsequent numerical modelling. Actually, there is no principle to conduct separation of long-period surface waves from its time history as accurately as possible, since the components which are dependent on frequencies always overlap in most cases. However, with reference to articles [14, 15], variations in the speed of wave trains in propagation media account for phase dispersion, which are referred to as “envelope delays” [13]. Specifically, the speeds of long-period surface waves increase with periods, as opposed to those of body waves, which decrease with periods.To illustrate phase dispersion in nature, we estimate the arrival times of phases which are dependent on periods by means of an array of narrowband Gaussian filters [27, 28]. Two acceleration time histories for CHB002EW and TKY023EW with epicentral distances 367 km and 377 km during the 2011 Tohuku, M9 earthquake, are considered representatives of long-period ground motions. Thus, we use them as examples to illustrate the proposed method in the following sections. As shown in Figure 1, phases with periods of 2, 5, 7, and 10 sec are individually extracted from the two records by using Gaussian filters. The variations of phases are shown in the form of arrival times corresponding to envelope peaks. It is seen that the arrival times of the long-period phases tend to arrive later with periods decreasing, while the 2 sec phase arrives at nearly the same time as that for PGA. This general trend coincides with the results from the analysis of long-period motions in the Po Plain of Italy [27].Figure 1
Phases dependent on periods from 2, 5, 7, and 10 sec are attained from the acceleration traces of (a) CHB002EW and (b) TKY023EW, by using a set of Gaussian filters.
(a)(b)The principle of determining long-period components is defined herein. It is conditioned that the truncated time series should match well with its original record in terms of long-period characteristics. Following this line, the steps are as follows.First, the maximum intervaltmc of zero-level crossings for acceleration time series is calculated such that the period p=2tmc of the earliest-arrival long-period component is determined based on phase dispersion illustrated above; second, the phase with the period of p is extracted from its acceleration time series by using the designed narrowband Gaussian filter; finally, a reference time point tl which is used to determine the LALP surface waves should be satisfied with the condition, tl∈tpga,tp. At this point, we could possibly search for an appropriate tl, after which the time segment consisting of sufficient surface waves enables the velocity response spectrum of the time segment to show an agreement with that of its original record around dominant long periods. Note that tp is the first wave pack of the filtered phase shown in Figure 2 and tpga is the arrival time for PGA.Figure 2
Illustration of determination of long-period components for (a) CHB002EW and (b) TKY023EW.
(a)(b)The detailed explanation of the condition fortl is that the tp is often treated as an arrival time for the filtered phase with a period of p [15, 27]; however, due to phase dispersion and scatting effects, it is expected that the truncated time segment contains sufficient surface waves energy at the dominant period as much as possible, thus resulting in the upper bound for tp; meanwhile, it is reasonable that later-arriving surface waves are converted from body waves [14, 15], thus resulting in the lower bound for tpga.As a result, the examples are shown in Figure2. The tmc for CHB002EW is 3.04 sec, and the corresponding p=2×3.04=6.08sec; a Gaussian filter with a centre period of p is employed to extract the phase shown in Figure 2(a); then, as shown in Figure 2(b), a time segment mainly consisting of long-period surface waves is truncated from its acceleration time series with reference to the tl, and this is conditioned by its velocity spectrum matching well with that of the original record (illustrated in Figure 3). Similarly, the process for TKY023EW is illustrated in Figure 2(b). The identified parameters of tl for these records are listed in Table 1.Figure 3
Comparison of truncated time series with corresponding original records in terms of velocity response spectrum. (a) CHB002EW and (b) TKY023EW.
(a)(b)Table 1
Parameters of determination and PSD for truncated long-period waves.
Recordtl (sec)αβpqCHB002EW150.58275.410.6718.1411.47TKY023EW138.86107.070.3512.019.96SIT003EW157.339999.034.628.963.59It is observed from Figure3 that the velocity response spectra of the truncated waves at long predominant periods show an agreement with those of the corresponding original records. In contrast, there are obvious deviations at short periods. In this regard, the truncated waves are capable of capturing the gross long-period characteristics of the original records.


## 2.2. Representation of Long-Period Components

As illustrated above, truncated waves with gross long-period properties are characterized by narrowband properties. Because of this, the simulation of long-period components consisting of LALP surface waves is generated on the basis of the stationary power spectrum density function (PSDF). Accordingly, an exponential form for modelling the PSDF of LALP surface waves is written as in equation (1), as is derived from the works regarding long-period waves:(1)Sω=αωpexp−β1ωq,where α is the scaling constant to modulate the amplitudes; ω denotes the angular frequency (rad/s); β, p, and q are coefficients related to the frequency domain.The derivative of the form with respect to angular frequencyω is given as follows:(2)dsωdω=α−pωp+1exp−βωq+αωpexp−βωq−β−qωq+1.It is easy to know that the maximum value forsω is attained when dsω/dω equals zero, and the predominant angular frequency is determined. Then, equation (2) can be rewritten as(3)αpωp+1exp−βωq=αpωpexp−βωqβqωq+1.Solving this equation, the predominant frequency is shown as follows:(4)fd=1q2πβqp.For equation (4), coefficients β, p, and q collectively control the predominant frequency, and the variability of the PSDF is available by modulating the relative proportion among three coefficients. As shown in Figure 4, the normalized PSDs with various parameters of β, p, and q demonstrate the ability of the function to capture narrow bandwidth where the greatest portion of energy is concentrated. At the same time, it is found that the bandwidth becomes wider when the predominant frequency shifts towards higher frequencies.Figure 4
Normalized power spectrum density with various parameters ofβ, ρ, and q.An actual PSD is estimated from LALP velocity surface waves using the modified periodogram of the Welch technique with Gaussian smoothing. Identification of the parameters is performed by fitting the analytical function equation (1) to the actual PSD in a nonlinear least-square sense. Thus, the fitting formula is presented as follows:(5)α^,β^,p^,q^=argminα,β,p,q∑i=1nPSDω−sω;α,β,p,q2.As a consequence, the parameters ofα, β, p, and q for the records CHB002EW and TKY023EW are tabulated in Table 1. Comparisons of the resulting PSDs with the actual PSDs are illustrated in Figures 5(a) and 5(b), respectively. In general, it is clear that the analytical functions perform well in modelling actual PSDs considering the narrowband properties.Figure 5
PSD for truncated LALP velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)For the temporal domain, we use Hilbert transform to obtain instantaneous amplitudes of the velocity time series. The envelopeelt for LALP velocity surface waves in time series is calculated as follows:(6)elt=vt2+Hvt2,in which Hvt is the Hilbert transform, vt is the velocity time series, and t is the time instant.As the instantaneous amplitudes are shown in Figure6, the multiple local maxima increase the difficulty in selecting an appropriate envelope model. To reduce deviations, the third level envelope of LALP velocity surface waves is attained by using the multimodal point method [29]. Then, the double-exponent function in the following equation is used to reflect the change in amplitudes with time:(7)e¯lt=c1expλ1t+c2expλ2t,where c1, c2, λ1, and λ2 are the fitted coefficients and e¯l denotes the fitted envelope.Figure 6
Instantaneous amplitudes for truncated velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)Correspondingly, the fitted envelopese¯l presented in Figure 6 are attained as a result of the fitting process, and fitted parameters are listed in Table 2.Table 2
Envelope parameters for truncated long-period velocity surface waves.
Recordc1λ1c2λ2CHB002EW5.83−0.00313.21−0.0030TKY023EW12.80−0.1112.90−0.0062SIT003EW13.88−0.01110.55−0.011It is observed that the fitted envelopese¯l are capable of capturing the trend of amplitudes with time evolution, and thus, the function can adapt to the different shapes by modulating the parameters.Considering the characteristics of PSD for LALP velocity time series, the time series is treated as a summation of cosine waves with narrowband frequencies in the sense of spectral representation, as defined in the following equation:(8)xt=2e¯l∑n=1Nsω;α,β,p,q2πΔf1/2cos2πfn+ϕn,where e¯l is the fitted envelope of the LALP velocity time series; sω;α,β,p,q denotes the power spectral density function with parameters of α,β,p,q; fn∈fl,fu is the frequency band; fl and fu are the lower and upper cutoff frequencies of the PSDF, respectively; the frequency step is Δf=fu−fl/N; fu is dependent on sampling rate of the time series; and fl is 0.01 Hz; each discrete frequency is fn=fl+nΔf; ∅n is a set of random phase angles distributed uniformly over the range of 0,2π; and N is the number of time steps.Note that the acceleration time history should be performed in baseline correction through a fourth-order Butterworth filter with bandwidth ranging from 0.1 to 1 Hz. This enables the removal of spurious drifts and the retention of the long-period properties. Consequently, it is seen from Figure7 that the realizations present stochastic properties in time domain, but predominant frequencies controlling the LALP surface waves cannot vary over time.Figure 7
Simulated long-period surface waves in the form of velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)Comparisons of the simulations with truncated long-period surface waves in terms of 5% damping spectra are shown in Figure8. It is found that the simulated spectra share a similar predominant period with the truncated waves. This is consistent with the results in Figure 5. Meanwhile, the truncated waves almost fluctuate between the median plus or minus one standard deviation ±σ. Because of this, the truncated waves are possibly treated as one sample of the realizations on the basis of stochastic properties.Figure 8
Response spectra for an ensemble of realizations of truncated long-period surface waves from two records of CHB002EW and TKY023EW (solid blue line). Of 100 simulated ground motions (thin gray lines), their median (thick red line) and median plus and minus one standard deviation±σ (dashed lines) are involved; acceleration spectra on the left side; velocity spectra on the right side. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 3. Wavelet Packets

Wavelet packets analysis is applied to the decomposition of a time seriesxt into a set of wavelet coefficients cj,ki. This means that energy is localized in both time and frequency domains. The basic expression in the form of convolution is defined as follows:(9)cj,ki=∫−∞∞xtψj,kitdt,where ψj,kit is the ith set of the modified mother wavelet, which is determined by both the jth scale parameter that dilates in frequency domain and the kth translation parameter that shifts through the time axis. The Meyer wavelet is chosen as the mother wavelet herein, as proposed in similar work [10].Conversely, a time seriesxt of 2N reconstructed from the wavelet coefficients on the basis of inverse wavelet packets transform is given as follows:(10)xt=∑i=12j∑k=12N−jcj,kiψj,kit.On the other hand, the trade-off between frequency and time resolution should be taken into account because the resolution is in accordance with the Heisenberg uncertainty principle. This means that the high resolution of frequency information comes at the expense of low resolution in temporal domain. Accordingly, the time resolutiondtw and frequency resolution dfw are expressed in equations (11) and (12), respectively, after the jth decomposition step. The product of dtw and dfw is a constant in equation (13):(11)dtw=2jdt,(12)dfw=fN2j,(13)dtw⋅dfw=12,where j is the scale parameter, dt is the time interval of a time series, and fN is the Nyquist frequency determined by the sampling rate of the time series. Compared with the time resolution, detailed information related to frequency domain is of interest in engineering practice [11]. This facilitates simulations in sense of compatibility with a target spectrum.To illustrate above, the time history for CHB002EW with the sampling intervaldt=0.01 is decomposed into 512 frequency rows at the 9th decomposition depth. Correspondingly, the frequency interval descends to 0.0976 Hz, while the corresponding time interval between the centres of adjacent wavelet packets enlarges to dtw=29∗0.01=5.12sec.As a result, Figure9 illustrates the squared wavelet coefficients cj,ki2 related to the time-frequency distribution. For frequency domain, the majority of coefficients at frequencies ranging from 4 to 7 Hz correspond to the peak values of the acceleration response spectrum in Figure 9(c); for temporal domain, the majority of coefficients occurring at approximately 100 sec correspond to a strong shocking phase during which the cumulative energy curve shows a sharp increase in Figure 9(d). This suggests that the wavelet packets transform has the potential to capture the characteristics pertaining to intensity and spectrum, each of which is well correlated with the time and frequency domains, respectively.Figure 9
Relationships pertaining to time and frequency domains are illustrated through wavelet packets decomposition for the record CHB002EW: (a) acceleration time history; (b) distribution of squared wavelet coefficients; (c) acceleration response spectrum; (d) cumulative energy curve.However, it should be kept in mind that the capability of wavelet packets to detect the components at longer periods cannot be the same as that for high-frequency components, because of shifting and scaling based on a power of 2. To illustrate this, the resulting frequency intervalf∈0.0977,0.1953 leads to a jump from 5.12 to 10.23 sec in the form of period. This means that wavelet analysis is not suitable for discerning long-period components in frequency domain. Instead, it specializes in localization of energy at high frequencies.At this point, high-frequency components are generated for the purpose of compatibility with a given spectrum in frequency domain and cumulative energy in time domain, as referred from the method proposed by Huang and Wang [11].


## 4. Formation of High-Frequency Components Using Wavelet Packets

### 4.1. Generation of Initial Seed Motion

Prior to the generation of high-frequency components, an initial seed motion should be created, because variability of the initial seed motion could influence nonstationary properties of desired motions [11, 30]. Thus, the initial seed motions for high-frequency components are created as a result of the modulation of the Gaussian random process with respect to time and frequency domains.


#### 4.1.1. Modulation in Time Domain

By visual observation of the acceleration waveforms shown in Figure2, a piecewise function is employed to modulate the Gaussian random process in time domain, and its mathematical expression is written as follows:(14)eht=0,t≤t0,t−t0tm−t0γ,t0<t≤tm,exp−ηt−tm,t>tm,where tm denotes the time corresponding to PGA; γ and η are parameters controlling the ascent and descent stage of the time segment before and after tm, respectively; t0 indicates the onset of a ground motion. Correspondingly, as shown in Figure 10, a set of normalized envelopes regarding various parameters of tm, γ, and η demonstrates the variability of waveform shapes.Figure 10
Piecewise function with various parameters oftm, γ, and η.Furthermore, a normalized cumulative energy model available from integration of the piecewise function with respect to time is given as follows:(15)ht=∫eh2tdt∫0Teh2t2dt=12η+tm−t0/2γ+10,t≤t0,tm−t02γ+1t−t0tm−t02γ+1,t0<t≤tm,−12ηexp−2ηt−tm+12η+tm−t02γ+1,t>tm.Thus, the parameters of the envelope function are attained by minimizing the integrated differences between the normalized actual energy curve and piecewise energy model, that is,(16)t^0,γ^,η^=argmint0,γ,η∫0t∫0th2tdt−∫0tha2tdt2dt,(17)hat=∫0tx2tdt∫0Dx2tdt,where xt is the acceleration time series, D is the duration of an acceleration time series, and hat denotes the normalized cumulative energy. The fitted parameters are listed in Table 3.Table 3
Parameters controlling envelope of motions.
Recordt0γηCHB002EW0.0114.450.061TKY023EW22.611.960.030SIT003EW0.0482.680.019It is found from Figure11 that the fitted curves are in close agreement with their actual curves at the strong shaking phases, during which the curves appear to be a vertical line. Meanwhile, discrepancies occur at approximately 140 sec, after which LALP surface waves arrive. Although the amplitudes of LALP surface waves are relatively small, in distinct contrast to those of body waves, they have a low attenuation rate. This results from the long-period properties of surface waves.Figure 11
Normalized cumulative energy curves for records and modulated Gaussian process. (a) CHB002EW and (b) TKY023EW.
(a)(b)


#### 4.1.2. Modification in Frequency Domain

Wavelet coefficients in time-frequency domain comply with a bivariate lognormal distribution, which is used to present the variability of ground motions [10]. To simplify this, a set of shaped-like lognormal functions are used to describe the power spectrum density for each column (time instant). Thus, a set of squared wavelet coefficients derived from the modulated Gaussian process are modified by multiplying each fitted function to simulate the variability of targeted motions.The shaped-like lognormal function of thekth column with respect to frequency f is given as follows:(18)Lkfuk,σk=1fσk2πexp−lnf−uk22σk2,where the parameters of μk and σk for the kth column are attained by fitting the assigned function to squared coefficients at the kth column in the sense of power energy and Lk represents the power spectrum at the kth column.As a result, normalized squared coefficients for specific columns of CHB002EW and TKY023EW are shown in Figure12, together with the fitted functions. The fitted models generally reflect the distribution of the power spectra because their energy tends to concentrate at low frequencies as time columns evolve.Figure 12
Comparison of the fitted models with actual records in terms of normalized power spectrum at specific columns. (a) CHB002EW and (b) TKY023EW.
(a)(b)To establish a parametric relationship between the identified parameters and the sequence of columns (time instants), a generic expression built on the superposition of Gaussian functions is defined as follows:(19)λkk=∑isiexp−k−εiψi2,where subscript i indicates the ith Gaussian function; si is the scaling factor; εi indicates the location for peak value; ψi denotes the width of “bell” Gaussian function; and λk represents the variables of μk and σk in the kth column, respectively. The superposition of two and three Gaussian functions are used to model the variables of μk and σk, respectively. The fitted results are listed in Tables 4 and 5.Table 4
Parameters controllingμ for each column.
Records1ε1ψ1s2ε2ψ2CHB002EW2.281.4511.881.5219.996.88TKY023EW1.52−13.0423.900.4719.262.95SIT003EW31.65−46.6327.610.8019.563.24Table 5
Parameters controlling varianceσ for each column.
Records1ε1ψ1s2ε2ψ2s3ε3ψ3CHB002EW1554.10145.9833.68−0.7321.164.381.2424.8721.47TKY023EW0.690.4710.110.3916.335.971.0089.0465.94SIT003EW1555.50154.6521.78−8.267.7841.059.075.2645.48As illustrated in Figure13, it is found that the parameter μ decreases with an increase in the sequence of columns in general. This is consistent with the fact that high-frequency components mainly appear at beginning, thus accounting for great contributions to the intensity of a ground motion; at the same time, low-frequency components gradually dominate over high-frequency components as time evolves. In contrast, the trend for the parameter σ appears to be complicated. It is possible that variations in bandwidth among time instants result in the variability of the time-frequency domain. This trend is similar to that depicted by using the lognormal probability distribution in paper [30].Figure 13
Variations of parametersμk and σk with respect to time instant. (a) CHB002EW and (b) TKY023EW.
(a)(b)It is shown from Figure14 that the stationary property of the Gaussian random process is transformed into nonstationary property of the initial seed motions. Furthermore, the properties of the initial seed motions are similar to those of their target motions, which are illustrated in Figure 15; note that this result is detailed later.Figure 14
Illustration of changes in the initial seed motions considering the time-frequency domain. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 15
Comparison of the simulations with the records with respect to wavelet coefficients. (a) CHB002EW and (b) TKY023EW.
(a)(b)


### 4.2. Simulation of High-Frequency Components Compatible with Target Motions

The methodology related to compatibility is that wavelet coefficients are scaled up or down among rows on the basis of the linear relationship between amplitudes of frequency components and an elastic response of a single degree of freedom system; meanwhile, coefficients are scaled among columns based on the squared relationship between amplitudes and cumulative energy [11].


#### 4.2.1. Compatibility with Frequency Response

Specifically, the scale factorcn+1j,ki is determined from the ratio of the targeted acceleration spectrum Satargetfi to the simulated spectrum Sasimulatednfi at frequency fi=idfw. Then, the new time series sn+1t is readily available from the inverse wavelet transform after the (n+1)th iteration. The iteration process is performed until the mean error is satisfied with the assigned threshold. Accordingly, the equations related to adjustment in frequency domain are expressed as follows:(20)cn+1ji=cnj,kiSatargetfiSasimulatednfi,i=1,2,…,2j,(21)sn+1t=∑i=12j∑k=12N−jcn+1j,kiψj,kit.It is essential to locate the instant timetk, at which the 5% damping single degree of the freedom oscillator with natural frequency fi presents the maximum response. This enables us to merely modify the amplitudes with centre frequency fi at the kth column instead of modifying the amplitudes with centre frequency fi over entire columns. In addition, the instant time tk should be within 0,tl. This is in accordance with the result of phase dispersion that the high-frequency components arrive earlier than long-period components.As suggested by the limitation of wavelet packets, there is low resolution at long periods. This means that the procedure cannot attain adequate accuracy during long periods ranging from 5 to 10 sec. Because of this, the frequency range in this process is defined within [0.3, 25 Hz], during which wavelet packets take advantage of high resolution to make simulations compatible with the targeted spectrum as accurately as possible.


#### 4.2.2. Compatibility with Cumulative Energy

The ratiown+1t for modification in time domain is calculated from the squared root of the incremental energy curve at the time interval tk−1,tk, and its expression is given as follows:(22)wn+1tk=Htargettk−Htargettk−1Hnsimulatedtk−Hnsimulatedtk−1,k=1,2,…,N2j.As mentioned in equation (13), a small frequency interval dfw=dt/2j+1 gives rise to a large time interval dtw=2jdt. Because of this, the ratios wn+1t are interpolated by using a cubic spline to obtain the ratios wn+1t′ at each time interval dt. Then, a new time series sn+2t is generated after multiplying the ratios wn+1t′. Its expression is written in equation (23). Note that the new series sn+2t should be baseline corrected through a high-pass filter with a low cutoff frequency of 0.1:(23)sn+2t=wn+1t′sn+1t.Consequently, it is found from Figure16 that the simulated acceleration spectra present close agreement with targeted motions, except for the longer period range (5–10 sec); for cumulative energy curves, the simulations nearly overlap with the records, suggesting that they share the same acceleration waveform. In addition, as illustrated in Figure 15, the simulations appear to be the similar distribution of wavelet coefficients to the records, since the majority of the large coefficients of the simulations are consistent with that of the records.Figure 16
Comparison of simulated high-frequency components with the records with respect to response spectrum and cumulative energy curve. (a) CHB002EW and (b) TKY023EW.
(a)(b)As suggested by Huang [11], two mean lognormal errors are employed to evaluate the accuracy of compatibility with the acceleration spectrum and cumulative energy curve, respectively. After 20 iterations in Figure 17, the errors for spectrum and cumulative energy curve are almost below 0.02 and 0.001, respectively.Figure 17
Mean squared error for each iteration.Accordingly, two samples of simulated high-frequency components are displayed in the form of the acceleration and velocity time series in Figure18. In addition, the velocity spectra attained from their acceleration time series are presented in Figure 19. The velocity spectra of simulated motions are similar to their original records at low to intermediate periods; however, the spectra of the original records at longer periods are significantly larger than those of the simulated motions. This means that the simulated motions are characterized by high-frequency components relative to their original records.Figure 18
Simulations of acceleration and velocity time series by means of wavelet packets. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 19
Comparison of original records with the simulations in terms of velocity spectrum. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 4.1. Generation of Initial Seed Motion

Prior to the generation of high-frequency components, an initial seed motion should be created, because variability of the initial seed motion could influence nonstationary properties of desired motions [11, 30]. Thus, the initial seed motions for high-frequency components are created as a result of the modulation of the Gaussian random process with respect to time and frequency domains.


### 4.1.1. Modulation in Time Domain

By visual observation of the acceleration waveforms shown in Figure2, a piecewise function is employed to modulate the Gaussian random process in time domain, and its mathematical expression is written as follows:(14)eht=0,t≤t0,t−t0tm−t0γ,t0<t≤tm,exp−ηt−tm,t>tm,where tm denotes the time corresponding to PGA; γ and η are parameters controlling the ascent and descent stage of the time segment before and after tm, respectively; t0 indicates the onset of a ground motion. Correspondingly, as shown in Figure 10, a set of normalized envelopes regarding various parameters of tm, γ, and η demonstrates the variability of waveform shapes.Figure 10
Piecewise function with various parameters oftm, γ, and η.Furthermore, a normalized cumulative energy model available from integration of the piecewise function with respect to time is given as follows:(15)ht=∫eh2tdt∫0Teh2t2dt=12η+tm−t0/2γ+10,t≤t0,tm−t02γ+1t−t0tm−t02γ+1,t0<t≤tm,−12ηexp−2ηt−tm+12η+tm−t02γ+1,t>tm.Thus, the parameters of the envelope function are attained by minimizing the integrated differences between the normalized actual energy curve and piecewise energy model, that is,(16)t^0,γ^,η^=argmint0,γ,η∫0t∫0th2tdt−∫0tha2tdt2dt,(17)hat=∫0tx2tdt∫0Dx2tdt,where xt is the acceleration time series, D is the duration of an acceleration time series, and hat denotes the normalized cumulative energy. The fitted parameters are listed in Table 3.Table 3
Parameters controlling envelope of motions.
Recordt0γηCHB002EW0.0114.450.061TKY023EW22.611.960.030SIT003EW0.0482.680.019It is found from Figure11 that the fitted curves are in close agreement with their actual curves at the strong shaking phases, during which the curves appear to be a vertical line. Meanwhile, discrepancies occur at approximately 140 sec, after which LALP surface waves arrive. Although the amplitudes of LALP surface waves are relatively small, in distinct contrast to those of body waves, they have a low attenuation rate. This results from the long-period properties of surface waves.Figure 11
Normalized cumulative energy curves for records and modulated Gaussian process. (a) CHB002EW and (b) TKY023EW.
(a)(b)


### 4.1.2. Modification in Frequency Domain

Wavelet coefficients in time-frequency domain comply with a bivariate lognormal distribution, which is used to present the variability of ground motions [10]. To simplify this, a set of shaped-like lognormal functions are used to describe the power spectrum density for each column (time instant). Thus, a set of squared wavelet coefficients derived from the modulated Gaussian process are modified by multiplying each fitted function to simulate the variability of targeted motions.The shaped-like lognormal function of thekth column with respect to frequency f is given as follows:(18)Lkfuk,σk=1fσk2πexp−lnf−uk22σk2,where the parameters of μk and σk for the kth column are attained by fitting the assigned function to squared coefficients at the kth column in the sense of power energy and Lk represents the power spectrum at the kth column.As a result, normalized squared coefficients for specific columns of CHB002EW and TKY023EW are shown in Figure12, together with the fitted functions. The fitted models generally reflect the distribution of the power spectra because their energy tends to concentrate at low frequencies as time columns evolve.Figure 12
Comparison of the fitted models with actual records in terms of normalized power spectrum at specific columns. (a) CHB002EW and (b) TKY023EW.
(a)(b)To establish a parametric relationship between the identified parameters and the sequence of columns (time instants), a generic expression built on the superposition of Gaussian functions is defined as follows:(19)λkk=∑isiexp−k−εiψi2,where subscript i indicates the ith Gaussian function; si is the scaling factor; εi indicates the location for peak value; ψi denotes the width of “bell” Gaussian function; and λk represents the variables of μk and σk in the kth column, respectively. The superposition of two and three Gaussian functions are used to model the variables of μk and σk, respectively. The fitted results are listed in Tables 4 and 5.Table 4
Parameters controllingμ for each column.
Records1ε1ψ1s2ε2ψ2CHB002EW2.281.4511.881.5219.996.88TKY023EW1.52−13.0423.900.4719.262.95SIT003EW31.65−46.6327.610.8019.563.24Table 5
Parameters controlling varianceσ for each column.
Records1ε1ψ1s2ε2ψ2s3ε3ψ3CHB002EW1554.10145.9833.68−0.7321.164.381.2424.8721.47TKY023EW0.690.4710.110.3916.335.971.0089.0465.94SIT003EW1555.50154.6521.78−8.267.7841.059.075.2645.48As illustrated in Figure13, it is found that the parameter μ decreases with an increase in the sequence of columns in general. This is consistent with the fact that high-frequency components mainly appear at beginning, thus accounting for great contributions to the intensity of a ground motion; at the same time, low-frequency components gradually dominate over high-frequency components as time evolves. In contrast, the trend for the parameter σ appears to be complicated. It is possible that variations in bandwidth among time instants result in the variability of the time-frequency domain. This trend is similar to that depicted by using the lognormal probability distribution in paper [30].Figure 13
Variations of parametersμk and σk with respect to time instant. (a) CHB002EW and (b) TKY023EW.
(a)(b)It is shown from Figure14 that the stationary property of the Gaussian random process is transformed into nonstationary property of the initial seed motions. Furthermore, the properties of the initial seed motions are similar to those of their target motions, which are illustrated in Figure 15; note that this result is detailed later.Figure 14
Illustration of changes in the initial seed motions considering the time-frequency domain. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 15
Comparison of the simulations with the records with respect to wavelet coefficients. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 4.1.1. Modulation in Time Domain

By visual observation of the acceleration waveforms shown in Figure2, a piecewise function is employed to modulate the Gaussian random process in time domain, and its mathematical expression is written as follows:(14)eht=0,t≤t0,t−t0tm−t0γ,t0<t≤tm,exp−ηt−tm,t>tm,where tm denotes the time corresponding to PGA; γ and η are parameters controlling the ascent and descent stage of the time segment before and after tm, respectively; t0 indicates the onset of a ground motion. Correspondingly, as shown in Figure 10, a set of normalized envelopes regarding various parameters of tm, γ, and η demonstrates the variability of waveform shapes.Figure 10
Piecewise function with various parameters oftm, γ, and η.Furthermore, a normalized cumulative energy model available from integration of the piecewise function with respect to time is given as follows:(15)ht=∫eh2tdt∫0Teh2t2dt=12η+tm−t0/2γ+10,t≤t0,tm−t02γ+1t−t0tm−t02γ+1,t0<t≤tm,−12ηexp−2ηt−tm+12η+tm−t02γ+1,t>tm.Thus, the parameters of the envelope function are attained by minimizing the integrated differences between the normalized actual energy curve and piecewise energy model, that is,(16)t^0,γ^,η^=argmint0,γ,η∫0t∫0th2tdt−∫0tha2tdt2dt,(17)hat=∫0tx2tdt∫0Dx2tdt,where xt is the acceleration time series, D is the duration of an acceleration time series, and hat denotes the normalized cumulative energy. The fitted parameters are listed in Table 3.Table 3
Parameters controlling envelope of motions.
Recordt0γηCHB002EW0.0114.450.061TKY023EW22.611.960.030SIT003EW0.0482.680.019It is found from Figure11 that the fitted curves are in close agreement with their actual curves at the strong shaking phases, during which the curves appear to be a vertical line. Meanwhile, discrepancies occur at approximately 140 sec, after which LALP surface waves arrive. Although the amplitudes of LALP surface waves are relatively small, in distinct contrast to those of body waves, they have a low attenuation rate. This results from the long-period properties of surface waves.Figure 11
Normalized cumulative energy curves for records and modulated Gaussian process. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 4.1.2. Modification in Frequency Domain

Wavelet coefficients in time-frequency domain comply with a bivariate lognormal distribution, which is used to present the variability of ground motions [10]. To simplify this, a set of shaped-like lognormal functions are used to describe the power spectrum density for each column (time instant). Thus, a set of squared wavelet coefficients derived from the modulated Gaussian process are modified by multiplying each fitted function to simulate the variability of targeted motions.The shaped-like lognormal function of thekth column with respect to frequency f is given as follows:(18)Lkfuk,σk=1fσk2πexp−lnf−uk22σk2,where the parameters of μk and σk for the kth column are attained by fitting the assigned function to squared coefficients at the kth column in the sense of power energy and Lk represents the power spectrum at the kth column.As a result, normalized squared coefficients for specific columns of CHB002EW and TKY023EW are shown in Figure12, together with the fitted functions. The fitted models generally reflect the distribution of the power spectra because their energy tends to concentrate at low frequencies as time columns evolve.Figure 12
Comparison of the fitted models with actual records in terms of normalized power spectrum at specific columns. (a) CHB002EW and (b) TKY023EW.
(a)(b)To establish a parametric relationship between the identified parameters and the sequence of columns (time instants), a generic expression built on the superposition of Gaussian functions is defined as follows:(19)λkk=∑isiexp−k−εiψi2,where subscript i indicates the ith Gaussian function; si is the scaling factor; εi indicates the location for peak value; ψi denotes the width of “bell” Gaussian function; and λk represents the variables of μk and σk in the kth column, respectively. The superposition of two and three Gaussian functions are used to model the variables of μk and σk, respectively. The fitted results are listed in Tables 4 and 5.Table 4
Parameters controllingμ for each column.
Records1ε1ψ1s2ε2ψ2CHB002EW2.281.4511.881.5219.996.88TKY023EW1.52−13.0423.900.4719.262.95SIT003EW31.65−46.6327.610.8019.563.24Table 5
Parameters controlling varianceσ for each column.
Records1ε1ψ1s2ε2ψ2s3ε3ψ3CHB002EW1554.10145.9833.68−0.7321.164.381.2424.8721.47TKY023EW0.690.4710.110.3916.335.971.0089.0465.94SIT003EW1555.50154.6521.78−8.267.7841.059.075.2645.48As illustrated in Figure13, it is found that the parameter μ decreases with an increase in the sequence of columns in general. This is consistent with the fact that high-frequency components mainly appear at beginning, thus accounting for great contributions to the intensity of a ground motion; at the same time, low-frequency components gradually dominate over high-frequency components as time evolves. In contrast, the trend for the parameter σ appears to be complicated. It is possible that variations in bandwidth among time instants result in the variability of the time-frequency domain. This trend is similar to that depicted by using the lognormal probability distribution in paper [30].Figure 13
Variations of parametersμk and σk with respect to time instant. (a) CHB002EW and (b) TKY023EW.
(a)(b)It is shown from Figure14 that the stationary property of the Gaussian random process is transformed into nonstationary property of the initial seed motions. Furthermore, the properties of the initial seed motions are similar to those of their target motions, which are illustrated in Figure 15; note that this result is detailed later.Figure 14
Illustration of changes in the initial seed motions considering the time-frequency domain. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 15
Comparison of the simulations with the records with respect to wavelet coefficients. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 4.2. Simulation of High-Frequency Components Compatible with Target Motions

The methodology related to compatibility is that wavelet coefficients are scaled up or down among rows on the basis of the linear relationship between amplitudes of frequency components and an elastic response of a single degree of freedom system; meanwhile, coefficients are scaled among columns based on the squared relationship between amplitudes and cumulative energy [11].


### 4.2.1. Compatibility with Frequency Response

Specifically, the scale factorcn+1j,ki is determined from the ratio of the targeted acceleration spectrum Satargetfi to the simulated spectrum Sasimulatednfi at frequency fi=idfw. Then, the new time series sn+1t is readily available from the inverse wavelet transform after the (n+1)th iteration. The iteration process is performed until the mean error is satisfied with the assigned threshold. Accordingly, the equations related to adjustment in frequency domain are expressed as follows:(20)cn+1ji=cnj,kiSatargetfiSasimulatednfi,i=1,2,…,2j,(21)sn+1t=∑i=12j∑k=12N−jcn+1j,kiψj,kit.It is essential to locate the instant timetk, at which the 5% damping single degree of the freedom oscillator with natural frequency fi presents the maximum response. This enables us to merely modify the amplitudes with centre frequency fi at the kth column instead of modifying the amplitudes with centre frequency fi over entire columns. In addition, the instant time tk should be within 0,tl. This is in accordance with the result of phase dispersion that the high-frequency components arrive earlier than long-period components.As suggested by the limitation of wavelet packets, there is low resolution at long periods. This means that the procedure cannot attain adequate accuracy during long periods ranging from 5 to 10 sec. Because of this, the frequency range in this process is defined within [0.3, 25 Hz], during which wavelet packets take advantage of high resolution to make simulations compatible with the targeted spectrum as accurately as possible.


### 4.2.2. Compatibility with Cumulative Energy

The ratiown+1t for modification in time domain is calculated from the squared root of the incremental energy curve at the time interval tk−1,tk, and its expression is given as follows:(22)wn+1tk=Htargettk−Htargettk−1Hnsimulatedtk−Hnsimulatedtk−1,k=1,2,…,N2j.As mentioned in equation (13), a small frequency interval dfw=dt/2j+1 gives rise to a large time interval dtw=2jdt. Because of this, the ratios wn+1t are interpolated by using a cubic spline to obtain the ratios wn+1t′ at each time interval dt. Then, a new time series sn+2t is generated after multiplying the ratios wn+1t′. Its expression is written in equation (23). Note that the new series sn+2t should be baseline corrected through a high-pass filter with a low cutoff frequency of 0.1:(23)sn+2t=wn+1t′sn+1t.Consequently, it is found from Figure16 that the simulated acceleration spectra present close agreement with targeted motions, except for the longer period range (5–10 sec); for cumulative energy curves, the simulations nearly overlap with the records, suggesting that they share the same acceleration waveform. In addition, as illustrated in Figure 15, the simulations appear to be the similar distribution of wavelet coefficients to the records, since the majority of the large coefficients of the simulations are consistent with that of the records.Figure 16
Comparison of simulated high-frequency components with the records with respect to response spectrum and cumulative energy curve. (a) CHB002EW and (b) TKY023EW.
(a)(b)As suggested by Huang [11], two mean lognormal errors are employed to evaluate the accuracy of compatibility with the acceleration spectrum and cumulative energy curve, respectively. After 20 iterations in Figure 17, the errors for spectrum and cumulative energy curve are almost below 0.02 and 0.001, respectively.Figure 17
Mean squared error for each iteration.Accordingly, two samples of simulated high-frequency components are displayed in the form of the acceleration and velocity time series in Figure18. In addition, the velocity spectra attained from their acceleration time series are presented in Figure 19. The velocity spectra of simulated motions are similar to their original records at low to intermediate periods; however, the spectra of the original records at longer periods are significantly larger than those of the simulated motions. This means that the simulated motions are characterized by high-frequency components relative to their original records.Figure 18
Simulations of acceleration and velocity time series by means of wavelet packets. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 19
Comparison of original records with the simulations in terms of velocity spectrum. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 4.2.1. Compatibility with Frequency Response

Specifically, the scale factorcn+1j,ki is determined from the ratio of the targeted acceleration spectrum Satargetfi to the simulated spectrum Sasimulatednfi at frequency fi=idfw. Then, the new time series sn+1t is readily available from the inverse wavelet transform after the (n+1)th iteration. The iteration process is performed until the mean error is satisfied with the assigned threshold. Accordingly, the equations related to adjustment in frequency domain are expressed as follows:(20)cn+1ji=cnj,kiSatargetfiSasimulatednfi,i=1,2,…,2j,(21)sn+1t=∑i=12j∑k=12N−jcn+1j,kiψj,kit.It is essential to locate the instant timetk, at which the 5% damping single degree of the freedom oscillator with natural frequency fi presents the maximum response. This enables us to merely modify the amplitudes with centre frequency fi at the kth column instead of modifying the amplitudes with centre frequency fi over entire columns. In addition, the instant time tk should be within 0,tl. This is in accordance with the result of phase dispersion that the high-frequency components arrive earlier than long-period components.As suggested by the limitation of wavelet packets, there is low resolution at long periods. This means that the procedure cannot attain adequate accuracy during long periods ranging from 5 to 10 sec. Because of this, the frequency range in this process is defined within [0.3, 25 Hz], during which wavelet packets take advantage of high resolution to make simulations compatible with the targeted spectrum as accurately as possible.


## 4.2.2. Compatibility with Cumulative Energy

The ratiown+1t for modification in time domain is calculated from the squared root of the incremental energy curve at the time interval tk−1,tk, and its expression is given as follows:(22)wn+1tk=Htargettk−Htargettk−1Hnsimulatedtk−Hnsimulatedtk−1,k=1,2,…,N2j.As mentioned in equation (13), a small frequency interval dfw=dt/2j+1 gives rise to a large time interval dtw=2jdt. Because of this, the ratios wn+1t are interpolated by using a cubic spline to obtain the ratios wn+1t′ at each time interval dt. Then, a new time series sn+2t is generated after multiplying the ratios wn+1t′. Its expression is written in equation (23). Note that the new series sn+2t should be baseline corrected through a high-pass filter with a low cutoff frequency of 0.1:(23)sn+2t=wn+1t′sn+1t.Consequently, it is found from Figure16 that the simulated acceleration spectra present close agreement with targeted motions, except for the longer period range (5–10 sec); for cumulative energy curves, the simulations nearly overlap with the records, suggesting that they share the same acceleration waveform. In addition, as illustrated in Figure 15, the simulations appear to be the similar distribution of wavelet coefficients to the records, since the majority of the large coefficients of the simulations are consistent with that of the records.Figure 16
Comparison of simulated high-frequency components with the records with respect to response spectrum and cumulative energy curve. (a) CHB002EW and (b) TKY023EW.
(a)(b)As suggested by Huang [11], two mean lognormal errors are employed to evaluate the accuracy of compatibility with the acceleration spectrum and cumulative energy curve, respectively. After 20 iterations in Figure 17, the errors for spectrum and cumulative energy curve are almost below 0.02 and 0.001, respectively.Figure 17
Mean squared error for each iteration.Accordingly, two samples of simulated high-frequency components are displayed in the form of the acceleration and velocity time series in Figure18. In addition, the velocity spectra attained from their acceleration time series are presented in Figure 19. The velocity spectra of simulated motions are similar to their original records at low to intermediate periods; however, the spectra of the original records at longer periods are significantly larger than those of the simulated motions. This means that the simulated motions are characterized by high-frequency components relative to their original records.Figure 18
Simulations of acceleration and velocity time series by means of wavelet packets. (a) CHB002EW and (b) TKY023EW.
(a)(b)Figure 19
Comparison of original records with the simulations in terms of velocity spectrum. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 5. Combination and Verification

After the long-period and high-frequency components are simulated individually, they are incorporated into each other in time domain. Accordingly, Figure20 shows four samples of the simulations for the records CHB002EW and TKY023EW. From visual inspection, four realizations are similar to their actual records in terms of acceleration waveforms, since their strong shaking process occurs at approximately 100 sec, which is consistent with the actual records seen in Figure 2. In addition, the obvious presence of LALP surface waves in simulated velocity waveforms highlights their long-period properties. However, it is noted that the components with an estimated period of 20 sec are absent from the simulated velocity series. This results from the low resolution of wavelet packets in discerning components with periods larger than 10 sec. However, it seems that such long-period components are not of interest in engineering practice because the resonant response of structures is mainly attributed to the components with periods ranging from 5 to 10 sec [31].Figure 20
Samples of the simulations in the form of acceleration and velocity time series. (a) CHB002EW and (b) TKY023EW.
(a)(b)Then, comparisons of four realizations with actual records in terms of the cumulative energy curve and Fourier amplitude spectrum are illustrated in Figure21. It is clear that each sample matches well with the actual records considering the cumulative energy curve. In addition, the Fourier amplitude spectra of the realizations are close to their actual records over frequencies ranging from 0.1 Hz to 20 Hz. However, distinct discrepancies appear at frequencies below 0.1 Hz, resulting from the limitations of wavelet packets.Figure 21
Comparison between four realizations and the recorded motions in terms of Fourier amplitudes and cumulative energy. (a) CHB002EW and (b) TKY023EW.
(a)(b)Furthermore, variations of the simulated long-period motions are mainly dependent on the simulated long-period components generated by means of spectral representation, when the simulated long-period motions are verified in sense of compatibility with a target response spectrum.In this regard, 100 realizations of long-period components are individually incorporated with one simulated high-frequency motion which consists of high-frequency components in order to form a set of simulated long-period motions, for which 100 response spectra with 5% damping are attained. As shown in Figure22, the median spectra of the simulated motions are closely aligned with the target spectra. Two dominant periods for the high-frequency and long-period components, each of which are 0.2 and 5.3 sec for CHB002EW and 0.7 and 6.6 sec for TKY023EW, are captured well by the median simulations. This suggests that the simulated long-period motions replicate the properties of the actual motions. Meanwhile, it is observed that the actual spectra at loner periods from 4 to 10 sec mostly fall within the bound of the standard deviation ±σ, and thus, this means that the actual motions could be one of the ensembles of the simulated motions in the stochastic sense.Figure 22
Response spectra (solid blue line) for the records CHB002EW and TKY023EW. Of 100 simulated long-period ground motions (thin gray lines), their median (thick red line) and median plus and minus one standard deviation±σ (dashed lines) are involved; acceleration spectra on the left side; velocity spectra on the right side. (a) CHB002EW and (b) TKY023EW.
(a)(b)


## 6. Summary

In this paper, a method for generating the long-period ground motions based on the combination of spectral representation and wavelet packets is proposed. A flowchart regarding this method is illustrated in Figure23. It is clear that the proposed method is divided into three parts: (I) spectral representation for long-period components, (II) simulation of high-frequency components using wavelet packets, and (III) superposition of the former into the latter in time domain. To depict clearly these procedures, an illustrative example for SIT003EW, which is another record in Tohuku earthquake M9, is presented in Figure 24 as follows.Figure 23
Schematic plot for generating long-period ground motions.Figure 24
Example to illustrate the proposed method.For the first phase, the LALP velocity surface waves are truncated from the real motion with reference totl shown in Figure 24(a), conditioned by the truncated segment showing an agreement with its original record displayed in Figure 24(b); the power spectrum density for the truncated velocity segment is estimated shown in Figure 24(c); instantaneous amplitudes of the velocity time series are obtained by multimodal point method [29]; because of these, the simulated LALP velocity surface waves are generated by means of spectral representation shown in Figure 24(d).For the second phase, an initial seed motion is generated from the modification process, in which a Gaussian random process is modified in time domain by using a piecewise function and in frequency domain by using a set of fitted power spectrum functions; accordingly, the adjusted seed motion is presented in Figure24(e); then, the iterative adjustments result in the high-frequency motion matching well with the target acceleration spectrum at short periods and cumulative energy curve, and the corresponding product is shown in Figure 24(f).For the final phase, the LALP acceleration time series is attained from the differentiation of the resulting velocity series, and then it is superimposed into the high-frequency motion; consequently, one sample for the simulated long-period motions by the proposed method is shown in the form of acceleration and velocity in Figures24(g) and 24(h), respectively. Furthermore, 100 realizations are generated to verify the matching well to the target spectra, as illustrated in Figures 24(i) and 24(j).It should be pointed out that the proposed method towards representative records is presented herein to demonstrate the feasibility for generating long-period motions. In addition, parameter values that describe stochastic property, LALP surface waves, and high-frequency components are determined individually for each record. Because of this, the methodology can be applied to other set of the motions with similar properties.Finally, it is noted that the proposed method is suitable for simulation of long-period ground motions with narrowband long-period properties. This results from the analytical model for the power spectrum density. In addition, provided that the later-arriving long-period components are not capable of representing the gross long-period properties, the criterion for determination of long-period components may appear to be less effective.On the other hand, future studies appear to investigate a relationship between the identified parameters and seismic variables and extend the method to simulate long-period motion based on specific seismic environment, thus improving the engineering applicability.



---

*Source: 2381080-2020-05-31.xml*