黄天立，尚旭强，程顺，徐兆恒，楼梦麟

(1.中南大学土木工程学院，湖南长沙，410075；2.同济大学土木工程防灾国家重点实验室，上海，200092)

地震动是工程结构地震响应分析和抗震设计的基础。由于记录的具有相同统计特征的实际地震动时程不足，在结构地震响应分析和抗震设计中仍需采用人工模拟生成的地震动时程[1−2]。从20世纪70年代至今，国内外学者提出了多种非平稳地震动的模拟方法，大致可分为4大类，即基于时域调制的方法、基于演变谱的方法、基于自回归滑移平均(ARMA)模型的方法和基于傅里叶相位差谱的方法等[3−9]。需指出的是，这些方法都有一定的局限性和适用范围，如：基于时域调制生成的地震动常使得结构地震响应的计算结果具有较大的误差；基于演变谱的方法合成人工地震动时，难以采用合理、准确的时频谱表示地震动的时频特性；基于ARMA 模型的方法其参数估算比较繁琐，且存在计算量大和不能保证收敛等缺点；基于傅里叶相位差谱的方法仅能用于基岩场地的地震动模拟。

地震动具有明显的强度和频率非平稳特性，对结构和土层的非线性响应具有显著影响[10−14]。为使模拟生成的地震动时程较好地保持实测地震动记录的时频特征，结合各种信号时频分析方法，WANG 等[10,15−18]提出了一些时频完全非平稳地震动的模拟方法。WANG 等[10]采用线调频基信号近似方法模拟生成了具有完全非平稳的地震动，并考虑了地震动频率非平稳特征对单自由度结构和框架结构非线性地震响应的影响。GIARALIS等[15]提出了采用小波变换模拟生成反应谱兼容的时频完全非平稳地震动的方法，并模拟生成了与欧洲抗震规范反应谱兼容的地震动时程。樊剑等[16]采用S变换提出了与抗震规范功率谱相容的非平稳地震动合成方法，发现人工合成的地震波功率谱与目标谱较符合，且生成的地震波较好地继承了实测地震波的非平稳特性。ZHANG 等[17]提出了仅采用Hilbert 变换模拟生成非平稳地震动的方法。何浩祥等[18]在基于小波包分解建立时变功率谱的基础上，建立了包括地震动频域谱和时域谱模型的统一时频谱，并基于统一时频谱人工合成地震动的方法。

Hilbert-Huang变换(HHT)方法具有良好的时频局部化特性，可以更客观地估计地震的时变功率谱，可用于模拟具有时频完全非平稳特性的地震动时程。WEN等[19]采用HHT变换估计得到天然地震动样本的Hilbert 时变功率谱，建立了非平稳随机地震动的模拟方法，并且将该方法应用于多分量随机地震动的模拟[20]。LI等[21]基于Hilbert时变功率谱建立了三向地震动的合成方法，通过HHT 变换调整合成地震动的时频特性，生成了与目标反应谱兼容的地震动时程。

由于HHT 方法中经验模式分解(EMD)得到的固有模式分量(IMF)之间不完全正交，由此获得的Hilbert 时频谱存在能量泄漏问题，进而会影响Hilbert 时变功率谱的估计结果[22−23]。胡灿阳等[23]采用正交化HHT 方法[24]，解决了Hilbert 时变功率谱估计中存在的能量泄漏问题，将实测地震记录估计得到的Hilbert 时变功率谱作为目标谱，提出了基于正交Hilbert 时变功率谱模拟单分量非平稳随机地震动的方法。陈清军等[25]结合演变谱理论和正交HHT 法的优点，提出了一种能同时反映强度和频率非平稳特性的类谐和长周期地震动合成方法。

本文将基于正交Hilbert 时变功率谱的单分量完全非平稳地震动模拟方法进一步拓展，提出了基于正交Hilbert 时变功率谱的多分量完全非平稳地震动模拟方法，该方法在确定描述各分量间相关性的互协方差函数时，假定各分量采用一组相同的随机相位偏差。分别以EL Centro 地震动加速度时程记录南北分量和Northridge 地震Newhall 台站实测三分量地震动加速度时程记录为单、多分量地震动时程模拟对象，模拟生成了具有完全非平稳特征的单、多分量地震动时程；采用Hilbert时频谱、累积归一化Arias 强度指标、加速度反应谱、Renyi熵和分形维数指标验证模拟生成的地震动时程的有效性。

