圆环对撞问题的理想刚塑性小变形分析及讨论

张铁光; 梁成; 梁卫民

摘要: 本文对对经受到撞击的圆环的刚塑性动力响应进行了小变形分析,计算结果与大变形分析作了对比。指出二者之间的差别主要是由于动力学方程不同引起的。

关键词:

Abstract: This paper presents an analysis of the small rigid-plastic deformation of an unsupported thin ring subjected to two equal and opposite impacts along a diameter by two rigid masses. Numerical examples are shown in order to make a comparison with the results on the basis of large deformation formulation^[4] of the same problem, It is found that the set of dynamic equations of small deformation solution is different from that of large deformation solution^[4].

Keywords:

历史地震资料表明，一次大地震发生后，存在再发生多次地震的可能性，后续发生的这一系列地震被称为余震，主震与余震共同组成序列型地震动，在其作用下，已损伤结构会发生二次损伤，甚至造成倒塌。最典型的是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谱模型等)分别生成主震和余震，或通过地震动衰减关系挑选主震和余震。此类方法实质上是主震随机模型，伴随的余震是由主余震之间的确定性关系而唯一确定的。此外，很多主震随机模型本身基于平稳性假定，与地震动实际非平稳特性不符。事实上，主震和余震在地震强度、频谱和持时等方面均密切相关，仅通过单一的地震动强度参数来反映二者的关系，显然不能准确全面地刻画出余震的特征。此外，考虑地震动发生机理的多因素特征，余震同样具有相当的随机性。因此，对主震和余震间相关性的研究，以及对地震动随机性和非平稳性的合理表征，是建立序列型地震动模型的关键难题。

