引言

重要工程结构的抗震设计需要一定数量的且满足一定地震环境及场地条件的地震动加速度时程作为结构动力反应分析的输入；基于随机振动理论的结构地震反应分析也需要满足同一集系特征的地震动时程作为输入． 由于强震观测记录有限，不能满足上述地震工程的理论研究及实际应用中诸多方面的需求，因此，以随机过程理论为基础、利用数值方法合成满足给定地震动工程特性的人工地震动，即地震动模拟或仿真，一直是地震工程领域理论研究和工程应用中的热点问题(李英民，赖明，1998a，b，c)． 地震动模拟应尽可能全面地反映实际天然地震动的工程特性，其合成的人工地震动应尽可能全面地模拟实际天然地震动对工程结构的影响效应．

由于地震震源的复杂破裂过程、地震波在地球介质中传播的频散现象以及地表土层非线性的力学行为，地震动通常具有非平稳特性，而且地震动的非平稳特性对工程结构的地震反应具有重要的影响(Kubo，1993；Hasgür，1995；杨红等，2001；Wang，2002；陈健云等，2003；曹晖，林学鹏，2006)． 因此，研究人员提出了多种非平稳地震动的模拟方法，具体可分为4大类，即基于时域调制的方法(Scanlan，Sachs，1974; Rezaeian，Der Kiureghian，2008)、基于渐进谱的方法(Yeh，Wen，1990；Conte，Peng，1997；梁建文，2005a，b；Liang et al，2007；张翠然，陈厚群，2008)、基于自回归滑移平均(ARMA)模型的方法(Chang et al，1982；Cakmak et al，1985；Deodatis，Shinozuka，1988；Conte et al，1992)和基于傅里叶相位差谱的方法(杨庆山等，2001)． 各种方法均有其局限性和适用范围(李英民，赖明，1998a，b，c；Zhang et al，2012)． 基于希尔伯特变换的方法(Zhang et al，2012)则通过建立给定天然地震动瞬时相位的随机过程模型，实现对非平稳天然地震动的模拟． Zhang et al(2012)利用希尔伯特谱和Morlet小波谱，验证了该方法生成的地震动样本具有给定的天然地震动的非平稳特性，其时域波形的形态和能量的时频分布特征与天然地震动相同． 为了进一步研究该方法在模拟天然地震动工程特性方面的能力，本文将基于单自由度体系的弹塑性动力响应，深入验证该方法在模拟实际天然地震动峰值特性及其对工程结构影响效应方面的效果．

1.   基本方法

基于希尔伯特变换的非平稳地震动模拟方法基本流程如下：

1)给定天然地震动记录x₀(t)．

2)利用希尔伯特变换求出时程x₀(t)的瞬时幅值 a₀(t)和解卷后的瞬时相位θ₀(t)如下：

式中P表示取广义积分的Cauchy主值．

3)利用傅里叶变换求出瞬时相位θ₀(t)的傅里叶幅值谱A_θ，0(f)和相位谱Φ_θ，0(f)：

式中Re和Im分别表示取复数的实部与虚部，i表示虚数单位，i²=-1．

4)利用汉宁(Hanning)窗对Φ_θ，0(f)进行低通滤波，

式中Hann[x(f)]表示对函数x(f)进行低通滤波处理后得到的低频分量．

5)假定Φ_θ，n(f)为随机过程Ψ(f)的一个样本，并且假定随机过程Ψ(f)为确定性函数A(f)与随机过程Θ(f)乘积的形式，即

其中确定性函数的形式如下：

式中Θ(f)假定为区间[-1, 1]均匀分布的随机频率过程，且满足： ① 对任意频率f，Θ(f)为区间[-1, 1]均匀分布的随机变量；② 若f₁≠f₂，则Θ(f₁)与Θ(f₂)为相互独立的随机变量．

6)按照下式构造与原始天然地震动x₀(t)对应的随机过程x(t)：

其中

7)设定不同的随机数，利用式(12)即可生成一系列与原始地震动具有相同非平稳特性的地震动样本x_s(t)．

