基于Copula 理论的序列型地震动随机模型

申家旭，陈 隽，丁 国

(1. 同济大学土木工程学院，上海 200092；2. 华东建筑设计研究总院，上海 200002)

历史地震资料表明，一次大地震发生后，存在再发生多次地震的可能性，后续发生的这一系列地震被称为余震，主震与余震共同组成序列型地震动，在其作用下，已损伤结构会发生二次损伤，甚至造成倒塌。最典型的是2010 年的新西兰基督城地震，7.0 级的主震使城市中的小部分建筑遭受了严重破坏，而数月之后的6.3 级的强余震使已损伤的建筑进一步受损，导致了大量建筑破坏和倒塌，造成了比主震严重的多的经济损失和人员伤亡[1]。在我国1999 年集集地震以及2008 年汶川地震中，序列型地震动不仅造成了大量的结构破坏，而且加大了灾后救援的难度，使人民的生命财产受到了严重损失。研究序列型地震动的特性和生成模型，成为工程结构抗震研究中一个重要的现实问题。

地震动的合理数学模型是工程结构抗震分析的基础。1947 年，Housner[2]指出工程结构在地震作用下的响应主要受两个因素的影响：一是结构自身的动力特性；二是地震动的特性。根据当时为数不多的地震动记录，Housner[2]认为实际地震动具有很强的随机性，并提出了最早的随机地震动模型——随机独立脉冲模型。近几十年来，世界各国陆续建立了覆盖面广的强震观测台网，极大地丰富了序列型地震动的实测数据，为人们认识和研究地震动提供了良好的数据基础。基于大量的观测记录，不同的地震动模型相继提出，其中具有代表性的有Kanai-Tajimi 模型[3−4]、谱表现方法[5−6]、演变功率谱模型[7]等。这些模型从不同角度表现了地震动的特性，并得到了广泛的应用。对于序列型地震动的模拟，早期的研究主要以重复主震或随机组合主余震记录[8−10]等方式来实现。然而，多位研究者已指出此类方法会过高地估计结构响应[11−13]。现行方法则主要集中在建立主震与余震地震动强度参数[14−18](如震级、PGA等)之间的统计学关系，通过随机地震动模型(如Kanai-Tajimi 谱模型等)分别生成主震和余震，或通过地震动衰减关系挑选主震和余震。此类方法实质上是主震随机模型，伴随的余震是由主余震之间的确定性关系而唯一确定的。此外，很多主震随机模型本身基于平稳性假定，与地震动实际非平稳特性不符。事实上，主震和余震在地震强度、频谱和持时等方面均密切相关，仅通过单一的地震动强度参数来反映二者的关系，显然不能准确全面地刻画出余震的特征。此外，考虑地震动发生机理的多因素特征，余震同样具有相当的随机性。因此，对主震和余震间相关性的研究，以及对地震动随机性和非平稳性的合理表征，是建立序列型地震动模型的关键难题。

聚焦于城市中某一区域(100km2～101km2)的抗震分析，覆盖面再广的观测台网中也仅会有少量的台站可以记录到由震源传递至此区域的地震动。而受到地震波传播过程中的行波效应、场地效应等因素的影响，区域中任一点的地震动均存在显著的空间变异性[19]，仅靠台站得到的数条单点地震动难以准确表达此区域中每栋建筑所承受的地震动。因此，对于城市区域建筑群的抗震研究，建立合适的序列型地震动场模型以反映城市区域中每一点处的地震动也是亟待解决的关键问题。

在一场大地震中，余震往往会发生非常多次，然而，多数余震震级较小，对结构安全影响较小。在结构抗震分析中，学者们普遍关注主震与最大余震的组合，即“主震+最大余震”型序列型地震动。本文的研究对象即为“主震+最大余震”型序列型地震动。

