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

文章编号：1000-4750(2021)01-0027-13
基于Copula 理论的序列型地震动随机模型
申家旭1
，陈 隽1
，丁 国
2
(1. 同济大学土木工程学院，上海 200092；2. 华东建筑设计研究总院，上海 200002)
摘 要：通过考虑序列型地震动中主震和余震间的空间相关性，将已有地震动物理随机模型推广至序列型地震动随机模型。

依据序列型地震动的产生与传播机制，将其表示为包含16个随机变量的Fourier 模型。

考虑主震与余震的空间相关性，基于Copula 理论，将模型更准确地表示为6个二维随机变量和2个一维随机变量；为描述地震动在局部场地上传播的过程，建立了序列型地震动的随机场模型；根据1038组实测序列型地震动及5组地震动台阵数据，对模型中的物理随机变量进行参数识别和统计分析，结合Copula 函数，给出了各参数的概率分布。

与实测地震动对比表明：该文提出的序列型地震动随机模型，能够较为真实地再现序列型地震动的空间相关性以及地震动在局部场地的行波效应。

关键词：Copula 理论；序列型地震动；空间相关性；参数统计；联合分布
中图分类号：P315.9 文献标志码：A doi: 10.6052/j.issn.1000-4750.2020.03.0128
A STOCHASTIC MODEL FOR SEQUENTIAL GROUND MOTIONS
BASED ON THE COPULA THEORY
SHEN Jia-xu 1
, CHEN Jun 1
, DING Guo
2
(1. College of Civil Engineering, Tongji University, Shanghai 200092, China;2. East China Architectural Design & Research Institute, Shanghai 200002, China)
Abstract: The existing ground motion physical stochastic model is extended to a stochastic model for sequential ground motions by considering the spatial correlation between the mainshock and aftershock. Considering the physical mechanism of ground motion generation and propagation, the sequential ground motions are represented by the Fourier model of 16 physical random variables. Meanwhile, considering the spatial correlation between the mainshock and aftershock, the model is more accurately represented by six two-dimensional random variables and two one-dimensional random variables based on the Copula theory. To describe the process of the ground motion propagation on local sites, a random field model of sequential ground motions is established. Based on 1038 pairs of measured sequential ground motions and five sets of ground motion array data, the probability distribution of each parameter is given by parameter identification and statistical analyses of physical random variables in the model combined with the Copula function. The comparison between the measured and the simulated ground motions shows that the proposed stochastic model for sequential ground motions can realistically reproduce both the spatial correlation of sequential ground motions and the wave traveling effect of ground motions at local sites.Key words: Copula theory; sequential ground motion; spatial correlation; parameter statistics; joint distribution
历史地震资料表明，一次大地震发生后，存在再发生多次地震的可能性，后续发生的这一系列地震被称为余震，主震与余震共同组成序列型地震动，在其作用下，已损伤结构会发生二次损伤，甚至造成倒塌。

最典型的是2010年的新西兰基督城地震，7.0级的主震使城市中的小部分建筑
收稿日期：2020-03-02；修改日期：2020-06-19
基金项目：国家自然科学基金重点项目(U1711264)；新疆大学天山学者计划项目
通讯作者：陈　隽(1972−)，男，河南人，教授，博士，博导，从事工程结构抗震、土木工程大数据方面的研究(E-mail: ******************.cn ).作者简介：申家旭(1995−)，男，山东人，博士生，从事工程结构抗震研究(E-mail: ****************);
丁　国(1990−)，男，山东人，硕士，从事工程结构抗震研究(E-mail: ***************).
第 38 卷第 1 期Vol.38 No.1工 程 力 学2021年
1 月
Jan.
2021
ENGINEERING MECHANICS
27
遭受了严重破坏，而数月之后的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谱模型等)分别生成主震和余震，或通过地震动衰减关系挑选主震和余震。

此类方法实质上是主震随机模型，伴随的余震是由主余震之间的确定性关系而唯一确定的。

此外，很多主震随机模型本身基于平稳性假定，与地震动实际非平稳特性不符。

