地震动全非平稳时频谱统计模型及其在结构分析中的应用

何浩祥，丁佳伟，闫维明

（北京工业大学工程抗震与结构诊治北京市重点实验室，北京100124）

引言

地震动是复杂的非平稳随机过程，其随机性主要来自于震源机制、传播途径和局部场地条件，且时域和频域的非平稳性都比较明显。地震动的基本特性可以通过振幅、频谱和持时三要素来描述。该三要素从时域频域两方面刻画了地震动特征，是影响结构响应的主要因素。通常采用强度包络函数来描述地震动时域幅值的变化特征，常用的数学模型主要分为单峰函数型［1-2］和分段函数型［3］两类。单峰模型不适用于描述多峰型地震动波形，分段连续模型对于描述冲击型地震记录不恰当，并且在分界点处不光滑。由于不同地震动持时和峰值出现时刻各不相同，直接对大量地震动进行统计分析难以获得简洁直观的模型，因此上述模型都是根据部分典型地震动波形和观察经验确定的，并不直接具有明确的统计特征。如何确定具有统计意义的地震动强度包络函数模型具有重要的研究意义。此外，地震动在频域上也具有显著的非平稳特性，通常采用功率谱来描述其频谱特性，以表征地震动在不同频率成分上的能量分布。常用的多种功率谱模型均以Kanai-Tajimi［4］谱为代表的过滤白噪声型功率谱为基础进行适当调整以抑制低频或高频成分［5-6］。但以上功率谱模型均为光滑的单峰曲线，与真实地震动功率谱差别较大，对非平稳的表征不够充分和细致。综上所述，传统时域或频域理论模型不能充分地反映地震动时频域的非平稳性。

时频分析是一种将时域信号映射到时频空间描述和表达的方法，在信号局部或瞬态特征刻画方面具有独特的优势。由现代时频分析方法获取的时变功率谱能够从时间、频率和幅值三个方面来共同描述非平稳随机过程的时频局部特性［7-9］。目前计算时变功率谱的方法主要有短时Fourier变换、WVD（Wigner-Ville Distribution）分解、经验模态分解以及小波变换和小波包变换等［10-11］。虽然时频分析方法众多，但大部分方法在处理地震动信号时并不能充分满足时域和频域的边缘条件，仍然不够精确。文献［12］通过对比分析确认了具有合适基函数的小波包分解方法可准确描述地震动的时频变化特征，是建立时变功率谱（演变谱）准确而有效的工具。然而，目前尚缺乏反映各类场地条件影响的具有统计意义的时变功率谱模型或精准的经验公式，这制约了时变功率谱在地震动特性分析中的深度应用。

此外，目前的结构随机振动分析理论虽然通常根据演变谱理论进行求解，但一般都进行了简化处理。采用的地震动模型通常为均匀调制模型和传统非均匀调制模型，前者没有考虑频域非平稳性而后者对非频域平稳性的反映不精确、不充分，影响了结构地震随机振动分析的精度。

有鉴于此，本文提出了基于能量持时的地震动强度统计模型建立方法，获得了不同场地类型下的地震动强度包络统计模型。在此模型基础上提出了多峰型功率谱模型和全非平稳时变功率谱模型，并根据非平稳信号的边缘条件和能量非负性验证了其精确性，为建立精确的地震动时变功率谱统一模型并进行精确的结构随机响应分析提供坚实基础。结构随机振动分析算例结果表明：采用传统的随机振动分析方法的结果偏保守且不够准确，而基于精确时变功率谱（演变谱）计算得到的结构响应能够准确而细致地反映时频两域非平稳性的影响，适宜在地震动分析和结构动力分析中推广。

1 基于能量持时的地震动强度统计模型

地震动在时域上呈明显的强度非平稳性，可采用随时间变化的加速度强度包络函数来描述振幅变化过程。目前强度包络函数模型主要分为单峰型函数模型和多段式强度包络函数模型［3］。前者的典型公式包括：

式中t为时间；I0，a，b和c均为特征参数。

在多段式强度包络函数模型中，三段式模型最为常用，其表达式为

式中a，b，c均为模型参数，t1和t2分别为控制强震平稳持续段的首末时刻。