综上，本文从地震动产生的物理机制出发，在已有的工程场地地震动模型的基础上，考虑主震与余震在产生、传播机制上的相关性，利用Copula 理论研究主震和余震各参数相关性，进而考虑局部场地上不同点的空间相关性，形成序列型地震动随机场模型。

1 Copula 理论

Copula 理论最早由Sklar[20]提出，即任意一个多元联合分布都可以分解为相应的多个边缘分布和一个Copula 函数，该Copula 函数唯一地确定了变量间的相关性。后经众多研究者的不断发展，Copula 理论已经成为了解决高维随机变量联合概率分布问题的有效手段。Copula 理论的优势在于将变量边缘分布函数的估计与Copula 函数的选择分开独立进行，并引入描述变量间非线性关系的秩相关系数，能够更加准确地描述变量的相关性。Copula 理论最早应用于金融、保险等领域，随后在水文学中得到广泛应用。相对来说，其在地震工程中的应用还处于起步阶段。仅有少量研究者将Copula 理论引入地震动与结构抗震的研究中。Goda 和Salami[21]基于Copula 函数建立了地震作用时结构最大变形和残余残余变形之间的关系；朱瑞广和吕大刚[22]利用Copula 函数分析了主震和余震34 个地震动强度参数的相关性，并基于Copula 理论提出了余震的条件均值谱模型[17]。

Copula 理论为确定两个或多个随机变量间的联合概率分布函数提供了一条有效且便捷的途径。前已述及，Copula 理论中引入了描述变量相关性的相关系数，其中应用较为广泛的为Pearson线性相关系数，Kendall 秩相关系数以及Spearman秩相关系数。本文选取能够反映变量间非线性关系的Kendall 秩相关系数作为变量间相关性的度量。由Kendall 秩相关系数及参数的边缘分布，即可根据Copula 函数建立变量间的联合概率分布函数。以二维随机变量为例，其联合概率分布可表示为：

式中：F1(x1)、F2(x2)为变量x1和x2的边缘累积分布函数；u1=F1(x1)，u2=F2(x2)，均为[0,1]间的均匀分布；C(u1,u2)为Copula 分布函数，其变量区间为[0,1]；θ 为与Copula 分布函数对应的相关参数，可由对应的Copula 分布函数及Kendall 秩相关系数求出。

为确定序列型地震动中主震与余震的相关性，本文采用Copula 理论定量分析并确定了主余震物理参数间的相关关系，并给出其显式表达。

2 序列型地震动随机场物理模型

主震与余震地震动参数之间的关系一直是地震学领域的研究重点。在地震动研究的早期，描述主震与余震关系的地震学三大定律[23−25]相继提出。随后的数十年里，大量的研究者们对其进行了修正和改进，力求真实地反映主余震的关系。而在地震工程学领域，研究者们更多地关注影响结构安全性的参数(如PGA、PGV、Sa 等)，着眼于将地震动参数与结构时程分析或谱分析相结合，研究主余震之间卓越周期、谱加速度等参数的相关性，以指导工程实践。如前所述，现有方法一般通过主余震地震动参数的统计关系来建立序列型地震动模型。由于统计关系本身的经验性特征，不能准确地体现主余震相关性的物理机制，且其本质是主震与余震一一对应的确定性模型，这与主震过后往往出现大量随机性余震的事实不符。更为关键的是，大量的地震事件表明：主震与余震在震源及传播途径上并不完全相同，主震断层的二次破裂及临近断层的破裂均会产生余震，因此，将主震与余震视为完全相关并不能准确表达出余震的自身特性。

2.1 序列型地震动物理随机函数模型

为克服描述单次地震的现象学随机地震动模型(如Kanai-Tajimi 谱模型，演变功率谱模型等)的经验统计的局限性，Wang 和Li[26]从地震动“震源-传播途径-局部场地”的物理机制出发，在随机Fourier 谱和随机函数模型的基础上分别对震源、传播途径、局部场地进行建模，给出了随机地震动的Fourier 幅值谱和相位谱模型。模型中影响地震动随机性的关键物理因素被抽象为随机变量，由此得到了地表一点处的地震动模型。宋萌[27]在此基础上优化了Fourier 相位谱模型，并将其中的经验参数改为了随机参数。基于上述两项工作，场地一点处的地震动在时域范围内可以表示为：

