时频非平稳地震动模型及结构随机振动应用*

何浩祥， 范少勇， 闫维明

(北京工业大学工程抗震与结构诊治北京市重点实验室 北京,100124)

引 言

地震动是由震源释放的地震波引起的地面运动，具有复杂而强烈的时间-空间随机性和非平稳性[1-2]。即使在同一次地震过程中，不同区域地震动的时频特性也有所差别。对地震动时频两域的非平稳性进行深入研究能够充分了解地震动的演变机制及随机性，并可为地震动模拟提供重要的评价依据。此外，关于地震动特性的研究是土木工程抗震设计与动力分析的基础，建立能够全面反映地震动非平稳性的模型可以更准确地计算结构随机响应，为结构精细化随机振动分析和可靠性评价提供坚实的基础。

由于通常无法在一次地震中同一区域取得多次地震动记录，传统的随机理论通常假定每一条地震动记录均符合各态历经假定，即认为可用时间的平均来代替集合的平均[2-3]，这导致具体地震动的确定性特征和随机成分没有被准确分离，且地震动时频域非平稳性被显著弱化。如果随之进行的结构随机振动分析对非平稳性考虑不充分，则结构随机响应将更不准确。在传统的随机振动研究中通常采用均匀调制方法反映地震动的时域非平稳性，即将真实地震动等效为一个零均值的平稳随机过程(一般为白噪声)和一个随时间渐变的确定性强度包络函数的Hadamard乘积(各元素点乘)[1-3]。常用的强度包络函数主要有连续型和分段型两类[2,4]。采用均匀调制方法得到的地震动时程包络一般均为单峰型曲线，并不能准确刻画地震动幅值变化的丰富性和复杂性，其频域分布具有明显的平稳性，不能反映地震动频域的强非平稳特征和能量分布特点。对于多阶自振周期与地震动卓越周期接近的结构，采用均匀调制地震动进行随机振动分析有可能严重影响结果的准确性。

为了解决上述不足，Saragoni等[5]提出通过时间段内的Gamma函数调制过滤高斯白噪声的方法，并提出了一种时域和频域完全非平稳随机地震动模型(包含强度和频率非平稳)。Lin等[6]提出用散粒噪声、过滤白噪声或三角级数叠加信号等替换均匀调制方法中的白噪声。Sanaz[7]提出采用强度包络函数和时变过滤白噪声模型共同模拟地震动时频域非平稳特性。在具体应用以上方法时存在难以精确确定频带能量分布的局限，虽然在人工地震动生成中易于实现，而在真实地震动随机分析中难以应用。梁建文[8]在演化谱理论的基础上，采用一个谱表示方法来模拟频率非平稳地震动过程，但此方法依赖目标谱和权重函数。Conte等[9]构建了基于sigma振子过程的完全非平稳模型，并可由真实地震动演变功率谱经自适应最小二乘法拟合确定模型参数，虽然其合成精度相对较高但计算复杂，且对时频两域的非平稳细节模拟仍不充分。此外，部分研究者建议在计算结构响应的统计特征时直接将相应的输入白噪声功率谱改为传统的地震动功率谱工程学模型[10]，如Kanai-Tajimi模型及其系列修正模型等，但由于相关的功率谱模型对频率非平稳性的表征依然不够精细[11-12]，且该处理方法在理论上也不严谨，因此不宜作为有效方法。

因此，针对上述问题开展新型时频域非平稳地震动模型和相关随机振动分析方法具有重要的理论与工程意义。笔者从地震动的传播和演变机制以及随机特性上进行探究，提出能够更准确全面反映地震动时频域非平稳特性的强度包络模型，并改进传统的均匀调制方法，进而建立反映时频非平稳关联特征的结构随机振动分析方法。

1 非均匀调制地震动模型

描述地震动的物理和数学模型众多[4,13]，从应用性的角度可分为地震动工程应用模型和地震动随机过程模型两种，前者强调地震动模型及其合成信号的特性应满足指定的反应谱或功率谱的需求，而后者更注重于地震动时频非平稳性的描述，并主要应用于结构随机振动分析领域。笔者的研究主要针对地震动随机过程模型开展。

在传统的均匀调制地震动随机模型中，并不考虑频率非平稳性，并通常将地震动表示为

A(t)=G(t)n(t)

(1)

其中：A(t)为非平稳地面加速度过程；G(t)为强度包络函数；n(t)为一零均值的平稳随机过程。