传统的强度包络函数模型都是依据一定的典型地震记录或经验建立的，并未从波形统计特征方面进行确切验证，具有一定的局限性，需要根据大量数据建立具有统计特征的强度包络曲线。在获取地震动强度包络统计值时，由于不同地震记录下的峰值出现时刻和平稳持时段均不同，如果直接将每条地震波的绝对值的最大值对齐，叠加后的总波形平稳段将交叉混叠，并不具有真实性和代表性，也无法应用传统的强度包络函数来表示。

实际上，强度包络函数的定义中隐含着对累积能量及持时的规定，能量持时定义为地震动加速度能量累积达到两个不同阈值之间的时间差，其中能量累积过程由地震动加速度平方的积分表示［13-15］。能量持时可按下式计算

式中H为Heaviside函数，Ab和Ae分别为第一个和第二个阈值，A（t）为归一化的Arias强度，计算公式如下

式中a（τ）为地震动加速度记录，Td为地震动总持续时间。

基于强度包络函数模型实际上是对累积能量持时的表达这一观点，本文提出基于能量持时的地震动强度统计模型建立方法，可以将大量地震记录的强震平稳持时段集中汇集在一起并使强度模型具有统计特征。该方法的具体步骤为：（1）将各类场地大量的地震记录进行归一化处理；（2）通过公式（4）和（5）分别计算不同场地类型下每条地震记录累积能量的起始阈值点tb及终止阈值点te；（3）计算每条地震记录起始与终止阈值点的中心时刻tm；（4）确定各类场地中所有地震记录中心时刻tm的最大值tmax，并以此为基准将各条地震记录的中心时刻tm与tmax对齐，视为同一时段发生的地震记录；（5）由于各地震记录的时程长度不同，将地震记录的起点时刻补零找齐，使地震记录时程具有相同长度，并构成整体叠合记录；（6）对上述叠合记录取绝对值并归一化，最终得到具有统计意义的等效波，之后可用拟合方法确定偏差最小的强度包络函数。上述方法的特点在于：可将地震动的主要能量（强震平稳段）集中在有限时域范围内，为评价能量统计特征提供良好的基本条件；避免了不同地震记录的持时不同造成的不便；能量对齐之后地震动峰值仍然在适当的时域范围内，在此基础上进行峰值分析、平滑处理和建立强度包络更有统计意义。

为了验证上述方法的有效性，针对四类场地按照如下的要求从太平洋地震工程研究中心（PEER）数据库中选波：1）所选地震震级应大于3.5级；2）地震波的PGA大于0.05g小于0.8g；3）震中距大于10 km小于240 km；4）平均土层剪切波速符合规范的要求；5）一次地震发生最多采用2条不同地点的地震记录；6）均采用水平向的记录。分别选取了世界范围内的100条典型的地震记录进行统计分析并建立强度包络模型。分析结果证明5%到75%的累积能量持时段包含了地震的强震的主体［13］，是震动最强和破坏力最强的部分，因此通过式（4）和（5）计算每条地震记录的5%和75%的能量阈值。以Ⅱ类场地和Ⅳ类场地的统计结果为例，如图1所示。可以看出，按照本文方法获得的基于能量持时的地震动时程曲线既保留了每条地震记录的变化趋势，又保证了总体的模型特征，细致而全面地表征了不同场地条件下的地震波强度演变特征。该方法克服了最大峰值对齐方法带来的波形多峰和交错现象。同时，相关结果表明基于能量持时的地震动强度统计波形确实呈现出与传统强度包络函数模型一致的变化趋势。分别采用单峰型函数模型和三段式函数模型对四类场地的统计地震记录进行拟合，并计算拟合误差，结果如表1所示。可见误差均很小，采用两种模型均具有较高的精度，因而传统模型确实具有良好的普适性。根据本文方法获得的强度时程曲线和相关包络函数模型准确有效，并具有统计意义，可以为今后的地震动时频域非平稳性分析和结构动力分析提供良好的数据支持。

图1 典型场地地震记录强度包络函数模型Fig.1 Envelope function model of seismic record intensity for typical sites

表1 各模型与四类场地地震记录的相对误差/%Tab.1 Relative error between each model and four types of site seismic records/%