1 基于正交化HHT 方法的时变功率谱估计

1.1 经验模式分解

为使基于Hilbert 变换获得的信号瞬时频率定义具有意义，HUANG 等[26]提出了固有模式函数(IMF)的概念，并给出了一种将任意信号x(t)分解成IMF 分量的方法，即经验模式分解(EMD)方法。EMD 方法通过一种被称为“筛分”的过程，对信号逐步进行分解，最后得到一系列IMF 分量cj(t)(j=1,2,…,n)和1个趋势项rn(t)之和，即

1.2 正交化经验模式分解

楼梦麟等[22]指出，经验模式分解(EMD)所得的固有模式分量(IMF)之间是不完全正交的，由此获得的Hilbert 时频谱存在能量泄漏问题，进而会影响Hilbert时变功率谱的估计。HUANG等[24]进一步提出了针对IMF 分量进行正交化处理的方法，即正交化HHT 方法，解决了Hilbert 时变功率谱估计中的能量泄漏问题，并可应用于地震加速度时程的谱分析。

通过对EMD 分解所得IMF 分量进行正交化处理，可得到完全正交的各阶IMF 分量，详细步骤见文献[24]。式(1)中信号x(t)可被分解成为n个正交IMF 分量c*j(t)(j=1,2,…,n)和1 个趋势项rn*(t)之和，即

式中：c*j(t)(j=1,2,…,n)之间是完全正交的。

1.3 Hilbert时变功率谱估计

对式(2)中任意一个正交IMF 分量(t)进行Hilbert变换，并构造解析信号zj(t)：

信号x(t)可表示为如下形式：

式中：Re(⋅)表示取实部。

将信号x(t)的振幅表示在时间−频率平面上，可得到Hilbert 幅值谱H(t,ω)，称为Hilbert 时频谱，即

式中：ωj(t)为第j个频率点。

Hilbert 时频谱可由EMD 分解得到的IMF 分量或由正交EMD 分解得到的正交IMF 分量计算得到，其描述了信号幅值随时间和频率的变化规律。

考虑信号x(t)和解析信号的关系，可从能量角度估计得到反映信号时频能量特性的Hilbert 时变功率谱f(t,ω)，即

基于正交IMF 分量确定Hilbert 时频谱，并通过式(8)确定Hilbert 时变功率谱，不仅避免了能量泄漏，同时也使采用正交HHT 方法估计得到的正交Hilbert时变功率谱更准确。

2 基于正交Hilbert时变功率谱的单分量非平稳地震动模拟验证

2.1 模拟方法

基于三角级数法模拟地震动，陈清军等[25]提出了结合演变谱和正交HHT方法的地震动模拟方法，其表达式如下：

式中：f(t,ωk)为采用式(8)估计得到的正交Hilbert时变功率谱；Δω=(ωk-ωk-1)为正交Hilbert 时变功率谱的频率间隔；Φk为[0,2π]区间均匀分布、相互独立的随机相位角。该方法模拟得到的地震动时程同时具有较强的强度和频率非平稳特性。

采用正交HHT 方法估计得到的Hilbert 时变功率谱能很好地描绘地震动能量和频率随时间的变化特征，然而，模拟生成的地震动时程的Hilbert时频谱不相同。对于地震动，往往只能实测到1个或者几个时程，因此，可将某实测地震动时程的Hilbert时频谱Hm(t,ω)视为模拟生成的地震动所有时程的Hilbert时频谱的期望值E[H(t,ω)]，即

假设单变量随机过程X(t)可以通过引入1个随机变量来表示：

式中：Φj为独立的随机相位角，均匀分布于区间[0,2π]内。随机过程X(t)的均值μX、协方差KXX(t1,t2)和方差σ2x(t)分别表示为：

由此可见，模拟生成的所有地震动时程，其均值与实测记录的均值相同，其协方差则主要取决于相位角的相关性和在每个瞬时频率的瞬时幅值，方差则由瞬时幅值的平方决定。