在上述模型中，强度包络函数通常是单峰正值函数，无法准确全面地刻画地震动多次波动的多峰现象，也不能体现地震动正负幅值的差异，因而不能充分反映地震动在时域上的随机性和非平稳性。更重要的是，强度包络函数并不能明显改变平稳随机过程在频域上的能量分布，因此上述模型不能反映频域非平稳特性，更无法直接转化为有效的速度和位移时程。增强上述模型对时频域非平稳性的表征能力需要同时对G(t)和n(t)的内涵和形式进行改进，相关的改进可以从地震动的传播和演变机制的研究成果中得到借鉴。

在地震发生时，震源所产生的尖脉冲经地层介质的吸收后，其波形将会被拉长，这种具有时延的单个反射波称为地震子波。在地震学和地震勘探中，褶积模型理论认为地震波形是地震子波动态叠加的结果[14-15]。相关研究表明：从地下许多反射界面发生反射时形成的地震子波，其振幅取决于反射界面反射系数的绝对值，极性的正负决定于反射系数的正负, 到达时间的先后取决于界面深度和覆盖层的波速。实际的地震记录A(t)可视为由地震子波w(t)和反射系数γ(t)褶积组成，相关公式如下

A(t)=w(t)γ(t)

(2)

比较式(1)和式(2)可发现，二者的数学形式一致但物理意义不同。在褶积模型中，地震子波是振荡衰减型低频随机信号，具有多峰性，幅值也呈现正负交替，这与均匀调制模型中的强度包络函数有明显区别。此外，褶积模型中的反射系数是零均值频域非平稳噪声信号，这与均匀调制模型中的白噪声也有所区别。地震子波和反射系数的非平稳特性及其褶积特征能较充分反映地震动的时频域非平稳性，这为均匀调制模型的改进提供了良好的借鉴。

在诸如褶积模型的地震动分解处理中，地震子波和反射系数依然是复杂随机信号，可以采用信号处理技术对其进行更深入的分析，而滤波方法是判断信号频域非平稳特征的主要手段之一。大量地震动频谱分析结果表明，地震动的频域能量主要分布在0～10Hz。由于地震动信号采样率较高，频谱分析中最高频率通常超过25Hz，因此可以通过滤波实现地震信号中重要成分和次要成分的分离。然而，采用诸如低通滤波器的理想滤波技术对信号的处理过于严格和刻板，缺乏灵活性和调控能力。相反，采用移动平均算法对信号进行平滑处理，相当于对信号进行了低通滤波，且相应的滤波器具有从通带到阻带平缓过渡的特点，因此既能分离局部随机波动的特征又能凸显整体变化规律。文献[16]提出地震动时域信号可分解为多峰平滑信号和非平稳噪声信号两部分。多峰平滑信号主要反映了地震产生机理和演变规律的基本特征，虽然由于知识不完备性难以阐释其细节变化，而非平稳噪声信号则反映了地震动传播过程中受到的扰动和强随机性。在此基础上，文献[16]建立由多峰平滑时域信号As(t)和非高斯时域噪声信号An(t)构成的地震动时域模型A(t)

(3)

非高斯时域噪声信号An(t)具有如下表达形式

An(t)=As(t)|αNn(t)|

(4)

其中：Nn(t)为归一化的高斯白噪声；α为其调幅系数。

由式(3)和式(4)可得到地震动时域模型的综合表达式为

(5)

对比式(2)和式(5)可发现，地震动时域模型A(t)与褶积模型在表达形式和物理内涵上均类似，因此可以改进均匀调制模型的不足，为精确细致反映地震动非平稳特性提供了必要的处理方法和分析工具。

上述地震动时域模型是由信号能量平衡条件来约束，虽然精度较高但在结构随机振动分析中不易处理，因此笔者建议将式(3)和式(5)改为如下的非均匀调制形式

(6)

其中：Acn(t)表示均值为1的白噪声。

虽然式(6)相当于忽略了原地震动时域模型A(t)的能量表达式中的交叉项，但由于As(t)和An(t)均为随机性较强的信号，其乘积相对较小，且可通过再次的平滑处理和调整调幅系数α来实现进一步的修正，因此依然能较准确地反映原有地震动信号的随机特性。

此外，多峰平滑时域信号As(t)可以通过以下两式计算得到

其中：A0(t)为原始地震动信号；As0(t)为初步平滑信号；S函数为平滑函数，在Matlab程序中可采用smooth命令实现；n和m均为滑动窗宽度或需要考虑的数据点数；sgn为符号函数。

由于As(t)需要反映地震动本质特性，因此试算时可不断增加滑动窗宽直到其幅值稳定为止。在具体分析中，n和m建议取值区间分别为[150,200]和[5,30]，同时需要选取适当的调幅系数α,使As(t)的最大峰值与原波A(t)最大峰值相同。