2 多峰功率谱统计模型

大量地震动特性研究表明：地震动的强度非平稳性及其对结构动力分析的影响固然重要，但地震动的频域非平稳性更加复杂，对结构随机动力响应的影响也更深远。基于上文由统计分析得到的基于能量持时的强度时程曲线，本文将开展多峰功率谱统计模型的研究。

地震动的功率谱在频域上都是复杂多变的，具有强烈的非平稳性。长期以来，为了分析方便，研究者通常使用以Kanai-Tajimi功率谱及其修正谱为代表的过滤白噪声模型来描述地震动功率谱，这些谱模型均呈现平滑单峰形状，不足以精确表征真实功率谱多峰现象和非平稳性，导致地震动的非平稳特性过度简化且以功率谱模型为基础的结构随机振动分析的精度不足。为了解决上述问题，本文采用多峰功率谱函数［16］作为频域的表达，如下式所示

式中S（ω）为多峰功率谱函数，ω为地震动圆频率分量，n为等效土层数，ωgj和ζgj分别为第j个土层的固有频率和阻尼比，S0为白噪声激励的功率谱密度。上述模型采用Yamada等［17］提出的由多个二阶线性滤波器响应之和组成的表达式，反映了地震动通过场地内多个土层的效果，可以精细模拟多峰现象及复杂变化特征，比其他功率谱模型具有更高的精度和准确性，是建立时变功率谱模型的关键。对上文基于能量持时的地震动强度统计模型建立方法得到的各类场地100条地震波分别取均值，最终得到各类场地的典型归一化地震波。根据各类场地代表性归一化地震记录的真实功率谱确定多峰平滑功率谱的参数，Ⅰ类至Ⅳ类场地的等效土层数量n依次为18，12，14和16。通过最小平方和误差的多元函数的目标优化方法可识别功率谱函数模型中的ωgj，ζgj和S0等参数。

为了验证多峰功率谱函数模型的精度，首先计算得到各类场地的传统功率谱的谱值，并以此作为真实值，分别采用Kanai-Tajimi谱、Clough-Penzien谱和欧进萍谱对各类场地的结果进行基于满足最小平方和差的参数拟合，得到相应的等效功率谱，并与多峰功率谱进行对比，其中Ⅱ类场地的对比结果如图2所示。各功率谱与真实功率谱的相对误差和确定系数（R-square）如表2所示。

表2 各功率谱与真实功率谱的相对误差及确定系数Tab.2 Relative errors and R-square of power spectrum and real power spectrum

图2 不同功率谱精度比较Fig.2 Accuracy comparison of different power spectra

从以上结果可以看出：传统的单峰功率谱与真实值之间存在较明显差异，仅能粗略地表征频域能量分布的基本趋势，并不能精确展现频域能量的非平稳性及多峰特征。此外，Kanai-Tajimi谱和欧进萍谱没有抑制低频能量的分布，这与实际不符，虽然Clough-Penzien谱在低频处的表征较为合理，但总体上仍偏大于真实功率谱。多峰功率谱与真实功率谱十分接近，全面细致地反映了地震动非平稳特征，且自动避免了单峰谱在低频和高频精度上的缺陷。因此，本文提出的多峰功率谱统计模型具有较高的精度，更适合在复杂地震动分析中使用。

3 时变功率谱（演变谱）的要求及特性

地震动在时域和频域均具有强非平稳性，因此研究地震动的时频全非平稳性更加全面，也更能反映地震动同时在时频两域的演变过程及其对结构响应的影响。利用上文由统计分析得到的基于能量持时的强度时程曲线，下文将开展时变功率谱（演变谱）的统计模型研究。

针对类似地震动这样的强非平稳随机过程Y（t），Priestley［18］在Wiener拓广的广义Fourier积分的基础上提出了演变谱模型，其意义在于描述演变随机过程在某一时刻邻近的“瞬间”平均能量在频域上的分布，从而可以用通常的“频率”与“能量”概念来描绘非平稳随机过程的谱特性。演变谱可以表示为

式中G（ω）为一个正交增量过程，A（t，ω）为t与ω的确定性的调制函数。根据正交增量的性质，Y（t）的协方差函数可以表示为