聚焦于城市中某一区域(10⁰ km²~10¹ km²)的抗震分析，覆盖面再广的观测台网中也仅会有少量的台站可以记录到由震源传递至此区域的地震动。而受到地震波传播过程中的行波效应、场地效应等因素的影响，区域中任一点的地震动均存在显著的空间变异性^[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函数建立变量间的联合概率分布函数。以二维随机变量为例，其联合概率分布可表示为：

$F({x_1},{x_2})=C[{F_1}({x_1}),{F_2}({x_2});\theta ]=C({u_1},{u_2};\theta)$

(1)

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

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

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

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

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

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

$\begin{split} a(t)= &- \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {A({\\text{λ}},\omega ,R)} \cdot \cos [\omega t +\\& {\varPhi }({\\text{λ}},\omega ,R)]{\rm{d}}\omega \end{split}$

(2)

式中：R为震中距；ω为圆频率； ${A}({\\text{λ}},\omega ,R)$ 为地震动加速度时程的Fourier幅值谱模型，其形式为：

$\begin{split} A({\\text{λ}},\omega ,R)= & \frac{{{A_0} \cdot \omega \cdot {{\rm{e}}^{ - K\omega R}}}}{{\sqrt {{\omega ^2} + {{({1 / \tau })}^2}} }} \cdot \\& \sqrt {\frac{{1 + 4{\xi _{\rm{g}}^2}{{({\omega / {{\omega _{\rm{g}}}}})}^2}}}{{{{[1 - {{({\omega / {{\omega _{\rm{g}}}}})}^2}]}^2} + 4{\xi _{\rm{g}}^2}{{({\omega / {{\omega _{\rm{g}}}}})}^2}}}} \end{split}$

(3)

${\varPhi }({\\text{λ}},\omega ,R)$ 为地震动加速度时程的Fourier相位谱模型，其形式为：

$\begin{split} {\varPhi }({\\text{λ}},\omega ,R)= & \arctan \left(\frac{1}{{\tau \omega }}\right) - R \cdot \ln [a\omega + 1000b +\\& 0.1323\sin (3.78\omega){\rm{ + }}c {\rm{cos}}(d\omega)] \end{split}$

(4)

式中： ${\\text{λ}}=[{A_0},\tau ,a,b,c,d,{\xi _{\rm{g}}},{\omega _{\rm{g}}}]$ 为基本物理参数向量，A₀为震源幅值系数，反映震源强度的影响，τ为Brune震源系数，反映震源断裂的时间过程特性，a、b、c和d为等效频率-波数关系的参数，ξ_g和ω_g分别为局部场地的等效阻尼比和等效卓越圆频率。

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

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

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

$\begin{split} a(t)= &- \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {{{A}}({{Z}},\omega ,R)} \cdot \cos [\omega t +\\& {{{\varPhi}} }({{Z}},\omega ,R)]{\rm{d}}\omega \end{split} \quad\qquad$

(5)

式中， ${{Z}}=[{{\\text{λ}}_{\rm{M}}},{{\\text{λ}}_{\rm{A}}}]$ 为主震与余震基本物理参数向量，即：

$\begin{split} {{Z}}=&[{A_0}_{\rm{M}},{\tau _{\rm{M}}},{a_{\rm{M}}},{b_{\rm{M}}},{c_{\rm{M}}},{d_{\rm{M}}},{\xi _{\rm{g}}}_{\rm{M}},{\omega _{\rm{g}}}_{\rm{M}},{A_0}_{\rm{A}},\\&{\tau _{\rm{A}}},{a_{\rm{A}}},{b_{\rm{A}}},{c_{\rm{A}}},{d_{\rm{A}}},{\xi _{\rm{g}}}_{\rm{A}},{\omega _{\rm{g}}}_{\rm{A}}] \end{split}$

(6)

考虑其参数间的相关性， ${{Z}}$ 可改写为：

${{Z}}=[{{{A}}_0},{\\text{τ}},{{a}},{{b}},{{c}},{{d}},{\xi _{\rm{g}}},{\omega _{\rm{g}}}]\;\;\;\qquad\qquad\qquad$

(7)

式中， ${{{A}}_0}{\text{、}}{\\text{τ}}{\text{、}}{{a}}{\text{、}}{{b}}{\text{、}}{{c}}{\text{、}}{{d}}$ 为描述主余震相关关系的六个二维随机变量。由于其不可观测性，参数值无法直接通过观测得到，可通过地震动反演由参数识别得到。

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

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

$\begin{split} a({r_l},t)=& - \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {A({\\text{λ}},\omega ,R)} \cdot \exp \left( { - \frac{{{\alpha _0}\omega {r_l}}}{2}} \right) \cdot\\& \cos \left[\omega t + {{\varPhi }}({\\text{λ}},\omega ,R) - \frac{{\omega {r_l}}}{{{c_{\rm{g}}}}}\right]{\rm{d}}\omega \end{split}$

(8)

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

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

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

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

$\begin{split} a({r_l},t)=& - \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {{{A}}({{Z}},\omega ,R)} \cdot \exp \left( { - \frac{{{\alpha _0}\omega {r_l}}}{2}} \right) \cdot \\& \cos \left[\omega t + {{{\varPhi}} }({{Z}},\omega ,R) - \frac{{\omega {r_l}}}{{{c_{\rm{g}}}}}\right]{\rm{d}}\omega \end{split}$

(9)

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土层厚等效剪切波速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)	>500	250~500	150~250	<150
等效剪切波速V_S30/(m/s)	>596	278~596	158~278	<158
地震动数量/组	130	636	238	34

下载: 导出CSV

| 显示表格

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

表 2 筛选得到的实测序列型地震动

Table 2. Selection of measured sequential ground motions

地震事件	主震			余震
地震事件	发生时间	Mw	数量	发生时间	Mw	数量
L′Aquila	2009.04.06 01:33	6.3	64	2009.04.06 17:47	5.6	64
Livermore	1980.01.24 19:00	5.8	12	1980.01.27 02:33	5.4	12
Mammoth Lakes	1980.05.25 16:34	6.1	6	1980.05.25 16:49	5.7	6
Chi-Chi	1999.09.20 17:47	7.6	480	1999.09.20 17:57	5.9	480
Chalfant Valley	1986.07.20 14:29	5.8	10	1986.07.21 14:42	6.2	10
Northridge	1994.01.17 12:31	6.7	62	1994.03.20 21:20	5.3	62
Umbria Marche	1997.09.26 09:40	6.0	18	1997.10.06 23:24	5.5	18
East Japan	2011.03.11 14:46	9.0	24	2011.03.11 15:15	7.7	24
Alaska	2002.10.23 03:30	6.7	38	2002.11.03 22:12	7.9	38
Imperial Valley	1979.10.15 23:16	6.5	28	1979.10.15 23:19	5.0	28
New Zealand	2010.09.03 16:35	7.0	154	2011.02.21 23:51	6.2	154
Whittier Narrows	1987.10.01 14:42	6.0	134	1987.10.04 10:59	5.3	134
Kumamoto	2016.04.14 21:26	6.2	8	2016.04.16 01:25	7.0	8
SMART1 Events	1981.01.29 04:51	5.9	54	−	−	−
	1983.09.21 19:20	6.5	70	−	−	−
	1985.06.12 17:22	5.8	66	−	−	−
	1986.05.20 05:25	6.3	72	−	−	−
	1986.11.14 21:20	7.3	72	−	−	−