限于篇幅，此模型更详细的理论推导可参考文献 [26 − 27]。

从地震动的物理机制出发可以发现，尽管主震与余震在地震动参数以及反应谱中表现出一定的相关性和差异性，但其均遵循“震源-传播途径-局部场地”的物理过程。而在这一物理过程中，其物理要素(即式(2)中的基本物理参数)具有不可完全观测和不可完全控制性，因而伴随主震而记录到的余震必然表现出显著的随机性。主震与余震各物理要素在空间上的相关性决定了主震与余震时程的相关性。因此，研究主震与余震各物理要素之间的相关性即可反映出主震与余震时程的相关性，这也是序列型地震动建模的关键。

考察式(2)，其中的λ 为含有八个物理参数的随机向量，可以看到，地表一点处的地震动加速度时程由λ 所决定。在地震动产生、传播的过程中，震源空间位置的相关性导致了描述震源、传播途径的物理参数是存在相关性的。地震动在局部场地的传播过程较为复杂，本文将局部场地的变化假定为线性，暂不考虑其非线性变化。因此，描述局部场地动力特性的物理参数ξg和ωg在主震和余震中可视为是并未发生改变的，而主震与余震的空间相关性寓于描述震源和传播途径的六个物理参数中。基于上述观点，地表一点处的序列型地震动加速度时程可表示为：

2.2 序列型地震动随机场物理模型

为描述地震动在局部场地上的空间相关性，Wang 和Li[28]基于地震动随机函数模型(即地表一点处地震动模型)，考虑地震动场双尺度建模，建立了工程场地地震动随机场的物理模型。缪惠全和李杰[29]在此基础上提出了工程场地地震动的相干函数模型，并与经典的Hao 模型、Harichandran-Vanmarcke 模型进行了对比，验证了其合理性。场模型可表述为：

式中：rl为场地任意一点沿波的传播方向上到场地中心的距离；α0为场地地震波衰减参数；cg为场地等效波速；a(rl, t)为场地内任意一点的加速度。

限于篇幅，随机场模型更详细的理论推导可参考文献 [28 − 29]。

该模型考虑了地震动的行波效应和场地效应，准确简洁地反映了地震动产生、传播过程中震源、传播途径和局部场地物理机制的作用。模型已应用于高层建筑[30]、超高层建筑[31]、大型冷却塔[32]等单体建筑的抗震可靠度分析和倒塌模拟，以及地下管网系统[33]等区域性生命线工程的抗震可靠性分析，有良好的工程应用基础。

第2.1 节，本文提出了地表一点处的序列型地震动模型。模型将同一地震事件中主震与余震视为两个随机过程，并进一步地考虑了两个随机过程中的随机变量在震源、传播途径中的相关性。在考虑使用随机场描述场地内各点序列型地震动的空间相关性时，由于余震震源位置存在不确定性，本文将地震动东-西分量与南-北分量视为相互独立的地震动时程记录，使用此方法处理，可将地表一点处主震与余震分量以相同方向表示。而由于本随机场刻画的是局部场地范围内地震动的行波效应和场地效应，因此，随机场的两个物理参数α0和cg在主震和余震中视为相同的物理量。基于此，在工程场地范围内，序列型地震动随机场物理模型可表示为：

3 序列型地震动的参数识别和统计

为确定式(9)中各物理参数的概率分布及主余震参数的联合概率分布，本节收集并筛选了不同场地类型的实测序列型地震动，并进行重采样及归一化处理；随后对主震和余震各物理参数分别统计，识别出各参数的边缘概率分布；最后根据各参数的相关系数及边缘概率分布，识别出最优Copula 函数，给出式(9)中各二维参数的联合概率分布。