选取实测EL Centro 地震动加速度时程南北分量为模拟对象，分别采用HHT和正交HHT方法估计Hilbert 时变功率谱，利用式(9)模拟生成人工地震动加速度时程，这里采用一组完全相同的随机相位角。实测EL Centro 地震动加速度时程南北分量以及分别基于Hilbert 时变功率谱和基于正交Hilbert时变功率谱模拟生成的1条地震动加速度时程如图1所示。从图1可以看出：基于Hilbert时变功率谱和基于正交Hilbert 时变功率谱模拟得到的地震动加速度时程都能很好地反映实测地震动加速度时程记录的强度非平稳特征。

图1 实测和模拟的单分量地震动加速度时程Fig.1 Recorded and simulated single-component earthquake ground motion acceleration time histories

2.2 单分量非平稳地震动时程模拟结果验证

为了演示基于Hilbert 时变功率谱和基于正交Hilbert 时变功率谱模拟生成的单分量地震动加速度时程的准确性和有效性，分别采用Hilbert 时频谱、累积归一化Arias 强度指标、加速度反应谱、Renyi熵和分形维数等指标予以验证。

2.2.1 Hilbert时频谱

实测EL Centro 地震动加速度时程记录南北分量、基于Hilbert 时变功率谱和基于正交Hilbert 时变功率谱模拟生成的单分量地震动加速度时程的Hilbert 时频谱见图2。从图2可以看出：这2 种方法模拟生成的地震动加速度时程，其Hilbert 时频谱都能有效地描述实测地震动时程的时频非平稳特性。

图2 单分量地震动的Hilbert时频谱Fig.2 Hilbert time-frequency amplitude spectra of singlecomponent earthquake ground motions

2.2.2 累积归一化Arias强度指标

累积归一化Arias 强度指标(CNAI)可用于表征地震动时程的累计能量释放特性，其计算式为[27]

式中：0 ≤t≤Td；Td为加速度时程的总持续时间；a(τ)为时间为τ时的加速度。

实测EL Centro 地震动加速度时程记录南北分量、基于Hilbert 时变功率谱和基于正交Hilbert 时变功率谱模拟生成的地震动加速度时程的累积归一化Arias强度指标见图3。从图3可以看出：基于正交Hilbert 时变功率谱模拟生成的地震动加速度时程的CNAI 与实测地震动时程的CNAI 基本一致，两者之间的差异可以忽略不计；而基于Hilbert 时变功率谱模拟生成的地震动加速度时程的CNAI 与实测地震动时程的CNAI 差别较大，相差40%左右。由此表明，基于Hilbert 时变功率谱生成的地震动时程与实测地震动时程相比存在能量泄露现象，这主要是HHT 方法分解得到的IMF分量之间不完全正交所造成的。累积归一化Arias强度指标表明，基于正交Hilbert 时变功率谱模拟生成的地震动加速度时程与实测地震动时程更接近，模拟结果更准确，能更好地表征实测地震动时程的能量释放特性。

图3 单分量地震动的累积归一化Arias强度指标图Fig.3 Cumulative normalized Arias intensities of singlecomponent earthquake ground motions

2.2.3 加速度反应谱

加速度反应谱可作为反映人工地震动加速度时程模拟效果的指标。当阻尼比为5%时，实测EL Centro 地震动加速度时程南北分量、基于Hilbert 时变功率谱和基于正交Hilbert 时变功率谱模拟生成的地震动加速度时程的弹性加速度反应谱如图4所示。从图4可以看出：在周期为0.5～2.0 s 时，基于正交Hilbert 时变功率谱模拟生成的地震动加速度时程的反应谱与实测地震动时程的反应谱之间存在一定差异；当周期大于2 s 时，两者之间吻合较好；在整个周期内，基于Hilbert 时变功率谱模拟生成的地震动加速度时程反应谱与实测地震动时程反应谱之间差别较大。由此表明，采用基于正交Hilbert 时变功率谱的人工地震动模拟方法可获得与实测地震动时程更接近的模拟结果。

图4 单分量地震动的加速度反应谱Fig.4 Pseudo-acceleration response spectra of singlecomponent earthquake ground motions

2.2.4 Renyi熵

Renyi熵是描述系统混乱度的指标，可用于描述信号的复杂度，本文将其作为验证人工地震动模拟方法有效性的一种指标。Renyi 熵Hα由下式确定[27]：