当t1=t2时，可得Y（t）的方差为

故定义非平稳随机过程Y（t）的演变谱密度SY（t，ω）如下

式中SX（ω）为零均值的平稳白噪声随机过程X（t）的功率谱，可表示为

式中S0为常数。若设S0=1，则非平稳随机过程Y（t）的演变谱密度SY（t，ω）可以简化为

由上述分析可知：对大量地震动进行分析，准确获取每一记录的时变功率谱（演变谱），并由此寻求和建立与实际地震动时频分布统计特征一致的非平稳地震动模型具有十分重要的理论和工程意义。为了确保相关演变谱密度函数的准确性，除了要满足上述条件外，还应通过信号能量守恒要求及边缘条件来校核。Cohen［5］根据Parseval定理定义了非平稳信号Y（t）时域和频域的边缘条件，表达式如下

式中SY（t，ω）为非平稳信号Y（t）的时变功率谱，也即演变谱。SFT（ω）为信号的Fourier变换幅值谱。从上式可以看出：信号某一特定时间的所有频率能量分布之和等于其瞬时能量；信号的某一特定频率的能量分布在全部时间上的总和等于能量密度频谱。精确的时变功率谱（演变谱）要满足上述的边缘条件。随着时频分析研究的进展，研究者提出了关于时频分布的其他的要求，如时频分布必须是实的，时频分布关于时间t和f的积分应该与非平稳信号的总能量E相等，即如下式

从上述的时频分布的要求可以看出，边缘分布的特性和非负性一起保证时频分布准确反映信号的谱能量、瞬时功率和总能量。边缘条件要求可以保证信号的总能量（如平均时间、平均频率、时宽和带宽）正确给定。非负性则可以进一步保证分布的期望是切合实际的物理解释。非负性和正确的边缘分布特性在一起可以保证时频分布的强有限支撑性。

国内外研究者相继建立了不同的地震动时变功率谱（演变谱）特性描述方法和模型。Okaya等［19］将平稳和瞬态波谱估计的多重滤波技术拓展到演变谱估计，但由于模型识别参数过程复杂，且精细程度与阻尼比经验取值有关，限制了其广泛应用。Nakayam等［20］利用窄带和低通滤波，直接生成复调制函数从而建立演变谱，但精度有限。Conte等［16］提出了全非平稳时变功率谱函数模型，其假设地震加速度为一个由多个均匀调制过程叠加而成的Sigma振子过程，时变功率谱为各过程的功率谱密度函数的叠加，具有多峰性。该方法在假定地震记录为平稳高斯过程的叠加存在较大的问题，且模型并不完全满足时频域边缘条件。林家浩等［21］、梁建文等［22］、胡灿阳等［23］和刘章军等［24］曾分别提出由强度包络函数和功率谱函数相乘而形成的时变功率谱模型，但其调节能力和多峰性有限，并不能充分反映时频域复杂的时变特征，同时也不满足边缘条件，在研究和工程应用方面存在不足。在进行结构随机振动分析时，以上模型所得的结构响应与真实结构的结果也存在较大差异。何浩祥等［12］提出采用时域和频域归一化边缘值的误差标准差来量化和评判时变功率谱的精度，并对短时傅里叶变换、WVD变换、小波变换和小波包变换等时频分析方法的精度进行了对比，结果表明离散Meyer小波函数的小波包变换具有较好的精度，且满足边缘条件及非负性要求，可用来获得地震动时变功率谱。然而，上述研究是针对具体地震动进行的，并没有根据大量地震动的分析结果建立具有统计意义的时变功率谱模型。为了进一步研究地震动的时变功率谱的基本规律和通用模型，本文对基于能量持时的地震动强度统计模型建立方法得到的各类场地典型波进行小波包分解，从而获得相应的时变功率谱，小波包基函数为9层离散Meyer小波函数，结果如图3所示。通过验算证明上述的时频谱严格满足边缘条件。与地震动时域能量时程或传统功率谱相比，真实时变功率谱的相对误差如表3所示。由表可见其误差非常小，可认为相关时变功率谱能够从时间、频率和幅值3个方面准确而全面地描述非平稳随机过程的时频特性。