3.1 实测序列型地震动的收集和处理

本文从太平洋地震工程研究中心[34](Pacific Earthquake Engineering Research Center, PEER)和KiK-net & K-NET 中收集得到了大量的序列型地震动记录。通过如下准则对其进行筛选：

1) 选取一次地震动事件中主震和伴随其发生的一系列余震中震级最大的一次余震；

2) 选取同一台站记录的主震和余震并组合；

3) 为避免土-结构的耦合作用，台站位于自由场地或较小建筑物的地面层；

4) 主震和余震的震级分别不小于5 级，排除不太可能对结构造成影响的地震；

5) 为了研究场地类别对结构响应及损伤状态的影响，将序列型地震动按场地类别进行分类，场地类别依据中国抗震规范划分为四类场地，按照20 m 土层厚等效剪切波速VS20和30 m 厚等效剪切波速VS30划分，划分标准如表1。由于KiKnet & K-NET 按照分层土记录土层信息，为便于分类，统一按照中国规范计算其VS20，计算依据参考郭锋[35]提出的VS20与VS30转换关系式。

表 1 场地划分标准及序列型地震动数量Table 1 Site classification standard and number of sequential ground motions

共筛选得到1038 组符合条件的序列型地震动，另外，选取了SMART1 台阵测得的五组地震动场记录以作地震动场参数的识别，其基本信息如表2 所示。本文主要研究地震动水平分量，故选择各台站记录中的E-W 分量和N-S 分量，将其视为相互独立的地震动处理。

表 2 筛选得到的实测序列型地震动Table 2 Selection of measured sequential ground motions

各台站测得的地震动时程采样频率存在差异，为便于分析，统一将地震动以0.02 s 的时间间隔进行重采样，频率限为25 Hz。重采样结果与原始地震动时程在时域和频域的差异均极其微小，对识别结果的影响可以忽略。另外，地震动的峰值统一调幅至0.1 g，作幅值归一化处理。

3.2 物理随机参数的识别与统计

利用实测序列型地震记录，针对第2.1 节和第2.2 节的模型参数的识别过程如图1 所示。

图 1 参数识别流程图Fig. 1 Flowchart of parameter identification

1) 对地震动时程作Fourier 变换，得到其Fourier幅值谱和Fourier 相位谱；

2) 由式(2)的幅值谱模型，采用最小二乘准则拟合，识别主震幅值谱模型中的物理参数A0M、τM、 ξgM和 ωgM；

3) 将 2)中的局部场地参数 ξgM和 ωgM代入余震幅值谱模型中，识别余震模型中的参数A0A和 τA；

4) 利用真实地震动的相位谱，求得其相位差谱分布；

5) 将2)和3)中的震源系数τM和τA代入式(3)的相位谱模型中，采用最小二乘准则将模型的相位差谱与真实相位差谱分布拟合，利用遗传算法识别分别得到主震和余震相位谱模型中的参数aM、bM、cM、dM和 aA、bA、cA、dA。

6) 对地震动场参数作识别，得到α0和cg。

在各参数中，幅值谱对应的A0M、τM、ξgM和ωgM以使模型逼近真实地震动幅值谱为原则，识别方法采用最小二乘法，可表述为：

遗传算法的参数设置，如表3 所示。

表 3 遗传算法参数设置Table 3 Settings of Genetic Algorithm

对于α0和cg的识别，首先识别出场地中心点 C00 的各物理参数 A0、τ、ξg、ωg、a、b、c、d，由最小二乘法确定幅值谱中的参数α0。对于场地等效波速cg，可通过下式识别：

式中：Δrl为两样本点在平面波场中的坐标差；Δt 为两样本的时间延迟，可采用Boissières 和Vanmarcke[36]提出的方法确定。