式中：Wψ(a,b)为时频表达系数；α为Renyi 熵的次数，本文取α=3。

实测EL Centro 地震动加速度时程记录南北分量的Renyi熵、基于Hilbert时变功率谱、基于正交Hilbert 时变功率谱生成的50 条地震动加速度时程的Renyi 熵如图5所示。从图5可以看出：基于正交Hilbert 时变功率谱模拟生成的50 条地震动加速度时程的Renyi熵在3.50上下波动，与实测地震动加速度时程的Renyi熵3.42较接近；而基于Hilbert时变功率谱模拟生成的50 条地震动加速度时程的Renyi熵在3.20上下波动，与实测地震动加速度时程记录的Renyi 熵差别较大。此外，由图5还可见50条白噪声信号加速度时程的Renyi熵在5.00上下波动，远远高于实测地震动加速度时程的Renyi熵。

图5 单分量地震动的Renyi熵Fig.5 Renyi entropies of single-component earthquake ground motions

2.2.5 分形维数

分形维数是估计数据自由度和表征数据在数据空间中分布方式的指标，也可用于验证非平稳地震动模拟的有效性。分形维数表征分形的任何维度包括容量维度、相关维数和信息维度，其中，Kolmogorov 容量(盒计数维度)得到普遍应用。通常，分形维数增加表明信息熵增加。本文采用盒计数方法计算地震动加速度时程的分形维数[27]。

设平面内有一个长、宽分别为a和b的分形集M(a,b)，用边长为ε的正方形网格覆盖此分形集，则分形集M(a,b)与正方形网格相交的数量为M(a,b;ε)。定义分形集M(a,b)在标度为ε时的盒维数D为

对于地震动时程，盒维数D表明其在标度为ε时的不规则性，与地震动的频谱分布有关。地震动时程的分形维数越大，表明地震动高频分量越多；反之，表明地震动低频分量越多。本文计算地震动时程的盒维数时，取a=0.08，b=0.1。

实测EL Centro 地震动加速度时程南北分量的分形维数图、基于Hilbert 时变功率谱和正交Hilbert时变功率谱模拟生成的100条地震动加速度时程的分形维数图如图6所示。从图6可以看出：基于Hilbert 时变功率谱、正交Hilbert 时变功率谱模拟生成的地震动加速度时程，其分形维数图与实测地震动时程的分析维数图相差较小，这表明基于Hilbert 时变功率谱、正交Hilbert 时变功率谱模拟生成的地震动时程都能准确地描述实测地震动记录的信号特征。

图6 单分量地震动的分形维数Fig.6 Fractal dimensions of single-component earthquake ground motions

3 基于正交Hilbert时变功率谱的多分量非平稳地震动模拟方法及其验证

3.1 模拟方法

将基于Hilbert 时变功率谱的单分量非平稳地震动模拟方法推广到具有m个分量的多分量非平稳地震动模拟，其关键在于确定各分量之间的互协方差函数。

假设多分量非平稳随机过程Y(t)的第k个分量Yk(t)表示为如下形式：

由式(12)，(13)和式(14)可求得Yk(t)的均值、协方差和方差。多分量非平稳随机过程Y(t)的第p和第q个分量之间的互协方差KYpYq(t1,t2)可表示为[19]

式中：E(⋅)表示取期望；ajp，θjp和ϕjp分别为第p个分量的第j阶IMF 的幅值、相位角和随机相位偏移；akq，θkq和ϕkq分别为第q个分量的第k阶IMF的幅值、相位角和随机相位偏移；rnp和rnq分别为第p和第q个分量的趋势项。

互协方差取决于每一个分量的幅值、相位角和随机相位偏移。若2个分量之间的随机相位偏移ϕjp和ϕkq相互独立，则这2 个分量之间的互协方差为0。通过建立2 个分量之间的互协方差，可确定两分量之间的相关性。本文采用的方法是，使2个分量的两组随机相位偏移相关。2个分量相同序列的IMF的随机相位偏移关系为

式中：ψjq为常数，表示2 个分量相同序列的IMF的随机相位偏移是相关的，而不同序列的IMF 的随机相位偏移是相互独立的。多分量随机过程的均值、协方差和方差均值也可由式(12)，(13)和式(14)求得，每个分量的Hilbert时频谱保持不变。第p和第q个分量之间的互协方差为

