Authors: H. P. Hong, X. Z. Cui, D. Qiao
Categories: Research, Continuous wavelet transforms, Nonstationary and non-Gaussian process, S-transform, Seismic ground motions, Simulation, Wind velocity
Source: Journal of Infrastructure Preservation and Resilience
We proposed a new iterative power and amplitude correction (IPAC) algorithm to simulate nonstationary and non-Gaussian processes. The proposed algorithm is rooted in the concept of defining the stochastic processes in the transform domain, which is elaborated and extend. The algorithm extends the iterative amplitude adjusted Fourier transform algorithm for generating surrogate and the spectral correction algorithm for simulating stationary non-Gaussian process. The IPAC algorithm can be used with different popular transforms, such as the Fourier transform, S-transform, and continuous wavelet transforms. The targets for the simulation are the marginal probability distribution function of the process and the power spectral density function of the process that is defined based on the variables in the transform domain for the adopted transform. The algorithm is versatile and efficient. Its application is illustrated using several numerical examples.
Keywords: Simulation, Nonstationary and non-Gaussian process, S-transform, Continuous wavelet transforms, Seismic ground motions, Wind velocity
Observed time histories of the seismic ground motions [31], wind velocity [48], wave height [32], etc. fluctuate randomly in time and space. The time histories are used as the input to carry out the structural analysis. Since the available recorded time histories of the random phenomena are limited, their synthetics are generated and used in practice. The simulation is based on the theory of stochastic processes [7, 22, 25]. The simulated ground motions are used to assess seismic risk of infrastructures such as dams and portfolio of buildings (e.g., [3, 24]). The simulated winds are employed to assess the performance structures and infrastructure system [55, 60].
For stationary Gaussian processes, and evolutionary processes [38], the simulation can be carried out using the spectral representation method [23, 47], developed based on the ordinary Fourier transform (FT). A stationary process is defined by its power spectral density (PSD) function, and an evolutionary process is defined by the evolutionary PSD that is a function of an amplitude modulation function. The evolutionary process with time-dependent amplitude modulation is widely used in generating seismic ground motions [1, 10] and fluctuating wind velocity for high-intensity wind events [5, 17, 20, 21].
Masters and Gurley [27] proposed an iterative spectral correction algorithm to simulate the stationary non-Gaussian processes, where the spectral representation method is used in each iteration to generate the time history. They showed that their algorithm outperforms the SRM-based simulation techniques presented in Yamazaki and Shinozuka [58], Gurley and Kareem [15], Grigoriu [14], and Deodatis and Micaletti [11]. It is noted that an algorithm similar to the spectral correction algorithm, namely the iterative amplitude adjusted Fourier transform (IAAFT) algorithm, was proposed by Schreiber and Schmitz [43, 44] in the context of generating surrogate for statistical hypothesis testing. The use of the translation process for the stationary non-Gaussian process proposed in Grigoriu [14] was extended for the nonstationary processes by others, including Ferrante et al. [12] and Shields et al. [46].
The evolutionary PSD is often assessed using (time-dependent) windowed Fourier transform, such as the short-time Fourier transform [6]. The resolution of such a transform is controlled by the width of the window. As the width of the window increases, a better resolution is obtained at the low frequencies, and the resolution in time deteriorates. A good resolution in both time and frequency (i.e., scale) can be obtained by applying the continuous wavelet transforms (WT) [9, 35]. A procedure to estimate the evolutionary PSD by applying the continuous WT was proposed by Spanos and Failla [50]. However, an algorithm that directly applies the continuous WT to simulate the nonstationary stochastic processes with a prescribed wavelet spectrum or time-scale PSD was unavailable. Recently, an iterative algorithm was presented by Chavez and Cazelles [4] to generate surrogate for statistical hypothesis testing. We will point out, in the following sections, a potential weakness of the algorithm, as well as the link between this algorithm and an interesting way of defining nonstationary processes in the wavelet domain introduced by Maraun et al. [26]. The lack of continuous WT-based algorithm to simulate time histories is partly due to that the use of continuous WT does not lead to the decomposed signal to be represented by a set of orthogonal basis functions. Rather than using the continuous WT, the application of the discrete WT and wavelet packet transform that have orthogonal basis functions is presented in Gurley and Kareem [16], and Yamamoto and Baker [57]. The resolution obtained by using these discrete transforms is less refined than that obtained by using the continuous WT.
The phase information in WT is local, while the phase information in the Fourier transform refers to the harmonics at zero time [52]. Stockwell et al. [53] (see also [37]) developed the S-transform (ST) that provides the time-frequency representation of the analyzed signal. It is a hybrid of continuous WT and windowed FT. The S-transform provides frequency-dependent resolution. Similar to the continuous WT, ST does not lead to a decomposed signal to be represented by a set of orthogonal basis functions. Stockwell [52] proposed a discrete orthonormal S-transform. The simulation of the seismic ground motions by using the discrete orthonormal S-transform or combination with ST was presented in Cui and Hong [8] and Hong and Cui [19]. However, an algorithm by using ST alone to simulate nonstationary stochastic processes has not been reported in the literature.
There are other techniques used to simulate the nonstationary processes. These include the application of autoregressive moving-average (ARMA) [42], Karhunen–Loéve expansion [36, 49, 51], and polynomial chaos [41], and Hilbert-Huang transform [56]. The ARMA uses the recursive relations that connect the random quantity to be simulated at successive time increments. The application of the Karhunen–Loéve expansion uses the eigenfunctions that are obtained by solving Fredholm integral equation of the second kind. In the Hilbert-Huang transform, a random process is decomposed into intrinsic mode functions by the method of empirical mode decomposition. A review of these simulation procedures is beyond the present study since these techniques involve varieties of mathematical concepts and algorithms.
In the present study, we exam and extend the definition of the nonstationary processes in the transform domain. We proposed an iterative power and amplitude correction (IPAC) algorithm to simulate nonstationary and non-Gaussian processes. The algorithm could be viewed as an extension of IAAFT [43] and the spectral correction algorithm [27] and is rooted in the concept of defining the stochastic processes in the transform domain. In particular, we provide details of using the proposed algorithm with FT, ST and WT, where the energy distribution in the transform domain that satisfies energy preservation is prescribed, and the marginal probability distribution function of the process is given. We provide numerical examples to show the proposed algorithm and compare the simulated time histories obtained by using the ST-based and (continuous) WT-based approach.
This section summarizes some basic properties of FT [6, 30], ST [37, 53], and continuous WT [9]. Only the continuous WT, including its discretized form (which differs from the discrete wavelet transform), is used in the present study. Unless otherwise indicated, WT refers to the continuous WT and its discretized form in the following. The summary provides the basis for the proposed iterative simulation algorithm to be described in the next sections.
Let x(t) denote a realization of a stochastic process such as the ground motion record, X(t). FT of x(t), and its inverse (IFT) can be expressed as
and,
where FT(•) and IFT(•) denote the FT and IFT operations, the subscript associated with these operators indicates the domain or the index where the operation is carried out; x^f denotes FT of x(t); f is the frequency in Hz, x^f=x^∗−f, and * denotes the complex conjugate. A symbol or function with a circumflex is used to represent its FT throughout the present study. If x(t) is given in the discrete form x(jΔ~t), j = 0, …, N − 1, with a sampling time interval Δt~ and the duration T, T = NΔ~t~, the (discretized) FT pair is given by,
and,
where Δ~f~ = 1/T, and the operators FT(•) and IFT(•) that are used for continuous FT are used for discrete FT as well. It is considered implicitly in the following that the numerical calculations of x^pΔf and x(jΔ~t) are to be carried by using the fast Fourier transform (FFT) [30] for computational efficiency. Moreover, the notation {•}N~ is used for the collection of its argument of length N. For example, {x(jΔ~t)}N~ represents all x(jΔ~t~) for j = 0, …, N − 1.
ST of x(t) is defined as [37, 53],
where xSfτ is the ST coefficient, ST(·) denotes the S-transform of its argument, and τ is the center of the window function w(f, τ − t) defined as,
The parameter κ in Eq. (6) controls the effective width of the window in ST. It can be shown [53] that,
and,
where IST(·) is the inverse S-transform (IST). Using Eqs. (7) and (8), the discretized version of x(t) and xSfτ, represented by x(jΔ~t~) and xSqΔfpΔt pair, can be written as,
and,
indicating that the evaluation of the ST coefficients at (pΔ~f, qΔt) and its inverse at jΔt~ is based on FT.
WT is defined as [9, 35],
where xWsτ is the wavelet coefficient, the operator WT(•) denotes WT, ψ(⋅) is the mother wavelet and, s is the scaling or dilation factor, and τ is the translation or position parameter. Eq. (11) can be expressed as [9, 35],
to facilitate its computation by using FFT for signals given in the discretized form. If the admissibility condition 0 < Cψ < ∞ is satisfied, where Cψ=∫−∞∞1/f×ψ^f2df, x(t) can be obtained using the following inverse WT [9],
where IWT(•) is the inverse of WT(•). If ψ(t) = ψ ∗ (−t), Eq. (13) becomes,
Moreover, if the analytical wavelet – complex-valued wavelet function that its FT is null for negative frequency – is used, Eq. (13) can be expressed in Morlet formulation [45],
where C1ψ=∫0∞ψ^∗f/fdf and Re() denotes the real part of a complex number.
There are several well-known wavelet families [9, 33, 35], including Daubechies wavelets, generalized Morse wavelets, and Morlet wavelets.
Eqs. (14), (15) and (12) can be written in the following discretized form,
and,
where c0 and s0 are parameters for the numerical computation; K is the total number of scales considered for the numerical integration; TL = − Ll and TΔtU = LU define the lower and upper limits for the integral over time τ, and C1ψ=∫0∞ψ^∗f/fdf. In the following, we restrict ourselves to the real-valued signal and the analytical wavelets or wavelets with ψ(t) = ψ^∗^(−t).Δt
According to the spectral representation method [47] with the use of FFT [59], a sample of a Gaussian stationary process, x(t), can be simulated by transforming Gaussian white noise w(t) to the Fourier domain, multiplying it with an intensity function y^f, and transforming it back to the time domain. That is,
where eiθFwt=ηFTwt, in which the function η(C) = C/|C| is introduced to normalize the complex number C. Based on FT pair, x^f=yfeiθFwt. Since, by definition, the double-sided PSD function of the process x(t) with duration T, SFxf, is given by,
it indicates that given the target PSD function SFxf, one could define a stationary Gaussian process in the Fourier domain by assigning y^f=x^f=SFxfT. The samples of the process so defined can be obtained using,
and the expected PSD of the sampled signals equals the prescribed SFxf. The use of the definition given in Eq. (20) preserves the energy of x(t) according to Parseval’s theory.
We note that by assigning y^f equal to Mtx^f, Eq. (19) becomes,
which simulates a uniformly amplitude modulated evolutionary process [38]. Such a process has an evolutionary PSD function equals Mt2x^fx^∗f/T, and M(t) is the amplitude modulation function, which will be considered to be positive. However, the use of |y(f)| equal to Mtfx^f in Eq. (19) does not lend itself to be interpreted as a proper inverse Fourier transform because the modulation function depends on the frequency. This reduces the computational efficiency that otherwise can be gained by using FFT; it also makes the distinction between the modulation function and intensity function more blurred. We will concentrate only on the case where the modulation function is defined outside of the transform domain. However, the consideration of modulation that depends on variables in the transform domain could be a valid assumption.
Maraun et al. [26] emphasized the usefulness of using Eq. (19) to obtain samples of stationary Gaussian process, and extended it to define a class of nonstationary Gaussian processes in the wavelet domain by the wavelet multipliers yWsτ, indicating that an individual process is defined by its multipliers and a synthesizing wavelet pair. Samples of x(t) based on such a definition are then given as,
where eiθWwt=ηWTwt. We use the intensity function yWsτ and eiθWwt in Eq. (23) instead of using yWsτ and WT(w(t)) as suggested in Maraun et al. [26]. The use of eiθWwt instead of WT(w(t)) is aimed at not biasing the energy arising from the intensity function since [WT(w(t))][WT(w(t))]^∗^ is not a constant in the wavelet domain by using WT defined in Eq. (12). The use of yWsτ (as well as y^f in Eqs. (21) and (22)) is more restrictive than yWsτ but is adequate for the proposed algorithm in the following section since we are focused on real-valued signals. However, a negatively valued intensity and complex-valued intensity may be considered for other applications.
Similar to the use of M(t) in defining the uniformly modulated evolutionary process mentioned earlier, we include M(t) in Eq. (23),
to define a modulated and intensity function adjusted (MODIF) process. The intensity function gives time-scale characteristics of the process, and the modulation function provides additional time-varying characteristics of the process.
We further extend the concept of defining the MODIF process in the time-frequency domain according to ST, denoted as the S-domain, where samples of x(t) are given as,
where ySfτ is an intensity function in the S-domain, and eiθSwt=ηSTwt.
It is noted that besides the above-mentioned transforms, there are other transforms used for signal analysis and modeling; for example, the generalized Fourier family transforms [2]. Therefore, it is relevant and straightforward to conceptually generalize the approach in defining the MODIF processes in the transform domain if other transform pair is considered. The definitions lend themselves to an easily understandable and almost trivial algorithm to simulate stochastic
Step C) is separated from Steps A) and B) and is not affected by the selected transformation. A critical issue of applying the MODIF process with prescribed target energy distribution is that the energy distribution of the sampled signals for given intensity function may not be readily established, except for the case where FT is used (i.e., transforms with non-redundant representation). This is because unlike the FT, both WT and ST provide redundant representation. The redundant representation results in that, in general, yWsτeiθWwt and ySsτeiθSwt do not represent the proper coefficients of WT and ST, respectively. In other words, yWsτeiθWwt and ySsτeiθSwt are not equal to xWsτ=WTIWTyWsτeiθWwt and xSfτ=STISTySsτeiθSwt, respectively.
To see the impact of this inequality on the simulated MODIF process by using Eq. (24), we note that we can define the double-sided time-scale PSD (TSPSD) function of the simulated process x(t), SWxsτ, as,
The use of this definition leads to energy preservation since the integral of SWxsτ in the wavelet domain equals the integral of |x(t)|^2^ in the time domain (see Proposition 2.4.1 in Daubechies [9]). Consequently, even we assign yWsτ equals SWxsτCψs and M(t) = 1 for the simulation, the average energy of the sampled signals according to Eq. (24) will likely deviate from the specified target SWxsτ.
Consider that we simulate the MODIF process using Eq. (25). We can define the double-sided time-frequency PSD (TFPSD) function of the simulated process, SSxfτ, as,
since the use of this definition leads to energy preservation [18], where Dκ=∫−∞∞1/ζ×exp−2πκζ−12dζ. However, the average energy of the sampled signals by using Eq. (25) with ySfτ equal to xSfτ=SSxfτDκf and M(t) = 1 will likely deviate from the specified target.
In addition to the discussed energy distortion, the application of the MODIF process is likely to lead to the samples obtained from Eqs. (22), (24) and (25) to follow a marginal cumulative distribution function (CDF) that deviates from the prescribed marginal CDF of the zero-mean process FX, t(x(t)). An iterative process is proposed in the following sections to simulate the nonstationary and non-Gaussian with prescribed target PSD and CDF. The PSD functions that satisfy the energy preservation by considering the selected transform are used as the basis to describe the proposed algorithm to maintain consistency. Although this could become clumsy in some instances, it is useful in checking that a consistent transform pair is employed.
To develop the proposed iterative algorithm, we note that, given the observed {x(jΔ~t)}N, the IAAFT algorithm was proposed by Schreiber and Schmitz [43, 44] in the context of generating surrogates for statistical hypothesis testing. The algorithm repeatedly uses FT and IFT, and ranked data. This algorithm is explained using the ranking of x(jΔt~) in the following.
The PSD function SFxf of {x(jΔ~t)}N~ is calculated using Eq. (19) with possible smoothing. The objective of IAFFT is to generate surrogates that match the calculated SFxf and shuffled {x(jΔ~t)}N. A similar algorithm - the spectral correction algorithm - was independently designed by Masters and Gurley [27] to simulate non-Gaussian processes for the given target SFxf and target marginal CDF FX(x(t)). A subtle difference between these two algorithms is how the prescribed target PSD function and CDF are obtained or assigned. For example, {x(jΔt)}N~ is obtained through distribution mapping in the spectral correction algorithm. In IAAFT, {x(jΔ~t)}N~ is given and shuffled. This shuffling, in the spectral correction method, can be viewed as matching the prescribed probability distribution. Once {x(jΔ~t)}N~ is prescribed and SFxf is calculated, by letting {ξ(j)}N equal to the ascendingly sorted {x(jΔ~t)}N~, the steps of the IAAFT algorithm
Steps 1) and 2) are the same as Steps A) and B) described earlier that simulates a Gaussian process, except an additional ranking of xPC(jΔ~t) is carried out, which is equivalent to define the CDF as a preparation for the iteration. In general, Step 2) leads to xPC(jΔt) with the PSD correction but may deviate from the target CDF assigned by {ξ(j)}N, and Step 3) leads to the sampled xAC(jΔt) with the amplitude correction (i.e., matching CDF assigned based on {ξ(j)}N) but may deviate from the target PSD. The iteration adjusts the PSD and CDF of the sampled time series to their corresponding targets. The tolerable differences between xPC(jΔt) and xAC(jΔt) can be used as the convergence criterion. Once convergence is achieved xPC(jΔt) or xAC(jΔt) can be used as the sampled time series. Note that since the time increment equals Δt, the corresponding Nyquist frequency for the considered signal equals 1/(2Δt~).
The IAAFT algorithm is designed for stationary processes. For the shuffling of {x(jΔ~t)}N~ to simulation stationary process, it is implicitly considered that the marginal CDF of x(t) at any given time remains to be the same. Also, the PSD function for the stationary process is time-independent. The IAAFT algorithm or the spectral correction method is not applicable to simulate nonstationary processes as they have time-varying PSD and CDF.
In this section, we describe the proposed iterative power and amplitude correction (IPAC) algorithm to simulate the time history {x(jΔ~t)}N~ of a zero-mean nonstationary non-Gaussian process. Since Δ~t~ is assigned, N can be determined based on the length of the signal to be simulated and vice versa, and the Nyquist frequency for the sampled signal equals 1/(2Δ~t). The proposed algorithm could be viewed as an extension to the IAAFT algorithm. For the simulation, it is considered that, for M(t) = 1, the PSD function of the process that is characterized based on FT, or ST, or WT is given, and the distribution type for the marginal CDF of x(t), FX, t(x(t)), is known. Moreover, it is considered that FX, t~(x(t)) can be completely defined by the zero-mean, the time-varying standard deviation, σ(t), and other prescribed distribution parameters if they are required (since, in some cases, a CDF with more than two parameters may be considered).
If FT is considered for a stationary process, the standard deviation σ(t) equals λFx which is time-independent, where λFx, equals the integral of SFxf over the frequency domain. Since SSxfτ provides the energy distribution over the time-frequency domain, the integral of SSxfτ over the frequency domain provides the energy distribution in the time domain, λSxτ,
and the integral of SSxfτ over the time domain provides the energy distribution in the frequency domain, SSxf. Analogously to the statistics for the stationary process, λSxτ represents the variance of x(τ), and σ(t) equals λSxt. Similarly, for the given SWxsτ, the time-varying variance λWxτ is given by,
and σ(t) equals λWxt. The integral of SWxsτ over the time domain provides the energy distribution in the scale domain, SWxs.
Let u(t) be a uniformly distributed random variable between 0 and 1 with its marginal CDF denoted as U(u(t)). The relation between u(t) and x(t) can be established based on the probability transformation, U(u(t)) = FX, t(x(t)). The steps in the IPAC algorithm in a pseudo-code form are shown in the flowchart depicted in Fig. 1 and are described as
Fig. 1 Iterative power and amplitude correction algorithm to simulate nonstatinary and non-Gaussian processes (Tf denotes the selected transform in this figure)
Sample {u(jΔ~t)}N~ based on a random number generation algorithm for a uniformly distributed random variable between 0 and 1. Assign {p(j)}N equal to the ascendingly sorted {u(jΔ~t)}N~, and the intensity function yTfs⇀p according to the considered transform pair (Tf(•), ITf(•)), where.
For (Tf(•), ITf(•)) = (FT(•), IFT(•)), yTfs⇀p=yFf=SFxfT, and σt=λFx which is time-independent,
For (Tf(•), ITf(•)) = (ST(•), IST(•)), yTfs⇀p=ySfτ=SSxfτDκf, and σt=λSxt,
For (Tf(•), ITf(•)) = (WT(•), IWT(•)), yTfs⇀p=yWsτ=SWxsτCψs, and σt=λWxt.
II)Iteration
II.1) Sample a sequence of Gaussian white noise, w(t), of length N, calculate eiϕTfs⇀p=ηTfwt;
II.2) Calculate xPCjΔt=ITfyTfs⇀peiϕTfs⇀p, pPCjΔt=FX,jΔtxPCjΔt, and find the rank of pPC(jΔ~t), denoted as rj~, for j = 0, …, N − 1;
II.3) Set xACjΔt=FX,jΔt−1prj, for j = 0, …, N − 1; and calculate eiϕTfs⇀p=ηTfxACjΔt;
II.4) Repeat Steps II.2) to II.3) until the convergence criterion is satisfied.
II.5) x(jΔ~t) = M(jΔt) × xAC(jΔt~).
The algorithm essentially simulates the MODIF process and iteratively corrects the PSD and CDF. The intensity function and the transform pair are used from Steps I) to II.4), while the modulation function only affects the assignment of the final results in Step II.5). Since xAC(jΔ~t) is used in Step II.5), the distribution match (i.e., matching samples of FX(X(t)), {x(jΔt)}N) is ensured by design. One could replace Step II.5) with x(jΔt) = M(jΔt) × xPC(jΔt~) without altering the results if a stringent convergence criterion is employed. As can be observed from the steps of the IPAC algorithm, the analysis, as well as the simulation, is carried out within the same transform pair. It avoids the need to map the obtained results from one type of transform into a different kind of transform (e.g., obtaining the spectrum using continuous WT and then transform it into evolutionary PSD). The algorithm emphasizes the ranking and matching of the probability values. Note that it may be attempting to replace the uniform distribution with the normal distribution for u(t). However, by doing so, it requires the use of the inverse distribution transformation in Steps II.2 and II.3) and increases computing demand.
The algorithm can be simplified if FX, t(x(t)) remains unchanged and only depends on σ(t), that is, the marginal probability distribution of z(t) = x(t)/σ(t), FZ(z(t)), is time-independent and z(t) has zero mean and unit variance. In such a case, we calculate ζjN=FZ−1pjN in Step I); we replace “ pPCjΔt=FX,jΔtxPCjΔt ” in Step II.2) and “ xACjΔt=FX,jΔt−1prj ” in Step II.3) with “ pPC(jΔ~t) = xPC(jΔt)/σ(jΔt) ” and “ xAC(jΔt) = ζ(rj)σ(jΔt) ”, respectively. This avoids the use of probability distribution function during the iteration to gain extra computational efficiency. This simplified version can also be used to generate surrogate for observed {x(jΔt)}N, which has the effect of the modulation function already removed. This is done by calculating {z(jΔt)}N~ = {x(jΔ~t)/σ(jΔt)}N, and letting {ζ(j)}N~ equal to the ascendingly sorted {z(jΔ~t)}N~ in Step I.1) (instead of ζjN=FZ−1pjN), where σ(jΔ~t) is to be calculated based on the PSD function estimated from {x(jΔt)}N~ by using a preferred transform.
The usefulness of surrogate in the context of wind engineering was presented in McCullough and Kareem [28]. The proposed algorithm, when used with WT to generate surrogate, differs from that given in Chavez and Cazelles [4] for testing time-localized coherence, in that the time-varying σ(jΔ~t) is neglected in their algorithm (i.e., the amplitude adjustment is based on {x(jΔt)}N~ rather than its normalized version in the IPAC algorithm). This is convenient and may likely speed up the convergence of the algorithm. However, the basis for the shuffling of {x(jΔ~t)}N~ is unclear if the marginal probability distribution of x(jΔ~t~) for a nonstationary process is assumed to be time-varying.
In this section, we illustrate the proposed algorithm by generating surrogate for a given ground motion record and for a given fluctuating component of wind velocity time history of a high-intensity wind event. We apply the algorithm to sample nonstatinary ground motions for prescribed target PSD, where the target is defined based on a set of ground motion records, and the CDF is assumed to be Gaussian and non-Gaussian. The test of the proposed algorithm for esoteric mathematical models is beyond the consideration of the present study.
Consider the record shown in Fig. 2. By applying ST with the window parameter κ = 1 (see Eq. (6)), the obtained TFPSD function is shown in Fig. 3a, and the calculated time-varying σ(t) is presented in Fig. 3b, showing that the TFPSD varies in time and frequency.
Fig. 2 Ground motions recorded at the CU station, UNAM, Mexico, for the Michoacán earthquake that occurred on September 19, 1985
Fig. 3 Results by using ST for the given record shown in Fig. 2: a TFPSD of the given record, b a generated surrogate and σ(t) of the given record, and c TFPSD of the generated surrogate shown in Fig. 3b
By applying the IPAC algorithm, a surrogate is simulated and shown in Fig. 3b. The TFPSD of this surrogate is depicted in Fig. 3c. The figure shows that the surrogate resembles the given record, and its TFPSD function resembles well that shown in Fig. 3a. As xAC(t) is used for generating surrogate (see Step II.5 in IPAC algorithm), the amplitude (or probability distribution) matching is certain, so no plot is provided. Additional test runs indicate that the convergence is usually achieved within five iterations, depending on the adopted convergence criterion. It was noted that the average TFPSD function from multiple generated surrogates tends to be smoother as compared to the TFPSD of the observed record, which is expected since the observed as well as a single sampled record are realizations of stochastic processes.
The analysis based on ST is repeated but by applying WT using the generalized morse wavelets (GMWs) [33],
where U(ω) is the Heaviside function, aβ, γ = 2(eγ/β)^β/γ^, and β and γ are model parameters. For GMW, Cψ = 2a~2β, γΓ(2β, γ)/(2γ) and Cβ, γ~Γ(β, γ)/(γ). The GMW is an analytical wavelet, and it was used to evaluate the coherence of the seismic ground motions [39]. For the numerical analysis, β = 3, γ = 20, c1ψ = a0 = 0.528, s0 = 2^0.1^ and K = 91 (see Eq.(17)) are employed since these values are suggested as the default values for the algorithm implemented in MATLAB (Version 2019a). The obtained TSPSD and σ(t) of the record are shown in Fig. 4a and b, respectively. A generated surrogate is also shown in Fig. 4b with its corresponding TSPSD function depicted in Fig. 4c. An inspection of the surrogate depicted in Fig. 4b and the original record presented in Fig. 2 indicates that they exhibit similar temporal trends. The TSPSD of the surrogate resembles that of the given record. Again, the convergence is achieved within a few iterations. A comparison of σ(t) shown in Figs. 3 and 4b indicates that they are almost identical. The minor differences between them are due to that ST and WT have different time-frequency (or time-scale) decomposition.
Fig. 4 Results by using WT for the given record shown in Fig. 2: a TSPSD of the given record, b a generated surrogate and σ(t) of the given record, and c TSPSD of the generated surrogate shown in Fig. 3b
Now, consider a wind record presented in Fig. 5a. For simplicity, the box window with a width of 32 samples is used to find the mean wind velocity of the time-varying wind record. By removing the mean, the fluctuating component of the wind is presented in Fig. 5b.
Fig. 5 Wind velocity record from Tower 4 and 10 m height of the rear-flank downdraft that occurred during the evening on June 4, 2002, near Lubbock, Texas [13, 34]: a the wind record, and b the fluctuating component of the record
By applying ST and WT, and following the same analyses that are carried out for the ground motion record shown in the previous section, the obtained results are presented in Figs. 6 and 7. In general, the observations that can be drawn from this example are similar to those presented in the previous section for the ground motion record.
Fig. 6 Results by using ST for the wind record shown in Fig. 5b: a TFPSD of the given record, b a generated surrogate and σ(t) of the given record, and c TFPSD of the generated surrogate shown in Fig. 6b
Fig. 7 Results by using WT for the wind record shown in Fig. 5b: a TSPSD of the given record, b a generated surrogate and σ(t) of the given record, and c TSPSD of the generated surrogate shown in Fig. 7b
Consider a set of 12 ground motion records oriented in the E-W direction for a seismic event that occurred on January 16, 1986, with a local magnitude of 6.1, focal depth of 10.2 km, and an epicentral distance of 25.2 km. These records are recorded by the LSST array in Lotung, Taiwan, where the separation between any two recording sites is less than 100 m, as shown in Fig. 8a. Three records from the 15 recording sites seem to be corrupted and are not considered. The record obtained from FA-1 site is illustrated in Fig. 8b. To minimizing the wave passage effect, first, each of the remaining 11 records is time-shifted with respect to the record presented in Fig. 8b such that the sum of the product of a considered record and that shown in Fig. 8b is maximized after the shift.
Fig. 8 a The LSST array station of selected records, b record at FA1–1 station, c average TFPSD by using ST based on 12 records, and d average TSPSD by using WT based 12 records
It is assumed that the average PSD of the considered 12 records could provide a good representation of the ground motions, at least for such type of seismic event. The calculated average TFPSD based on ST and the calculated average TSPSD based on WT are shown in Fig. 8c and d, respectively. The calculated σ(t) by using the average TFPSD and the average TSPSD presented in Fig. 8c and d are included in Fig. 8b. The obtained PSD and the standard deviation indicate the nonstationarity of the ground motions. σ(t) values obtained by using ST and WT are in very good agreement.
An assessment of the empirical probability distribution of the standardized variable z(t) = x(t)/σ(t) is carried out by considering all 12 records. The empirical distribution of z(t) by considering all 12 records is presented in Fig. 9, indicating that the empirical distribution can be fitted by a normal distribution only for the initial segment of the records. Moreover, the distribution shape is time-varying and non-Gaussian. The non-Gaussian behaviour of the ground motions is supported by the findings given in Radu and Grigoriu [40], indicating that the Gaussian assumption for the seismic ground motions records included in the PEER NGA-West dataset may not always appropriate. However, for this particular set of records, the tail of the distribution is shorter than that of the normal distribution, which differs from the longer tail behaviour suggested by Radu and Grigoriu [40].
Fig. 9 Empirical distributions of the normalized time series of the considered ground motions. a Time interval (0, 2), (2, 4), (4, 6); b (6, 8), (8,10), (10, 12); c (10, 15), (20, 25), (30, 35); d entire duration and the fitted GGD with β
0= 3.01 and β1= 1.54
For illustration purposes, it is assumed that the marginal probability density function of x(t) can be represented by the generalized Gaussian distribution (GGD) [29, 54],
where μ denotes the mean, β0 and β1 are positive model parameters, and Γ(·) denotes the gamma function. If β0 equals 2, it represents the normal distribution. For β0 > 2 and < 2, the distribution tail is lighter and heavier than that of normal distribution. The variance equals β12Γ3/β0/Γ1/β0, and the kurtosis coefficient equals Γ(5/β0)Γ(1/β0)/Γ^2^(3/β0).
By considering β0 = 2 and β1=2 (i.e., standard Gaussian), we use ST and the average target TFPSD function shown in Fig. 8c to sample the records using the IPAC algorithm. Since a comparison of two sampled records is irrelevant for a stochastic process, only a sampled record is illustrated in Fig. 10a. The average TFPSD function obtained from the 1000 sampled records is presented in Fig. 10b, and the calculated spectral acceleration (SA) for a damping ratio of 5% is shown in Fig. 10c for the 1000 sampled records. Similarly, we use WT and the average target TFPSD function shown in Fig. 8d to carry out the simulation. The obtained results are presented in Fig. 10d-f. The PSD functions shown in Fig. 10b and e are almost identical to their corresponding targets presented in Fig. 8c and d. The mean of SA values shown in Fig. 10c and f are in good agreement. The standard deviation of SA obtained by using ST smaller than that obtained by using WT.
Fig. 10 Results based on simulated nonstationary Gaussian records by using ST and WT (the results presented in a to c are for ST, and the results presented in d to f are for WT): a a sampled record based on ST, b average TFPSD of the 1000 sampled record, c SA from 1000 sampled records using ST, d a sampled record based on WT, e average TSPSD of the 1000 sampled record, f SA from 1000 sampled records using WT
We have tested the IPAC algorithm to simulate ground motions for additionally selected target PSD functions. It was observed that in some cases, when the initial or final segment of records has less than 0.5% of total energy, the algorithm may converge very slowly or may not converge. In such a case, it is suggested that such segments with negligible energy are to be truncated. This is in agreement with common practice in earthquake engineering.
To simulate the non-Gaussian process, we consider fX, t(x(t)) shown in Eq. (31) but with β0 = 3.01 and β1 = 1.54 (i.e., a kurtosis coefficient of 2.4) since their use fit the data adequately, as depicted in the last plot in Fig. 9. We repeat the analysis that is carried out for the results presented in Fig. 10. The obtained results for non-Gaussian ground motions are presented in Fig. 11. A comparison of the results shown in Figs. 10 and 11 indicates that the results follow the same trends. To assess the quantitative differences between the obtained SA based on Gaussian and non-Gaussian assumptions, we calculate the ratio of the mean of SA shown in Fig. 11 (i.e., non-Gaussian case) to its corresponding value shown in Fig. 10 (i.e., Gaussian case). We do the same for the standard deviation of SA. The obtained values are presented in Fig. 12, indicating that the mean and standard deviation of SA for the non-Gaussian case with a lighter tail are smaller than those for the Gaussian case for the vibration period less than about 0.5 s. The decrease in SA by considering non-Gaussian excitation is most noticeable for a shorter vibration period. This is because stiffer structures are more sensitive to peak acceleration values. In general, the ratio based on ST is smoother than that based on WT. Note that we refrained from discussing the ratio between the SA values obtained based on ST and WT since such a comparison could be misleading because the TSPSD and TSPSD used are based only on 12 records and from the same seismic event.
Fig. 11 Results based on simulated nonstationary non-Gaussian records by using ST and WT (the results presented in a to c are for ST, and the results presented in d to f are for WT): a a sampled record based on ST, b average TFPSD of the 1000 sampled record, c SA from 1000 sampled records using ST, d a sampled record based on WT, e average TSPSD of the 1000 sampled record, f SA from 1000 sampled records using WT
Fig. 12 Ratio of statistics of SA of simulated records with non-Gaussian and Gaussian distribution assumptions
We note that the simulation of a record (Gaussian or non-Gaussian) described in Figs. 10 and 11 is carried out by using a laptop with Intel(R) Core (TM) i7–7700 CPU @ 3.60GHz (6 core and 12 treads). The wall clock time for simulating a record (including I/O) is less than about 0.5 s if ST is used and is less than about 0.15 s if CWT is used. The difference in the computing time by using ST and CWT can be explained by noting that the octave scale is used in CWT.
We elaborated on the concept of defining a modulated and intensity function adjusted (MODIF) stochastic process in the transform domain. The definitions of the stochastic processes in the transform domain lend themselves to an easily understandable and almost trivial algorithm to simulate stochastic processes. As such a simulated signal may not lead to the prescribed target power spectral density function and marginal cumulative distribution function of the process, we proposed a new iterative algorithm, called iterative power and amplitude correction (IPAC) algorithm, so the sampled signal after the iteration have the prescribed properties. Besides simulating nonstationary and non-Gaussian processes, the proposed iterative algorithm can be used to generate surrogate. The algorithm can be used with Fourier transform, S-transform, and continuous wavelet transforms.
Practical illustrative numerical examples showed the applicability of the proposed algorithm by sampling surrogates for the ground motions and the fluctuating component of winds. The use of the IPAC algorithm to simulate nonstationary Gaussian and non-Gaussian ground motions based on S-transform (ST) and continuous wavelet transform (WT) is presented. The spectral accelerations are calculated using the simulated records. In general, the mean and standard deviation of SA of the simulated records based on ST and based on WT agree well despite the differences between ST and continuous WT and between the frequency-dependent window used in ST and the generalized Morse wavelet used in the continuous WT.
We thank W.X. Zhou for his constructive comments.
H.P. Hong: Conceptualization, design of the work, funding acquisition, analysis, interpretation of data, writing - original draft, review & editing. X.Z. Cui: Design of the work, analysis, interpretation of data, writing - original draft, review & editing, the creation of new software used in the work. Q. Da: Analysis, interpretation of data, writing - original draft. The authors read and approved the final manuscript.
The Natural Sciences and Engineering Research Council of Canada (RGPIN-2016-04814, for HPH). The China Scholarship Council (No.201707980007, for XZC).
The ground motion record shown in Fig. 2 is from http://aplicaciones.iingen.unam.mx/AcelerogramasRSM/Inicio.aspx; the 12 ground motion records from LSST array are from http://www.earth.sinica.edu.tw/. Some or all data, models or code generated or used during the study are available from the corresponding author by request. These include all the coefficients of FT, ST, and WT used in the present study and the simulated time histories.
The authors declare that they have no competing interests.
The ground motion record shown in Fig. 2 is from http://aplicaciones.iingen.unam.mx/AcelerogramasRSM/Inicio.aspx; the 12 ground motion records from LSST array are from http://www.earth.sinica.edu.tw/. Some or all data, models or code generated or used during the study are available from the corresponding author by request. These include all the coefficients of FT, ST, and WT used in the present study and the simulated time histories.