为了验证笔者提出的非均匀调制地震动模型的精确性，应用相关方法对Taft(N21E)地震加速度信号进行分解和合成，并与传统的均匀调制地震动模型进行比较，相应的结果如图1所示。对原信号进行平滑处理时，式(7)和式(8)滑动窗宽度n和m分别取为200和5，优化后的调幅系数α为2.25。此外，由于已经知道Taft波的具体信息，因此在均匀调制合成方法中没有采用常用的强度包络函数模型，而是通过直接对地震波绝对值进行窗宽为500的平滑处理并归一化获得。

图1 不同地震动模型加速度合成效果对比

由图1可看出，多峰平滑包络信号As(t)具有往复振荡的特点，且最大值不必等于1，这与传统的时域强度包络信号有明显的区别。此外，噪声信号Acn(t)的幅值相对较大。采用非均匀调制模型合成的地震动与真实地震动在时域上十分接近，而按照均匀调制模型合成的地震动则与真实值有明显差别。为了对比不同模型在频域上的精度，对不同信号进行了功率谱分析，如图2和图3所示。

图2 不同地震动模型的加速度功率谱对比

图3 非均匀调制模型各部分频谱分布

结果表明，由非均匀调制模型合成的地震动的频谱分布也与真实地震动接近，而均匀调制模型的频谱分布较均匀且不能反映卓越频率特征。从图3也可看出，多峰平滑包络信号充分保留了真实地震动的主要频谱成分，是非均匀调制地震动模型精确性的根本保障。为进一步验证非均匀调制地震动在低频信号合成方面的能力和精度，对Taft(N21E)地震位移信号进行合成，滑动窗宽度n和m分别取为150和26，优化后的调幅系数α为4.0，时域和频域结果如图4所示。由图可见，非均匀调制模型依然具有良好的精度，而均匀调制模型对低频能量的表征精度进一步下降，该结论与前人研究结果吻合[2]。

综上所述，由于平滑处理后的多峰平滑包络信号具有正负振荡特性以及足够的细节刻画能力，地震动的主要频谱成分被充分提取，加之噪声信号的强随机性，由二者合成的地震信号能够精确表征原信号在时域和频域的非平稳特征。因此，非均匀调制地震动模型既可以根据原波产生大量的地震动随机信号，从而为基于Monte Carlo方法的结构随机分析提供可靠的随机样本，也可以为基于功率谱的结构时域随机振动分析提供更精确的激励模式。

图4 地震动位移信号合成效果对比

2 全非平稳地震动模型

目前计算演变谱的方法虽然较丰富，但精度通常不高。小波包分解方法采用正交小波将信号分解成各尺度上的正交小波包分量，再对各分量用其相应的小波包基函数为窗函数进行时频变换，其分解机理明确，由此得到的时变功率谱从理论表达上是准确的，且分辨率高。文献[16]的研究表明，将小波包分解产生的Paley序子频带序列按升频进行重新排列，并采用正交性和紧支撑性的离散Meyer(dmey)小波作为小波包基函数可以更准确地描述模拟地震记录信号时频变化特征。

基于上述研究并参考非均匀调制地震模型的构建模式，笔者提出采用二维平滑技术的全非平稳地震动演变谱模型A(ω,t)，该模型由多峰平滑时频域谱As(t)和时频域非高斯噪声谱An(t)构成，表达式为

A(ω,t)=As(ω,t)+An(ω,t)=

As(ω,t)[1+αNn(t)]=As(ω,t)Acn(t)

(9)

其中：As(t)和Acn(t)的维数均为p×q，p为频率向量的维数，q为时域向量的维数。

为了验证上述全非平稳地震动演变谱模型的精确性，采用离散Meyer小波函数对Taft(N21E)地震加速度信号进行9层小波包分解，从而获得地震动演变谱，如图5(a)所示；利用稳健补偿最小二乘计算方法[19]对该演变谱进行二维平滑处理，总滑动窗宽度为2，获得的多峰平滑演变谱如图5(b)所示；生成调幅系数α为2.50的随机噪声谱，如图5(c)所示；最后按式(9)合成地震动演变谱，如图5(d)所示。为了便于观察，上述图采用谱绝对值表示原波频谱和多峰平滑演变谱等结果。

从结果可看出：平滑演变谱与原波演变谱在总体变化规律上类似，但幅值较小，合成后的演变谱与原波演变谱十分接近，能同时反映时频域的非平稳变化特征。分别计算真实演变谱和合成演变谱的时域和频域的边缘条件，即分别求演变谱在时域和频域的叠加值，结果如图6和图7所示。由图可知，全非平稳合成地震动的频域合成精度较高，而时域合成精度略低，主要原因在于目前的二维平滑算法不易实现时频两域同精度的协调平滑，该问题可以通过改进算法得到进一步解决。