表3 时变功率谱模型相对误差/%Tab.3 Relative error of time-varying power spectrum model/%

图3 不同场地下的地震动真实时变功率谱（演变谱）Fig.3 Actual time-varying power spectra（evolution spectra）of ground motions in different fields

经过对地震动真实时变功率谱的深度拟合分析，选用类似式（1）的适用性强的时域强度包络函数和式（6）的多峰功率谱，本文提出如下时变功率谱等效函数模型

式中a1，a2，b1，b2，c1和c2均为模型的拟合参数。四类场地具有不同的土层模拟值，四类不同场地的时变功率谱结果如图4所示。

从图3的结果可以看出：地震动时频功率谱是十分复杂多变的，时频分布具有强烈的突变性和随机性。随着场地类型的提高，地震动高频分量所占比例逐渐增长。结合图4的结果来看，本文提出的时变功率谱等效模型保留了真实时变功率谱的主体特征，剔除了部分噪声并可较精确地表现多峰现象和时频域的非平稳基本特性。

图4 不同场地下的地震动时变功率谱等效模型Fig.4 Equivalent model of time-varying power spectrum of ground motion in different fields

仅仅依靠时变功率谱本身并不能校核其精确性，为了验证该时变功率谱等效模型的精度，需要检验其是否满足边缘条件及非负性。利用式（13），将真实时变功率谱的时域或频域投影与时变功率谱等效模型的时域或频域投影作对比并进行误差分析，Ⅰ类至Ⅳ类场地的对比结果如图5所示。此外，与地震动时域能量时程或传统功率谱相比，等效时变功率谱的相对误差如表3所示。通过对比分析，可以认为等效模型误差在允许的范围内，因此本文提出的时变功率谱等效模型可以较精确地表征地震记录在时频两域的非平稳性，且符合时频分解的边缘条件特征及非负性，故等效模型具有工程应用价值。

图5 不同场地的时变功率谱时频域边缘值对比Fig.5 Comparison of time-varying power spectrum marginal values in time-frequency domain at different sites

应用本文提出时变功率谱函数模型，根据式（12）可以推得时变调制函数A（ω，t）的函数表达式如下

其中演变功率谱的强度调制函数参数如表4所示。

表4 时变功率谱模型强度调制函数参数取值Tab.4 Parameter values of intensity modulation function of time-varying power spectrum model

综上所述，本文提出的时变调制函数具有良好的精度和适应性，可以实现对时变功率谱的曲面的精细调节和多峰模拟，从而更加细致地刻画真实的时变功率谱，适合在研究和工程应用中推广。

4 基于演变谱的结构全非平稳动力分析

应用地震动时频功率谱（演变谱）模型可以在随机微分方程均方解的意义下求解结构响应时变过程及其概率特征［25-26］。设结构受时频域全非平稳随机激励引起的响应y（t）表示为

则任意两个响应量的互相关函数为

利用正交增量的性质和函数积分的性质得到：

若使t1=t2=t，式（19）右端的被积函数就是y1（t）与y2（t）的互功率谱函数，即

使t1=t2=t，双非平稳随机激励引起的响应y（t）的方差为

假设外部激励的自谱密度为Sxx（ω），由此可以构造一个虚拟的确定性的外部激励

此激励在t时刻产生的响应为

若令y1=y2=y，可以得到任意响应的自功率谱函数计算式如下

5 算例分析

为了进一步验证本文提出的时变功率谱函数模型在结构随机振动分析中的有效性，对一个七层的二维剪切模型框架结构进行随机振动分析。每层结构层高为3 m，每层均重5.42×108kg，结构阻尼比为5%，基本周期为0.42 s。为了便于比较各种非平稳模型和计算方法的区别，将Ⅱ类场地的典型地震波作为水平向地震动输入，按照前述方法计算结构地震响应。本文计算所用到的功率谱函数模型如图2所示，比较当采用真实功率谱（按传统方法计算得到的功率谱）、时变功率谱、等效时变功率谱、欧进萍功率谱、Clough-Penzien功率谱及Kanai-Tajimi谱的情况下，分别在非平稳强度包络曲线（原始地震动时域包络）与光滑强度包络曲线（对原始时域包络进行单峰函数拟合）下进行结构随机响应分析。考虑全非平稳时频谱情况下，功率谱为时变功率谱或等效时变功率谱，其结构位移标准差按式（22）计算，其中一层和七层的位移标准差曲线如图6和7所示。当功率谱为真实功率谱、欧进萍功率谱、Clough-Penzien功率谱及Kanai-Tajimi谱时，可采用一般的演变谱简化计算方法，即式（20）中的A（ω，t）被简化为时域强度包络曲线g（t），相应结果如图6和7所示。