下载: 导出CSV

| 显示表格

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

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

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

图 1 参数识别流程图

Figure 1. Flowchart of parameter identification

下载: 全尺寸图片幻灯片

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

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

3) 将2)中的局部场地参数ξ_gM和ω_gM代入余震幅值谱模型中，识别余震模型中的参数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) 对地震动场参数作识别，得到α₀和c_g。

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

$J=\min {\sum {\left| {A\left( {{{Z}},\omega ,R} \right) - \tilde A( \omega )} \right|} ^2}$

(10)

式中： ${{A}}({\\text{λ}},\omega ,R)$ 为本文提出的Fourier幅值谱模型； $\tilde A(\omega)$ 为真实地震动幅值谱值。

与幅值谱不同，相位谱往往类似于白噪声，难以根据真实地震动的相位谱识别其参数。然而，地震动相位差谱却有较为明显的特征，相位差一般符合正态分布或近似正态分布。基于此，相位谱对应的a、b、c、d由以使模型的相位差分布逼近真实地震动相位差分布为原则，采用遗传算法作为识别方法，目标函数可表述为：

$J=\min \sum {{{\left| {g\left( {\Delta {\phi _k}} \right) - \tilde g\left( {\Delta {\phi _k}} \right)} \right|}^2}}$

(11)

式中： $g(\Delta {\phi _k})$ 为本文提出的相位谱模型的相位差分布； $\tilde g(\Delta {\phi _k})$ 为真实地震动相位谱的相位差分布； $\Delta {\phi _k}={\phi _{k + 1}} - {\phi _k}$ 为频域上相邻两个相位角的差。

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

表 3 遗传算法参数设置

Table 3. Settings of Genetic Algorithm

参数类型	参数设置值
创建初始种群方法	二进制编码法
种群规模	200
适应度函数	排列尺度变换函数
选择操作函数	随机均匀分布选择法
变异函数	均匀变异函数
变异概率	0.04
交叉函数	离散重组交叉方式
交叉概率	0.8
终止条件	最大进化代数为500代
函数精度	10⁻⁸

下载: 导出CSV

| 显示表格

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

${c_{\rm{g}}}=\frac{{\Delta {r_l}}}{{\Delta t}}$

(12)

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

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

表 4 概率密度函数及其对应参数

Table 4. Probability density functions and parameters

随机变量	场地	主震			余震
随机变量	场地	函数类型	P₁	P₂	函数类型	P₁	P₂
A₀	Ⅰ	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.685	Weibull	0.669	0.479
	Ⅱ		0.661	0.359		0.634	0.380
	Ⅲ		0.726	0.389		0.848	0.392
	Ⅳ		0.811	0.483		0.738	0.160
a	Ⅰ	Lognormal	1.886	0.550	Lognormal	2.150	0.462
	Ⅱ		1.812	0.659		2.054	0.544
	Ⅲ		1.461	0.721		1.834	0.643
	Ⅳ		1.046	0.739		1.774	0.638
b	Ⅰ	Lognormal	2.012	0.498	Lognormal	2.312	0.458
	Ⅱ		1.933	0.557		2.220	0.528
	Ⅲ		1.746	0.631		1.940	0.587
	Ⅳ		1.421	0.514		1.790	0.637
c	Ⅰ	Weibull	1.272	1.231	Weibull	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.858
d	Ⅰ	Weibull	1.682	1.591	Weibull	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

下载: 导出CSV

| 显示表格

表 5 概率密度函数及其对应参数

Table 5. Probability density functions and parameters

随机变量	场地	函数类型	P₁	P₂
ξ_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.610	0.851
	Ⅳ		11.409	2.827