事实上，主震和余震在地震强度、频谱和持时等方面均密切相关，仅通过单一的地震动强度参数来反映二者的关系，显然不能准确全面地刻画出余震的特征。

此外，考虑地震动发生机理的多因素特征，余震同样具有相当的随机性。

因此，对主震和余震间相关性的研究，以及对地震动随机性和非平稳性的合理表征，是建立序列型地震动模型的关键难题。

聚焦于城市中某一区域(100 km2~101 km2)的抗震分析，覆盖面再广的观测台网中也仅会有少量的台站可以记录到由震源传递至此区域的地震动。

而受到地震波传播过程中的行波效应、场地效应等因素的影响，区域中任一点的地震动均存在显著的空间变异性[19]，仅靠台站得到的数条单点地震动难以准确表达此区域中每栋建筑所承受的地震动。

因此，对于城市区域建筑群的抗震研究，建立合适的序列型地震动场模型以反映城市区域中每一点处的地震动也是亟待解决的关键问题。

在一场大地震中，余震往往会发生非常多次，然而，多数余震震级较小，对结构安全影响较小。

在结构抗震分析中，学者们普遍关注主震与最大余震的组合，即“主震+最大余震”型序列型地震动。

本文的研究对象即为“主震+最大余震”型序列型地震动。

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

1 Copula理论
Copula理论最早由Sklar[20]提出，即任意一个多元联合分布都可以分解为相应的多个边缘分布和一个Copula函数，该Copula函数唯一地确定了变量间的相关性。

后经众多研究者的不断发展，Copula理论已经成为了解决高维随机变量联合概率分布问题的有效手段。

Copula理论的优势在于将变量边缘分布函数的估计与Copula函数的选择分开独立进行，并引入描述变量间非线性关系的秩相关系数，能够更加准确地描述变量的相关性。

Copula理论最早应用于金融、保险等领域，随后在水文学中得到广泛应用。

相对来说，其在地震工程中的应用还处于起步阶段。

仅有少量研究者将Copula理论引入地震动与结构抗震的研究中。

Goda和Salami[21]基于Copula函数建立了地震作用时结构最大变形和残余残余变形之间的关系；朱瑞广和吕大刚[22]利用Copula函数分析了主震和余震34个地震动强度参数的相关性，并基于Copula理论提出了余震的条件均值谱模型[17]。

Copula理论为确定两个或多个随机变量间的
28工 程 力 学
联合概率分布函数提供了一条有效且便捷的途径。

前已述及，Copula理论中引入了描述变量相关性的相关系数，其中应用较为广泛的为Pearson 线性相关系数，Kendall秩相关系数以及Spearman 秩相关系数。

本文选取能够反映变量间非线性关系的Kendall秩相关系数作为变量间相关性的度量。

由Kendall秩相关系数及参数的边缘分布，即可根据Copula函数建立变量间的联合概率分布函数。

以二维随机变量为例，其联合概率分布可表示为：
F(x1,x2)=C[F1(x1),F2(x2);θ]=C(u1,u2;θ)(1)式中：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相位谱模型，并将其中的经验参数改为了随机参数。

基于上述两项工作，场地一点处的地震动在时域范围内可以表示为：
a(t)=−
1
2π
+∞
−∞
A(λ,ω,R)·cos[ωt+
Φ(λ,ω,R)]dω(2)
A(λ,ω,R)
式中：R为震中距；ω为圆频率；为地震动加速度时程的Fourier幅值谱模型，其形式为：A(λ,ω,R)=
A0·ω·e−KωR
√
ω2+(1/τ)2
·
1+4ξ2g(ω/ωg)2
[1−(ω/ωg)2]2+4ξ2g(ω/ωg)2
(3)
Φ(λ,ω,R)为地震动加速度时程的Fourier相位谱模型，其形式为：
Φ(λ,ω,R)=arctan
(
1
τω
)
−R·ln[aω+1000b+
0.1323sin(3.78ω)+c cos(dω)](4)
λ=[A0,τ,a,b,c,d,ξg,ωg]
式中：为基本物理参数向量，A0为震源幅值系数，反映震源强度的影响，τ为Brune震源系数，反映震源断裂的时间过程特性，a、b、c和d为等效频率-波数关系的参数，ξg和ωg分别为局部场地的等效阻尼比和等效卓越圆频率。

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