式(21)取决于每个分量相同序列的IMF幅值函数、相位角之差和随机相位偏移之差。当相位角之差和随机相位偏移之差，即[θjp(t1)-θjq(t2)+ψjp-ψjq]很小时，余弦函数接近于1，互协方差取决于2个分量的IMF；当[θjp(t1)-θjq(t2)+ψjp-ψjq]接近于π/2 时，余弦函数接近于0，互协方差也接近于0。需指出的是，这仅是模拟多分量随机过程各分量之间相关性的一种方法。

为了模拟多分量随机过程，理论上需要知道各分量之间的互协方差，以此确定式(18)中各分量的随机相位偏移。然而，对于大多数实际多分量地震动记录而言，往往仅有1 个实测记录可以采用，因此，本文假设多分量随机过程各分量相同序号的IMF代表波形，则在同一时间，式(20)中ψjq为0。对于多分量随机过程中的每一分量的模拟，假设采用一组相同的随机相位偏移来模拟，因此，在(21)中，假设ψjp=ψjq=0。

基于这个假定，选择Northridge 地震中Newhall台站记录的垂直(z)和2个水平(x，y)方向分量地震动时程为模拟对象，分别基于Hilbert 时变功率谱和正交Hilbert时变功率谱，利用式(9)和(21)模拟生成多分量地震动加速度时程。

Northridge 地震中Newhall 台站实测的垂直(z)和2 个水平(x，y)方向的地震动加速度、基于Hilbert 时变功率谱和正交Hilbert 时变功率谱模拟生成的1 条三分量地震动加速度时程如图7所示。从图7可以看出：基于Hilbert 时变功率谱和正交Hilbert 时变功率谱模拟得到的三分量地震动时程均能客观地描述实测地震动时程记录，基于正交Hilbert 时变功率谱模拟得到的地震动时程与实测地震动时程差别更小，由此表明，对于多分量非平稳地震动时程的模拟，基于正交Hilbert 时变功率谱的模拟方法比基于Hilbert 时变功率谱的模拟方法更优，模拟生成的地震动时程也能更准确、全面地描述实测地震动时程。

图7 实测和模拟的三分量地震动加速度时程Fig.7 Recorded and simulated tri-directional earthquake ground motion acceleration time histories

3.2 多分量非平稳地震动时程模拟结果验证

为了演示基于Hilbert 时变功率谱和正交Hilbert 时变功率谱模拟生成的多分量地震动加速度时程的准确性和有效性，分别采用Hilbert 时频谱、加速度反应谱、Renyi熵和分形维数等指标予以验证。

3.2.1 Hilbert时频谱

Northridge 地震中Newhall 台站实测的垂直(z)和2 个水平(x，y)方向的地震动加速度时程、基于Hilbert 时变功率谱和正交Hilbert 时变功率谱模拟生成的1条三分量地震动加速度时程的Hilbert时频谱图如图8所示。

从图8可以看出：与模拟生成单分量非平稳地震动时程情况类似，基于Hilbert 时变功率谱和基于正交Hilbert时变功率谱这2种地震动模拟方法模拟生成的地震动时程都能有效地描述实测地震动时程的时频非平稳特性。基于Hilbert 时变功率谱模拟生成的地震动时程，其Hilbert 时频谱的形状虽然与实测地震动时程的Hilbert 时频谱类似，但其峰值与实测地震动时程的峰值相差较大；而基于正交Hilbert 时变功率谱模拟生成的地震动时程不仅在Hilbert 时频谱的形状上与实测地震动时程的Hilbert 时频谱吻合，而且其峰值与实测地震动时程的峰值相差很小。由此表明，基于正交Hilbert 时变功率谱的多分量地震动模拟方法改进了基于Hilbert 时变功率谱的多分量地震动模拟方法，特别是在反映实测地震动时程记录的时频能量特征方面，基于正交Hilbert 时变功率谱方法模拟生成多分量地震动时程能更好地描述实测地震动时程的时频能量特征。

图8 三分量地震动的Hilbert时频谱Fig.8 Hilbert time-frequency amplitude spectra of tridirectional earthquake ground motions