上述方法的具体论述详见Zhang et al(2012)． 通常该方法能够合理地模拟地震动的高频成分，但该方法在生成地震动样本时会引入一些额外的无意义的低频成分；另一方面，该方法假定地震动的低频成分为物理意义明确的确定性变量，进而通过建立高频成分的随机过程模型实现对非平稳地震动的模拟． 因此，利用该方法生成的地震动样本x_s(t)通常要保留原始天然地震动样本x₀(t)的频率低于截止频率f_c的低频成分． 应当指出的是，由于震源机制、地震波传播路径及局部场地条件等因素的差异，天然地震动中物理意义明确的低频成分对应的截止频率亦不尽相同． 在实际工作中，应根据期望保留的低频成分的频带范围，合理地设定截止频率f_c．

图1给出了原始El Centro波(El Centro台站获取的1940年美国Imperial Valley地震的南北分量记录，矩震级M_W6.95，震中距12.99 km)的加速度、速度和位移波形；利用上述方法生成的一个人工地震动样本的相应波形如图2所示，可以看出，该方法生成的人工地震动样本与原始天然地震动样本具有相似的加速度、速度和位移波形．

图  1  原始El Centro波的波形

Figure  1.  The waveform of original El Centro recordings

下载: 全尺寸图片 幻灯片

图  2  模拟El Centro波的一个地震动样本的波形

Figure  2.  The waveform of one ground-motion sample simulating the El Centro wave

下载: 全尺寸图片 幻灯片

利用该方法共生成出与图1所示天然地震动相应的1 000个人工地震动样本，针对每个样本确定其峰值加速度、峰值速度和峰值位移，进而确定每个地震动参数分布的直方图，如图3所示． 其中，横轴为地震动参数的取值区间(共40个)，纵轴为对应于每个取值区间的地震动样本数量，图中还标示出了原始天然地震动的峰值． 可以看出，该模拟方法生成的地震动样本的峰值大部分位于原始天然地震动的峰值附近． 因此，人工地震动样本能够较好地模拟实际天然地震动的峰值特性，包括峰值加速度、峰值速度和峰值位移特性．

图  3  模拟El Centro波的地震动样本峰值的分布(a)峰值加速度;(b)峰值速度;(c)峰值位移

Figure  3.  Distribution of the peak values of ground-motion samples simulating the El Centro wave(a)Peak ground acceleration；(b)Peak ground velocity;(c)Peak ground displacement

下载: 全尺寸图片 幻灯片

2.   方法验证

本文将利用理想单自由度体系的弹塑性动力响应，验证上述基于希尔伯特变换的非平稳地震动模拟方法在模拟天然地震动对结构影响效应方面的效果． 所涉及到的单自由度体系的动力响应参数如下所述．

单自由度体系在地面运动作用下动力响应的控制方程为

式中： m为振子的质量; c为体系的阻尼系数(采用黏滞阻尼理论)；f_s(u)为体系的恢复力，是振子相对地 面的位移u(即体系的变形)的函数; ü(t)，(t)和u(t)分别为振子相对地面的加速度、速度和位移；ü_g(t)是地面运动的加速度；振子的初始速度和初始位移均为零． 对于理想弹塑性体系而言，恢复力和变形之间的关系如图4所示，u_y为体系的屈服位移，f_y为体系的屈服强度． 在变形达到u_y之前，体系为线弹性体系，其刚度为k；当变形超过u_y，体系为理想弹塑性体系，受输入地震动过程的影响，体系刚度在不同时刻取k或0． 针对理想弹塑性体系，定义与其相应的线性体系(corresponding linear system，简写为CLS)(Chopra，2001)，该体系的质量和阻尼与弹塑性体系一致，其刚度为弹塑性体系屈服前的刚度． 弹塑性体系的“固有”周期和阻尼比分别定义为其CLS的“自振”周期和阻尼比，即

图  4  恢复力-变形关系

Figure  4.  Relationship between restoring force(f_s) and deformation(u)

下载: 全尺寸图片 幻灯片

设在某一给定地震动作用下，弹塑性体系CLS的最大变形反应为u₀，相应的恢复力为f₀，如图4所示． 定义弹塑性体系的屈服强度降低因子(Chopra，2001)为

在给定地震动作用下，R_y越大，体系的动力响应就越早进入到弹塑性状态． 记体系弹塑性变形的最大值为u_m，体系的延性系数定义为