图6 非平稳强度曲线包络底层与顶层部位移标准差比较Fig.6 Comparison of standard deviations between enveloped bottom and top layers of non-stationary strength curve

从以上结果可以明显的看出：根据时变功率谱和等效时变功率谱计算的结构响应与根据真实功率谱计算的结构非常接近，均可视为精准结果，可根据研究需求合理选用。基于传统过滤白噪声单峰模型功率谱的结构响应显著大于全非平稳条件下的结果，其原因在于单峰型光滑功率谱的谱值均大于真实功率谱，没有真实细致刻画功率谱的变化特征。此外，结构的响应与功率谱和强度包络形状密切相关，功率谱的幅值和形状直接影响结构位移响应的幅值，而强度包络函数形式直接影响结构位移响应曲线的变化趋势。因此，在结构随机振动分析时，若对结构响应的幅值（极值）精准度要求较高，需要选择更精确的多峰功率谱模型；若对响应时程的变化趋势精准度要求较高，需要着重选择变化细致的强度包络模型。

为了进一步的对比不同方法的差异，建立层数分别从1到25的25个二维剪切型钢筋混凝土框架结构并进行随机振动分析，假定各结构的无阻尼基本周期为0.06n（n为相应总层数），阻尼比均为5%。提取基于各方法的所有结构顶层位移标准差最大值，与真实功率谱的计算结果相对比如图8所示。

图8 具有不同周期的结构响应精度比较Fig.8 Comparisons of structural response accuracy with different periods

可以看出：对于不同的周期结构，时变功率谱与传统功率谱的计算结果均十分相近，可视为精确结果。等效时变功率谱的结果一般有±10%之内的偏差，是可以接受的。而基于传统单峰功率谱的结果均远大于真实功率谱的结果，最大偏差在60%以上。在结构设计和分析中，可认为根据以上单峰谱计算的结构响应是不精确或偏于保守的，不宜在此基础上再考虑安全系数。

综上所述，本文提出的时变功率谱函数模型及其等效模型，不仅满足时频谱的边缘条件要求，并且更加真实、准确、细致地对地震动全非平稳性进行了描述。经过在结构随机振动分析中与基于传统单峰功率谱的演变谱计算方法相比较，进一步验证了本文的时变功率谱函数模型的精确，且可以在工程实践中推广。

图7 光滑强度曲线包络下底层与顶层部位移标准差比较Fig.7 Comparison of standard deviation of bottom and top displacement with smooth strength curve

6 结论

传统的地震动时域强度包络函数均为经验函数，缺乏足够的统计意义。本文提出了基于能量持时的地震动强度统计模型建立方法，该方法保留了各条地震记录的变化趋势、主体特征和样本特性。经过验证，该方法适用于各类场地的大量地震记录的分析，为研究地震动强度的统计特性提供了基础。针对传统单峰模型功率谱偏于简化和保守导致精度不足的局限，多峰功率谱可以细致表征有效频率范围内功率谱幅值的变化趋势，适用于对地震分析精度要求较高的地震动研究及大型复杂结构设计。鉴于传统时频分析得到的时变功率谱缺乏精确的统计模型和解析表达式，本文基于多峰功率谱模型建立了时变功率谱统计模型及其等效模型，在时域和频域方面均满足能量的非负性及边缘条件的性质，且精度较高。经过进一步的地震动特性分析和结构随机动力分析验证，认为基于欧进萍谱、Clough-Penzien谱及Kanai-Tajimi谱的地震动特性和结构响应与采用真实功率谱的结果差距明显，而时变功率谱统计模型及其等效模型均具有足够的精度，可为今后的地震动特性研究和结构动力分析提供更好的技术支持。