从地震动的物理机制出发可以发现，尽管主震与余震在地震动参数以及反应谱中表现出一定的相关性和差异性，但其均遵循“震源-传播途径-局部场地”的物理过程。

而在这一物理过程中，其物理要素(即式(2)中的基本物理参数)具有不可完全观测和不可完全控制性，因而伴随主震而记录到的余震必然表现出显著的随机性。

主震与余
工 程 力 学29
震各物理要素在空间上的相关性决定了主震与余震时程的相关性。

因此，研究主震与余震各物理要素之间的相关性即可反映出主震与余震时程的相关性，这也是序列型地震动建模的关键。

考察式(2)，其中的λ为含有八个物理参数的随机向量，可以看到，地表一点处的地震动加速度时程由λ所决定。

在地震动产生、传播的过程中，震源空间位置的相关性导致了描述震源、传播途径的物理参数是存在相关性的。

地震动在局部场地的传播过程较为复杂，本文将局部场地的变化假定为线性，暂不考虑其非线性变化。

因此，描述局部场地动力特性的物理参数ξg和ωg在主震和余震中可视为是并未发生改变的，而主震与余震的空间相关性寓于描述震源和传播途径的六个物理参数中。

基于上述观点，地表一点处的序列型地震动加速度时程可表示为：
a(t)=−1
2π
+∞
−∞
A(Z,ω,R)·cos[ωt+
Φ(Z,ω,R)]dω(5)
Z=[λM,λA]
式中，为主震与余震基本物理参数向量，即：
Z=[A0M,τM,a M,b M,c M,d M,ξg M,ωg M,A0A,
τA,a A,b A,c A,d A,ξg A,ωg A](6)
Z
考虑其参数间的相关性，可改写为：
Z=[A0,τ,a,b,c,d,ξg,ωg](7) A0、τ、a、b、c、d
式中，为描述主余震相关关系的六个二维随机变量。

由于其不可观测性，参数值无法直接通过观测得到，可通过地震动反演由参数识别得到。

2.2 序列型地震动随机场物理模型
为描述地震动在局部场地上的空间相关性，Wang和Li[28]基于地震动随机函数模型(即地表一点处地震动模型)，考虑地震动场双尺度建模，建立了工程场地地震动随机场的物理模型。

缪惠全和李杰[29]在此基础上提出了工程场地地震动的相干函数模型，并与经典的Hao模型、Harichandran-Vanmarcke模型进行了对比，验证了其合理性。

场模型可表述为：
a(r l,t)=−1
2π
+∞
−∞
A(λ,ω,R)·exp
(
−
α0ωr l
2
)
·
cos [
ωt+Φ(λ,ω,R)−
ωr l
c g
]
dω(8)
式中：r l为场地任意一点沿波的传播方向上到场地中心的距离；α0为场地地震波衰减参数；c g为场地等效波速；a(r l, t)为场地内任意一点的加速度。

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

该模型考虑了地震动的行波效应和场地效应，准确简洁地反映了地震动产生、传播过程中震源、传播途径和局部场地物理机制的作用。

模型已应用于高层建筑[30]、超高层建筑[31]、大型冷却塔[32]等单体建筑的抗震可靠度分析和倒塌模拟，以及地下管网系统[33]等区域性生命线工程的抗震可靠性分析，有良好的工程应用基础。

第2.1节，本文提出了地表一点处的序列型地震动模型。

模型将同一地震事件中主震与余震视为两个随机过程，并进一步地考虑了两个随机过程中的随机变量在震源、传播途径中的相关性。

在考虑使用随机场描述场地内各点序列型地震动的空间相关性时，由于余震震源位置存在不确定性，本文将地震动东-西分量与南-北分量视为相互独立的地震动时程记录，使用此方法处理，可将地表一点处主震与余震分量以相同方向表示。