根据上述过程，对各组序列型地震动逐一识别，得到描述主震和余震的各物理参数在四类不同场地类型下的分布规律。考察其统计分布，采用贝叶斯信息准则(Bayesian Information Criterion，BIC) 准则[37]求得最优概率模型，获得主震和余震各参数的概率密度函数，对应的概率密度函数统计值如表4 和表5 所示。表中参数P1、参数P2的含义：当概率密度函数为Normal 时，分别为均值和标准差；概率密度函数为Lognormal 时，为对数的均值和标准差；概率密度函数为Weibull 时，为形状参数和尺度参数。需要指出的是，在对α0和cg进行识别前，应按照场地条件对地震动记录进行分组以满足工程需要。然而，用于观测局部场地的地震动台阵较少，易于获取的地震动场记录均来自于SMART1 台阵，其对应的场地为Ⅱ类场地。因此，本文给出的α0和cg的统计结果严格来说仅适用于Ⅱ类场地。

表 4 概率密度函数及其对应参数Table 4 Probability density functions and parameters

表 5 概率密度函数及其对应参数Table 5 Probability density functions and parameters

以参数A0为例，图2 给出了主震A0M和余震A0A的统计分布及其概率密度函数(Probability Density Function，PDF)曲线。其中BIC 定义为：

式中：f (xi; p, q, r)为概率分布函数的密度函数，xi(i=1, 2,···, N)为样本点数目，p, q, r 为分布参数；k1为分布参数的数目。

3.3 主余震的相关性分析

以参数A0为例，图3 给出了各类场地下主震A0M与余震A0A的关系，其中IV 类场地下收集到的地震动数量较少。分别计算各参数的Kendall 秩相关系数，A0、τ、a、b、c、d 的相关系数分别为：0.589、0.508、0.277、0.395、0.195、0.172。可知，主震与余震的参数间存在一定的相关性，本文通过引入Copula 理论对其进行描述，并给出主余震参数的联合概率分布，实现对主余震参数相关性的显式表达。

为此，根据第3.2 节中得到的各参数的边缘分布函数，利用二维Copula 函数建立主震和余震参数之间的联合分布。共选取了七种常用的Copula函数：Independent、Gaussian、t、Gumbel、Plackett、Frank、和Clayton 函数，其中，考虑到参数c 和d 的相关性较弱的情况，在备选Copula 中加入了Independent Copula 函数，代表了参数相互独立的情况。其余六种函数分别从不同角度描述了参数间的相关性。

图4 给出了确定两参数间联合概率分布函数的步骤：

1) 计算主震与余震参数的Kendall 秩相关系数τ；

2) 由主震和余震参数的边缘概率分布f1、f2和秩相关系数τ，求出备选Copula 函数的相关参数θ；

3) 利用BIC 作判别，遍历七种常用Copula 函数，求得最优Copula 函数；

图 2 A0M 和A0A 的边缘概率密度函数Fig. 2 Marginal PDFs of A0M and A0A

4) 由最优Copula 函数求得主余震参数的联合概率分布函数。

识别得到的各参数最优Copula 函数及相关参数，如表6 所示。

通过参数间的联合PDF，本文根据等概率变换原则模拟出了符合对应联合PDF 的散点图。同样以A0为例，图5 给出了4 类不同场地条件下，参数的真实数据与模拟数据的散点图。可见，模拟数据和真实数据的分布保持着较高的吻合度。因此，可以认为，本文利用Copula 函数得到的参数间的联合PDF 能够准确描述序列型地震动模型中各二维随机变量的取值和分布，即能够较好地描述主震和余震参数的相关性。

图 3 A0M 和 A0A 的散点图Fig. 3 Scatters of A0M and A0A

图 4 识别最优Copula 函数的流程图Fig. 4 Flowchart for identifying optimal Copula function

表 6 最优Copula 函数及相关参数Table 6 Optimal Copula functions and relevant parameters

图 5 实测数据与模拟数据对比Fig. 5 Comparisons of real data and simulation