针对给定的地震动加速度ü_g(t)，若给定弹塑性体系的固有周期T、阻尼比ζ及屈服强度降低因子R_y，则首先计算自振周期为T、阻尼比为ζ的CLS的最大相对位移反应u₀，然后根据式(17)确定弹塑性体系的屈服位移u_y=u₀/R_y，进而利用常规数值积分方法(如Newmark法)求解式(15)，即可确定与参数R_y相对应的体系弹塑性动力响应． 相应地，仿照弹性反应谱的定义，可以定义地震动ü_g(t)的弹塑性反应谱.

分别称为弹塑性绝对加速度反应谱、相对速度反应谱和相对位移反应谱，分别表示在给定屈服强度降低因子R_y的条件下，给定地震动引起的不同体系的弹塑性绝对加速度反应、相对速度反应和相对位移反应的最大值． 除上述体系最大响应之外，地震动对体系的输入能量以及体系自身的阻尼耗能和弹塑性滞回耗能也是描述地震动对结构影响效应的重要参数． 在式(15)两边分别乘以du并积分可得

上式等号左边3项依次为振子动能E_k、体系阻尼耗能E_d和体系弹塑性滞回耗能E_h，等号右边为地震动对体系总的输入能量E₁． 针对给定地震动，这些能量是体系自身特性(以参数T、ζ和R_y描述)的函数，且有

分别称为弹塑性输入能量谱、阻尼耗能谱和滞回耗能谱，其中T_d为地震动总的持续时间． 式(20)—(22)与式(24)—(26)定义了地震动相应于屈服强度因子R_y的弹塑性反应谱． 另外，地震动相应于延性系数μ的弹塑性需求谱定义如下． 对于给定的地震动加速度ü_g(t)，使得弹塑性体系(T，ζ)的延性系数达到μ所需的体系屈服位移为u_y，则可定义地震动的弹塑性需求谱为

针对给定地震动ü_g(t)、参数T、ζ和μ，可以采用迭代的方法确定上述弹塑性需求谱(Chopra，2001)． 弹塑性需求谱的工程意义在于，结构在给定地震动作用下，其塑性变形发展程度(即其破坏程度)达到指定要求(以μ描述)所需的结构屈服位移． 对于弹性体系，μ=1，不允许其发生塑性变形，那么其屈服位移就必须大于或等于在给定地震动作用下结构弹性相对位移反应的最大值． 此外，对于μ=1，即弹性体系，式(27)定义的地震动需求谱等同于其弹性相对位移反应谱，因为地震动在某一控制周期处的弹性相对位移反应谱可理解为在该地震动作用下，使得相应体系保持弹性所需的最小屈服位移． 利用式(27)定义的需求谱可定义伪速度需求谱如下：

定义了上述地震动的各类弹塑性谱之后，即可利用这些参数描述地震动对弹塑性体系的影响效应，进而利用其对上述模拟方法进行验证．

针对图1所示的原始El Centro波，利用基于希尔伯特变换的模拟方法生成与其具有相似非平稳特性的30个地震动样本，计算所有地震动样本与屈服强度因子R_y=3相对应的弹塑性绝对加速度反应谱、相对速度反应谱、相对位移反应谱、输入能量谱、阻尼耗能谱、滞回耗能谱，以及分别与延性系数μ=1(弹性)和μ=4相对应的伪速度需求谱(由式(28)定义)，具体计算结果如图5所示． 基于希尔伯特变换的地震动模拟方法将给定的天然地震动视为某一随机过程的一个样本或一次实现，通过建立地震动瞬时相位的随机过程模型实现对该地震动随机过程的求解． 从图5中可以看出，作为同一随机过程的不同样本，人工地震动样本能够很好地模拟实际天然地震动样本在弹塑性反应谱、输入能量谱、阻尼耗能谱、滞回耗能谱以及弹塑性需求谱方面的特征． 而上述几类谱均从不同角度描述了地震动对结构的影响效应，因此，基于希尔伯特变换的地震动模拟方法生成的人工地震动样本能够较好地模拟实际天然地震动样本对结构的影响效应．

图  5  原始El Centro波与人工地震动样本频谱的比较