3.2.2 加速度反应谱

当阻尼比为5%时，Northridge 地震中Newhall台站实测的垂直(z)和2 个水平(x，y)方向的地震动加速度时程、基于Hilbert 时变功率谱和正交Hilbert时变功率谱模拟生成的1条三分量地震动加速度时程的加速度反应谱如图9所示。

从图9可以看出：基于Hilbert时变功率谱模拟生成的三分量地震动加速度时程的加速度反应谱与基于正交Hilbert 时变功率谱模拟生成的三分量地震动加速度时程的加速度反应谱显著不同；基于正交Hilbert 时变功率谱模拟生成的三分量地震动加速度时程弹性反应谱与实测地震动时程弹性反应谱更接近。由此表明，基于正交Hilbert 时变功率谱的多分量非平稳地震动模拟方法，可获得与实测地震动时程更接近的模拟结果。

图9 三分量地震动的加速度反应谱Fig.9 Pseudo-acceleration response spectra of tri-directional earthquake ground motions

3.2.3 Renyi熵

Northridge 地震中Newhall 台站实测的垂直(z)和2 个水平(x，y)方向的地震动加速度时程和分别基于Hilbert 时变功率谱、正交Hilbert 时变功率谱模拟生成的50条三分量地震动加速度时程的Renyi熵如图10所示。从图10可以看出：基于Hilbert时变功率谱模拟生成的50 条地震动加速度时程的Renyi 熵低于实测地震动加速度时程的Renyi 熵；基于正交Hilbert 时变功率谱模拟生成的50 条地震动加速度时程的Renyi熵与实测地震动加速度时程的Renyi 熵更接近。由图10还可见：50 条白噪声信号加速度时程的Renyi熵远远高于实测地震动加速度时程的Renyi熵。

图10 三分量地震动的Renyi熵Fig.10 Renyi entropies of tri-directional earthquake ground motions

3.2.4 分形维数

Northridge 地震中Newhall 台站实测的垂直(z)和2 个水平(x，y)方向的地震动加速度时程记录的分形维数图、基于Hilbert 时变功率谱、正交Hilbert时变功率谱模拟生成的100条三分量地震动加速度时程的分形维数图如图11所示。表1所示为实测和基于Hilbert 时变功率谱、正交Hilbert 时变功率谱模拟得到的Northridge-Newhall 地震垂直(z)和2个水平(x，y)方向加速度时程的分形维数。

从图11和表1可以看出：基于Hilbert 时变功率谱和正交Hilbert时变功率谱这2种模拟方法得到的多分量非平稳地震动时程的分形维数与实测地震动时程的分形维数之间相差都比较小，这表明基于正交Hilbert 时变功率谱模拟生成的多分量地震动时程和基于Hilbert 时变功率谱模拟生成的多分量地震动时程都能准确地描述了实测地震动记录的信号特征。

表1 实测和模拟的Northridge-Newhall地震三分量加速度时程的分形维数Table 1 Fractal dimensions of measured and simulated Northridge-Newhall tri-directional earthquake acceleration time history

图11 三分量地震动的分形维数Fig.11 Fractal dimensions of tri-directional earthquake ground motions

4 结论

1)基于多分量随机过程中各分量过程的模拟可采用一组相同随机相位偏差的假定，确定了描述多分量非平稳地震动各分量间相关性的互协方差函数。将基于正交Hilbert 时变功率谱模拟单分量完全非平稳地震动模拟方法进一步拓展，提出了基于正交Hilbert 功率谱的多分量完全非平稳地震动模拟方法。

2)分别以EL Centro 地震动加速度时程记录南北分量和Northridge 地震Newhall 台站实测三分量地震动加速度时程记录为单、多分量地震动时程模拟对象，模拟生成了具有完全非平稳特征的单、多分量地震动时程。采用Hilbert 时频谱、累积归一化Arias强度指标、加速度反应谱、Renyi熵和分形维数等指标验证了模拟生成的地震动时程的有效性。

3)与基于Hilbert 时变功率谱生成的单、多分量地震动时程相比，基于正交Hilbert 时变功率谱生成的单、多分量非平稳地震动时程能更真实地描述实测地震动记录的强度和频率非平稳特性。