α₀	Ⅱ	Normal	1.014×10⁻⁴	2.214×10⁻³
c_g	Ⅱ	Lognormal	1.414	0.278

下载: 导出CSV

| 显示表格

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

下载: 全尺寸图片幻灯片

图 2 A_0M和A_0A的边缘概率密度函数

Figure 2. Marginal PDFs of A_0M and A_0A

下载: 全尺寸图片幻灯片

${\rm{BIC}}= - 2\sum\limits_{i=1}^N {\ln f({x_i};p,q,r)} + {k_1}\ln N$

(13)

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

3.3 主余震的相关性分析

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

图 3 A_0M和A_0A的散点图

Figure 3. Scatters of A_0M and A_0A

下载: 全尺寸图片幻灯片

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

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

图 4 识别最优Copula函数的流程图

Figure 4. Flowchart for identifying optimal Copula function

下载: 全尺寸图片幻灯片

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

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

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

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

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

表 6 最优Copula函数及相关参数

Table 6. Optimal Copula functions and relevant parameters

随机变量	场地	函数类型	相关参数θ
A₀	Ⅰ	Plackett	9.896
	Ⅱ		18.808
	Ⅲ		26.833
	Ⅳ		7.156
τ	Ⅰ	Clayton	2.214
	Ⅱ		2.601
	Ⅲ		1.287
	Ⅳ		1.951
a	Ⅰ	Plackett	2.710
	Ⅱ		3.592
	Ⅲ		2.840
	Ⅳ		3.057
b	Ⅰ	Clayton	0.916
	Ⅱ		1.230
	Ⅲ		1.262
	Ⅳ		1.442
c	Ⅰ	Clayton	0.467
	Ⅱ		0.493
	Ⅲ		0.407
	Ⅳ		1.400
d	Ⅰ	t (ν=4)	0.211
	Ⅱ		0.255
	Ⅲ		0.271
	Ⅳ		0.659

下载: 导出CSV

| 显示表格

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

图 5 实测数据与模拟数据对比

Figure 5. Comparisons of real data and simulation

下载: 全尺寸图片幻灯片

4 序列型地震动的验证

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

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

下载: 全尺寸图片幻灯片

图 6 实测地震动与模拟地震动的均值反应谱对比

Figure 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 实测单点时程与模拟时程的对比

Figure 7. Comparisons of measured and simulated records

下载: 全尺寸图片幻灯片

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

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

图 8 实测台阵时程与场模型模拟时程对比

Figure 8. Comparisons of measured and simulated records

下载: 全尺寸图片幻灯片

5 结论

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

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

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

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

期刊类型引用(8)

1.	高嘉鸿, 庞锐, 宰德志, 徐斌. 随机主余震作用下AP1000核电厂可靠度研究. 振动工程学报. 2025(07) 百度学术
2.	周洲, 韩淼, 吕大刚, 曾立静, 于晓辉. 考虑多损伤指标的RC框架结构主余震易损性分析. 工程力学. 2025(08) 本站查看
3.	金聪鹤，钱永久，徐望喜，黎璟. 考虑抗力非平稳劣化和荷载相关性的桥梁可靠性分析. 工程力学. 2024(02): 202-212 . 本站查看
4.	徐斌，荣卓，庞锐，周扬. 高堆石坝面板随机动力分析和抗震效果研究. 水力发电学报. 2023(01): 30-39 . 百度学术
5.	殷京科，李典庆，杜文琪. 主余震序列作用下边坡位移响应及地震动参数相关性分析. 工程力学. 2023(03): 44-53 . 本站查看
6.	刘懿俊，施秋华，谭飞，谢静峰，焦玉勇. 基于Copula理论的深圳市地面塌陷致灾因素的联合概率分析. 武汉大学学报(工学版). 2023(05): 549-556 . 百度学术
7.	朱志辉，刘禹兵，高雪萌，周高扬，余志武. 概率密度演化方程差分格式的计算精度及初值条件改进. 工程力学. 2022(11): 13-21 . 本站查看
8.	许泽坤，陈隽. 非线性结构地震响应的神经网络算法. 工程力学. 2021(09): 133-145 . 本站查看

圆环对撞问题的理想刚塑性小变形分析及讨论

计量