而由于本随机场刻画的是局部场地范围内地震动的行波效应和场地效应，因此，随机场的两个物理参数α0和c g在主震和余震中视为相同的物理量。

基于此，在工程场地范围内，序列型地震动随机场物理模型可表示为：
a(r l,t)=−
1
2π
+∞
−∞
A(Z,ω,R)·exp
(
−
α0ωr l
2
)
·
cos
[
ωt+Φ(Z,ω,R)−
ωr l
c g
]
dω(9) 3 序列型地震动的参数识别和统计
为确定式(9)中各物理参数的概率分布及主余震参数的联合概率分布，本节收集并筛选了不同场地类型的实测序列型地震动，并进行重采样及归一化处理；随后对主震和余震各物理参数分别统计，识别出各参数的边缘概率分布；最后根据各参数的相关系数及边缘概率分布，识别出最优Copula函数，给出式(9)中各二维参数的联合概率分布。

3.1 实测序列型地震动的收集和处理
本文从太平洋地震工程研究中心[34](Pacific Earthquake Engineering Research Center, PEER)和KiK-net & K-NET中收集得到了大量的序列型地震
30工 程 力 学
动记录。

通过如下准则对其进行筛选：
1) 选取一次地震动事件中主震和伴随其发生的一系列余震中震级最大的一次余震；
2) 选取同一台站记录的主震和余震并组合；
3) 为避免土-结构的耦合作用，台站位于自由场地或较小建筑物的地面层；
4) 主震和余震的震级分别不小于5级，排除不太可能对结构造成影响的地震；
5) 为了研究场地类别对结构响应及损伤状态的影响，将序列型地震动按场地类别进行分类，场地类别依据中国抗震规范划分为四类场地，按照20 m土层厚等效剪切波速V S20和30 m厚等效剪切波速V S30划分，划分标准如表1。

由于KiK-net & K-NET按照分层土记录土层信息，为便于分类，统一按照中国规范计算其V S20，计算依据参考郭锋[35]提出的V S20与V S30转换关系式。

表 1 场地划分标准及序列型地震动数量
Table 1 Site classification standard and number of
sequential ground motions
场地条件ⅠⅡⅢⅣ
等效剪切波速V S20/(m/s)>500250~500150~250<150
等效剪切波速V S30/(m/s)>596278~596158~278<158地震动数量/组130********
共筛选得到1038组符合条件的序列型地震动，另外，选取了SMART1台阵测得的五组地震动场记录以作地震动场参数的识别，其基本信息如表2所示。

本文主要研究地震动水平分量，故选择各台站记录中的E-W分量和N-S分量，将其视为相互独立的地震动处理。

表 2 筛选得到的实测序列型地震动
Table 2 Selection of measured sequential ground motions
地震事件
主震余震
发生时间Mw数量发生时间Mw数量L′Aquila2009.04.06 01:336.3642009.04.06 17:475.664 Livermore1980.01.24 19:005.8121980.01.27 02:335.412 Mammoth Lakes1980.05.25 16:346.161980.05.25 16:495.76 Chi-Chi1999.09.20 17:477.64801999.09.20 17:575.9480 Chalfant Valley1986.07.20 14:295.8101986.07.21 14:426.210 Northridge1994.01.17 12:316.7621994.03.20 21:205.362 Umbria Marche1997.09.26 09:406.0181997.10.06 23:245.518 East Japan2011.03.11 14:469.0242011.03.11 15:157.724
Alaska2002.10.23 03:306.7382002.11.03 22:127.938 Imperial Valley1979.10.15 23:166.5281979.10.15 23:195.028 New Zealand2010.09.03 16:357.01542011.02.21 23:516.2154 Whittier Narrows1987.10.01 14:426.01341987.10.04 10:595.3134 Kumamoto2016.04.14 21:266.282016.04.16 01:257.08
SMART1 Events
1981.01.29 04:515.954−−−
1983.09.21 19:206.570−−−
1985.06.12 17:225.866−−−
1986.05.20 05:256.372−−−
1986.11.14 21:207.372−−−各台站测得的地震动时程采样频率存在差异，为便于分析，统一将地震动以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代入
工 程 力 学31
余震幅值谱模型中，识别余震模型中的参数A 0A 和τA ；
4) 利用真实地震动的相位谱，求得其相位差谱分布；
5) 将2)和3)中的震源系数τM 和τA 代入式(3)的相位谱模型中，采用最小二乘准则将模型的相位差谱与真实相位差谱分布拟合，利用遗传算法识别分别得到主震和余震相位谱模型中的参数a M 、b M 、c M 、d M 和a A 、b A 、c A 、d A 。

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