图5 全非平稳地震谱分解与合成

图6 全非平稳地震动时域信号

此外，还可基于小波包分解分别获得由均匀调制模型和非均匀调制模型生成的地震波的演变谱，如图8所示。结果表明，非均匀调制模型地震波演变谱与原波演变谱吻合较好，其频域上的非平稳性也较准确细致。均匀调制模型地震波的时频分布与真实波有较大差异，在地震动模拟和结构分析时需谨慎使用。为了进一步比较以上不同地震动模型的合成精度，笔者选取了Ⅰ～Ⅳ场地的10条地震波进行相关分析，并以合成偏差率作为精度评判指标。时域合成偏差率定义为原波与合成波各点幅值平方差的累积和与原波幅值平方和的比值；频域合成偏差率定义为原波与合成波的功率谱幅值差的累积和与原波功率谱的比值。均匀调制模型地震动在时域和频域上的平均合成偏差率分别为0.292和0.304；非均匀调制地震动的平均合成偏差率分别为0.076和0.069；全非平稳地震动的平均合成偏差率分别为0.231和0.091。可见，非均匀调制地震动模型和全非平稳地震动模型的精度均比传统的均匀调制模型高，在地震动合成或结构随机振动分析中采用这两种模型将有效提升计算精度。

图7 全非平稳地震动功率谱

图8 均匀与非均匀调制模型时频谱对比

3 基于全非平稳谱模型随机振动分析

由于笔者提出的非均匀调制模型和全非平稳模型的形式已与均匀调制模型有明显差别，因此需要在传统随机振动分析理论的基础上进行适当推广。

若采用非平稳随机激励f(t)表示为前述多峰平滑时域信号和均值为1的白噪声的点乘形式，则在该改进均匀调制非平稳随机激励f(t)的作用下，初始静止的线性体系响应可表示为

(10)

其中：h为脉冲响应函数；τ为0～t的任一时刻。

由于多峰平滑包络的均值为零，故激励和体系响应均具有零均值，且其协方差函数为

Ky(t1,t2)=E[y(t1)y(t2)]=

Ry(t1,t2)-E[y(t1)]E[y(t2)]=Ry(t1,t2)

(11)

响应y(t)的协方差函数和自相关函数Ry(t)相等，且有

(12)

根据Wiener-Khinchin定理，对于宽平稳随机过程Acn(t)有

E[Acn(τ1)Acn(τ2)]=Rx(τ1-τ2)=

(13)

其中：Scn为非均匀调制模型中噪声Acn(t)的自功率谱密度函数。

因此，式(11)可表示为

(14)

其中

(15)

令t1=t2=t，则可得到响应y(t)的方差为

(16)

其中：σy为响应的标准差。

响应y(t)的自演变谱密度为

Sy(ω,t)=I*(ω,t)I(ω,t)Scn(ω)=

(17)

对于多自由度线性结构，在非均匀调制地震激励A(t)=As(t)Acn(t)的作用下，其随机动力方程为

(18)

令Y(t)=Φq(t)，其中Φ为结构振型向量，将其带入上式并在方程两侧同乘ΦT，则有

-γ*M*As(t)Acn(t)

(19)

其中：M*=ΦTMΦ;C*=ΦTCΦ;K*=ΦTKΦ;γ*=ΦTMI/M*。

考虑到结构参数矩阵与振型的正交性，上式可表示为n个单自由度方程

(20)

该系统的响应Y(t)的方差为

(21)

若由于结构自由度数较多或计算效率偏低，可考虑采用虚拟激励法计算结构随机响应[20-21]。基于虚拟激励法构造一个虚拟的外部激励，则式(18)变为

(22)

仍可采用振型分解法得到该系统的响应Y(t)的方差，有

(23)

(24)

上述分析方法是针对非均匀调制激励的，当采用完全非平稳地震模型时，需要按式(9)计算全非平稳地震动演变谱模型A(ω,t)，并按下式计算Ii

(25)

其他计算公式类同。由于本研究提出的非平稳地震动随机模型的精度主要取决于多峰平滑信号的时频域特性，加之实际结构的动力分析需要考虑地震动传播条件和场地特征的影响，因此对于具体结构的随机分析应选取相应场地的典型地震动或满足工程要求的人工生成地震波进行计算，并确保其多峰平滑信号具有代表性。

4 算例与分析