(a) Absolute acceleration response spectrum (R_y=3);(b) Relative velocity response spectrum (R_y=3);(c) Relative displacement response spectrum (R_y=3);(d) Input energy spectrum (R_y=3);(e) Damping dissipated energy spectrum (R_y=3);(ƒ) Hysteretic dissipated energy spectrum (R_y=3);(g) Pseudo-velocity dem and spectrum (μ=1);(h) Pseudo-velocity dem and spectrum (μ=4).Thick black lines denote spectra of original ground motion,and thin grey lines denote spectra of simulated ground-motion samples

Figure  5.  Comparison of the spectra of original El Centro wave with those of artificial ground-motion sample

下载: 全尺寸图片 幻灯片

图6给出了TCU052台站记录到的1999年台湾集集地震的南北分量地震动波形(矩震级M_W7.62，震中距39.58 km)． 该记录的一个最大特征是其速度波形曲线出现了一个显著的低频脉冲，同时导致地表出现了永久性的位移，如图6中位移曲线所示． 基于希尔伯特变换的模拟方法生成的与图6所示记录相应的一个地震动样本如图7所示． 该模拟方法主要针对随机特征明显的地震动高频成分进行随机建模，实现对天然地震动的模拟，在此过程中保留了具有较明确物理意义的地震动低频成分． 对于图6所示天然地震动而言，其速度和位移波形均受地震动低频成分控制，因此，即使该方法生成的地震动加速度与原始天然加速度波形在高频成分方面存在一定的差异，但速度和位移样本与原始天然地震动相似，如图6和图7所示． 图6所示的天然地震动样本与30个人工地震动样本不同类型弹塑性谱之间的比较如图8所示． 从该图中可以看出，人工地震动样本较好地模拟了实际天然地震动的工程特性．

图  6  原始TCU052台站记录到的1999年台湾集集地震的南北分量地震动波形

Figure  6.  Waveforms of original TCU052 recordings obtained in 1999 Taiwan Jiji(Chi-Chi)earthquake

下载: 全尺寸图片 幻灯片

图  7  模拟TCU052台站记录的一个地震动样本的波形

Figure  7.  Waveforms of one ground-motion sample simulating TCU052 recordings

下载: 全尺寸图片 幻灯片

图  8  原始TCU052记录频谱与人工地震动样本频谱的比较

(a) Absolute acceleration response spectrum (R_y=3);(b) Relative velocity response spectrum (R_y=3); (c) Relative displacement response spectrum (R_y=3);(d) Input energy spectrum (R_y=3);(e) Damping dissipated energy spectrum (R_y=3);(ƒ) Hysteretic dissipated energy spectrum (R_y=3);(g) Pseudo- velocity dem and spectrum (μ=1);(h) Pseudo-velocity dem and spectrum (μ=4).Thick black lines denote spectra of original ground motion,and thin grey lines denote spectra of simulated ground-motion samples

Figure  8.  Comparison of the spectra of original TCU052 recordings with those of artificial ground-motion sample

下载: 全尺寸图片 幻灯片

3.   结论

非平稳地震动模拟方法生成的地震动样本不仅需要在时域波形曲线的形态以及能量的时频分布特征等方面与原始天然地震动相似，而且对于工程而言，所生成的样本还需要如实地模拟实际天然地震动对结构的影响效应． 基于希尔伯特变换的非平稳地震动模拟方法通过对实际天然地震动瞬时相位的随机建模，实现了对天然地震动非平稳特性的模拟，并且最初的研究也论证了该方法在模拟天然地震动时域波形曲线的形态以及能量的时频分布特征方面的效果(Zhang et al，2012)． 作为对该方法的进一步验证，本文基于理想弹塑性单自由度体系的动力响应，验证了该方法在模拟天然地震动对结构影响效应方面的效果． 结果表明，该方法生成的地震动样本能够较好地模拟实际天然地震动的工程特性． 在地震工程的理论研究和工程实践中，可以利用该方法合成出与给定天然地震动具有相同特性的一系列地震动样本，进而采用蒙特卡洛(Monte Carlo)方法开展结构随机地震反应的数值模拟研究，从而为地震动工程特性的研究和结构可靠度分析提供基础．