在各参数中，幅值谱对应的A 0M 、τM 、ξgM 和ωgM 以使模型逼近真实地震动幅值谱为原则，识别方法采用最小二乘法，可表述为：
J =min
∑ A (Z ,ω,R )−˜A (ω) 2
(10)
A (λ,ω,R )˜A
(ω)式中：为本文提出的Fourier 幅值谱模
型；为真实地震动幅值谱值。

与幅值谱不同，相位谱往往类似于白噪声，难以根据真实地震动的相位谱识别其参数。

然而，地震动相位差谱却有较为明显的特征，相位差一般符合正态分布或近似正态分布。

基于此，
相位谱对应的a 、b 、c 、d 由以使模型的相位差分布逼近真实地震动相位差分布为原则，采用遗传算法作为识别方法，目标函数可表述为：
J =min ∑
|g (∆ϕk )−˜g
(∆ϕk )|2(11)
g (∆ϕk )˜g
(∆ϕk )∆ϕk =ϕk +1−ϕk 式中：为本文提出的相位谱模型的相位差
分布；为真实地震动相位谱的相位差分布；为频域上相邻两个相位角的差。

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

表 3 遗传算法参数设置
Table 3 Settings of Genetic Algorithm
参数类型参数设置值创建初始种群方法
二进制编码法
种群规模200
适应度函数排列尺度变换函数选择操作函数随机均匀分布选择法变异函数均匀变异函数
变异概率0.04
交叉函数离散重组交叉方式
交叉概率0.8
终止条件最大进化代数为500代
函数精度
10
−8
对于α0和c g 的识别，首先识别出场地中心点C00的各物理参数A 0、τ、ξg 、ωg 、a 、b 、c 、
d ，由最小二乘法确定幅值谱中的参数α0。

对于场地等效波速c g ，可通过下式识别：
c g =
∆r l ∆t
(12)
式中：Δr l 为两样本点在平面波场中的坐标差；Δt 为两样本的时间延迟，可采用Boissières 和Vanmarcke [36]
提出的方法确定。

根据上述过程，对各组序列型地震动逐一识别，得到描述主震和余震的各物理参数在四类不同场地类型下的分布规律。

考察其统计分布，采用贝叶斯信息准则(Bayesian Information Criterion ，BIC) 准则
[37]
求得最优概率模型，获得主震和余震
各参数的概率密度函数，对应的概率密度函数统计值如表4和表5所示。

表中参数P 1、参数P 2的含义：当概率密度函数为Normal 时，分别为均值和标准差；概率密度函数为Lognormal 时，为对数的均值和标准差；概率密度函数为Weibull 时，为形状参数和尺度参数。