为了进一步验证非平稳地震动模型及其在结构随机振动分析中的有效性，笔者建立了层数分别为1～20的20个二维剪切型钢筋混凝土框架结构并进行随机振动分析。为便于比较不同结构的动力响应规律，各个结构的层高均为3 m，层重均为5.2×108kg。假定各个结构的无阻尼基本周期为0.06n(n为相应总层数)，则可算出各层的等效剪切刚度，进而建立刚度矩阵。假定各结构阻尼比均为5%，则可构建Rayleigh阻尼矩阵，最终按式(16)建立结构动力学方程。将Taft波作为水平向地震动输入，按照随机振动方法计算结构随机地震响应。为了研究地震动合成模型对结构响应的影响，笔者基于均匀调制模型、非均匀调制模型和全非平稳模型分别对以上20个具有不同自振周期结构的标准差进行计算分析。基于均匀调制模型和非均匀调制模型的5层结构的1层、3层和5层的位移标准差如图9所示，其他典型结构顶层位移标准差结果依次如图10所示，5层结构的顶层位移演变谱绝对值如图11所示。

图9 结构典型楼层位移标准差

图10 典型结构顶层位移标准差比较

图11 结构的顶层位移演变谱绝对值

从以上结果可以看出，基于非均匀调制模型和全非平稳模型的结果较接近，且其波动规律与地震动幅值的变化规律相关，表明其随机响应能够较细致地反映地震动时域非平稳性对结构的影响。此外，当基本周期较小时基于上述两种模型的结构标准差随时间变化剧烈，随着周期增大其变化趋缓，反映了地震动频域非平稳性的影响。与此相反，由于均匀调制模型反映时频域非平稳性的能力均明显不足，相应的结果通常变化平缓且最大值偏小，不能充分反映地震动非平稳性对结构的影响。

为了进一步对比3种方法结果的差异并探究其原因，提取基于均匀调制模型和全非平稳模型的所有结构顶层位移标准差最大值，并与基于非均匀调制模型的结果相比获得最大值的比值，根据结构基本频率(基本周期倒数)列出上述结果，如图12所示。

图12 结构标准差最大值比值规律

可以发现，均匀调制模型最大值均小于非均匀调制模型最大值，且其变化规律与原地震动的功率谱外形相近。这是因为均匀调制模型地震动功率谱的频率能量分布是均匀的，而非均匀调制模型地震动的频谱特征能够充分反映地震动频域非平稳性，因而与真实功率谱吻合，而二者的差异仍然具有真实功率谱的特征，最终导致按照均匀调制模型计算短周期结构(基本周期小于0.3s)时很可能严重低估结构响应。此外，如果地震动的长周期成分较显著或结构属于大型复杂结构时，由均匀调制模型得到随机响应也可能偏低，不利于结构的安全设计。由全非平稳模型得到的结构响应与非均匀调制模型结果较接近，但由于其在时频两域的非平稳特征相对不够精确，因此建议作为参考方法。

综上，可认为按照传统均匀调制模型计算结构随机响应不能充分获得结构时频域的非平稳特性，结果偏于不安全，基于笔者提出的非均匀调制模型的结构随机响应计算方法能够较好地解决上述问题。基于全非平稳模型的结构响应计算方法在原理上更严格而全面，但由于现有多维平滑技术的协调性和精度略低，因此需要进行更深入的研究和改进。

5 结 论

1) 地震动在时域和频域均具有复杂的随机性和非平稳性，传统的地震动随机过程模型通常采用均匀调制模型，并不能全面细致地表征地震动时频非平稳性。笔者根据地震褶积模型理论以及之前的研究成果，提出由多峰平滑信号和非平稳噪声信号点乘而成的非均匀调制地震动模型。该模型能够在时域上充分反映地震波的波动和衰减特性，其频谱分布也与真实地震动十分接近，因此可以同时表征时频域非平稳性。在此基础上，根据演变谱原理和需求提出了以时间和频率为变量的全非平稳地震动建模方法。分析结果表明，笔者提出的两种模型在时域和频域的合成精度上均优于均匀调制模型。

2) 针对非均匀调制模型和全非平稳模型的形式和特点，在传统随机振动分析理论的基础上对结构随机振动响应求解方法进行了改进。通过算例证明了按照均匀调制模型计算结构随机响应不能充分获得结构时频域的非平稳特性，且结果偏于不安全，而利用笔者提出的两种模型可以获得更精确的结构响应概率特征和更细致的非平稳性演变过程。本研究主要针对具体地震动特性开展，在今后研究中可根据不同场地条件和工程需求获取具有典型非平稳特征和概率统计特征的地震波形，实现结构随机响应特征的分析和预测。