4 序列型地震动的验证

为验证序列型地震动模型的有效性，对实测的主余震地震动加速度时程进行模拟重构，并与实测地震动的均值反应谱进行对比验证。另外，选取了2018 年阿拉斯加地震记录到的地震动时程，对其进行模拟重构，进一步验证本文模型的合理性。

首先，对每条实测的地震动时程样本，通过式(5)识别出模型的各参数，再进行重构，最后得出模拟地震动。图6 对比了各类场地条件下主震和余震的均值加速度反应谱，阻尼比ξ=0.05。模型与实测地震动的均值加速度反应谱基本吻合，可见式(5)模型与表4～表6 的统计结果具有较高的可信度。

图 6 实测地震动与模拟地震动的均值反应谱对比Fig. 6 Comparisons of mean spectra of measured and simulated records

2018 年阿拉斯加地震发生于 2018 年 11 月30 日，其中，主震震级为Mw 7.0，余震震级为Mw 5.7。本文选取了AKK217 和AKK220 测站测得的序列型地震动，分别识别出主震和余震参数，依据式(5)进行重构，结果如图7 所示。对比可见，序列型地震动模型从样本时程的角度能和实测时程保持较高的一致性，可以较好地重构地震动时程。需要指出，模拟时程并不能很好地描述地震动的衰减阶段，衰减段持续时间较短，这是模型需要改进之处。

图 7 实测单点时程与模拟时程的对比Fig. 7 Comparisons of measured and simulated records

此外，本文通过对SMART-1 台阵的实测时程作模拟，验证地震动场模型的合理性，选取了Event 5 地震动事件中C00 测站和I06 测站的实测加速度时程，记录的加速度时程如图8(a)所示。地震动在局部场地传播过程中，地震波传播速度是有限的，因此不同测站记录到的地震动存在着相位差，在时程上表现为时间延迟现象。

本文首先通过式(9)识别得到Event 5 地震的物理参数，并模拟出SMART-1 局部场地原点即C00 的地震动时程。进而利用识别得到的物理参数，以及表5 中的α0和cg，模拟出沿传播方向与C00 的Δrl为190 m 的I06 测站时程。需要说明，模拟时程中，除α0和cg外，其余物理参数都相同。模拟结果如图8(b)所示。对比模拟得到的时程可见，其呈现出与实测时程相似的时间延迟现象。

图 8 实测台阵时程与场模型模拟时程对比Fig. 8 Comparisons of measured and simulated records

5 结论

通过考虑主震和余震的空间相关性，建立了序列型地震动随机单点模型，为反映地震动在局部场地的变异性，进一步地提出了序列型地震动随机场模型。利用1038 组序列型地震动的数据，识别得到了模型中各物理参数的统计分布。运用Copula 理论给出了主震和余震空间相关性的显式表达。通过将模拟地震动与实测单点地震动、地震动场作对比，验证了本文提出的模型的可行性。

模型将主震和余震的关系以概率分布的形式给出，并且在已知主震的情况下可以根据Copula理论得到余震参数的条件概率分布，保留了余震自身的随机性。

在实际应用中，确定工程场地条件后，从表4～表6 中确定各参数的边缘概率分布和Copula 函数确定参数的联合概率分布，据此生成各参数的随机样本，代入式(5)或式(9)中，即可模拟出序列型地震动，应用于单体建筑或区域性建筑群的抗震分析。

地震动的产生与传播过程极其复杂，如何准确描述地震的震源及传播途径一直是地震动建模的重点和难点。地震产生的机理复杂，传播途径中的影响地震动传播的因素较多，本文提出的模型在描述地震震源及传播途径方面做了较多简化；模型模拟出的时程在衰减段与实测时程有一定差异；另外，地震动在局部场地上的传播仍有其他影响因素，本文在此也作了一定简化，仅考虑了其行波效应。这些不足之处仍有待进一步研究。