需要指出的是，在对
表 4 概率密度函数及其对应参数
Table 4 Probability density functions and parameters
随机变量场地
主震
余震函数类型
P 1
P 2
函数类型
P 1
P 2
A 0
Ⅰ
Lognormal
−3.043 1.270
Lognormal
−4.660 1.082Ⅱ−2.852 1.259−4.537 1.300Ⅲ−1.797 0.936−3.427 1.005Ⅳ−1.606 0.805−2.196 0.997τ
ⅠWeibull
0.719
0.685Weibull 0.6690.479Ⅱ
0.6610.3590.6340.380Ⅲ0.7260.3890.8480.392Ⅳ0.8110.4830.7380.160a
ⅠLognormal
1.886
0.550Lognormal 2.1500.462Ⅱ
1.8120.659
2.0540.544Ⅲ 1.4610.721 1.8340.643Ⅳ 1.0460.739 1.7740.638b
ⅠLognormal
2.012
0.498Lognormal 2.3120.458Ⅱ
1.9330.557
2.2200.528Ⅲ 1.7460.631 1.9400.587Ⅳ 1.4210.514 1.7900.637c
ⅠWeibull
1.272
1.231Weibull 1.512 1.582Ⅱ
1.315 1.271 1.659 1.700Ⅲ 1.443 1.351 1.587 1.661Ⅳ 1.674 1.243 1.749 1.858d
Ⅰ
Weibull
1.682
1.591Weibull 1.729 1.782Ⅱ 1.565 1.615 1.844 1.788Ⅲ 1.625 1.607 1.895 1.835Ⅳ
1.631
1.684
1.896
2.016
32工 程 力 学
α0和c g进行识别前，应按照场地条件对地震动记录进行分组以满足工程需要。

然而，用于观测局部场地的地震动台阵较少，易于获取的地震动场记录均来自于SMART1台阵，其对应的场地为Ⅱ类场地。

因此，本文给出的α0和c g的统计结果严格来说仅适用于Ⅱ类场地。

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

其中BIC定义为：
BIC=−2
N
∑
i=1
ln f(x i;p,q,r)+k1ln N(13)
式中：f (x i; p, q, r)为概率分布函数的密度函数，x i (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函数；
表 5 概率密度函数及其对应参数
Table 5 Probability density functions and parameters 随机变量场地函数类型P1P2
ξg Ⅰ
Weibull
0.524 2.688Ⅱ 0.484 2.709Ⅲ 0.494 2.595Ⅳ 0.356 1.928
ωg Ⅰ
Weibull
11.261 1.130Ⅱ12.043 1.224Ⅲ 5.6100.851Ⅳ11.409 2.827
α0ⅡNormal 1.014×10−4 2.214×10−3
c gⅡLognormal 1.4140.278
工 程 力 学33
4) 由最优Copula函数求得主余震参数的联合概率分布函数。

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

通过参数间的联合PDF，本文根据等概率变换原则模拟出了符合对应联合PDF的散点图。

同样以A0为例，图5给出了4类不同场地条件下，参数的真实数据与模拟数据的散点图。

可见，模拟数据和真实数据的分布保持着较高的吻合度。

因此，可以认为，本文利用Copula函数得到的参数间的联合PDF能够准确描述序列型地震动模型
图 2 A0M和A0A的边缘概率密度函数
Fig. 2 Marginal PDFs of A0M and A0A
图 3 A0M和A0A的散点图
Fig. 3 Scatters of A0M and A0A 34工 程 力 学
中各二维随机变量的取值和分布，即能够较好地描述主震和余震参数的相关性。

图 5 实测数据与模拟数据对比
Fig. 5 Comparisons of real data and simulation
4 序列型地震动的验证
为验证序列型地震动模型的有效性，对实测
表 6 最优Copula 函数及相关参数
Table 6 Optimal Copula functions and relevant parameters
随机变量场地函数类型相关参数θA 0
Ⅰ
Plackett
9.896
Ⅱ18.808Ⅲ26.833Ⅳ 7.156τ
ⅠClayton
2.214
Ⅱ
2.601Ⅲ 1.287Ⅳ 1.951a
ⅠPlackett
2.710
Ⅱ
3.592Ⅲ 2.840Ⅳ 3.057b
ⅠClayton
0.916
Ⅱ
1.230Ⅲ 1.262Ⅳ 1.442c
ⅠClayton
0.467
Ⅱ
0.493Ⅲ 0.407Ⅳ 1.400
d
Ⅰ
t (ν=4)
0.211
Ⅱ 0.255Ⅲ 0.271Ⅳ
0.659
图 4 识别最优Copula 函数的流程图
Fig. 4 Flowchart for identifying optimal Copula function
工 程 力 学35。