在工程地震研究中如何得到场地相关的地震动一直受到研究者的广泛关注。人工合成地震动可为地震动记录匮乏地区的地震危险性分析和结构抗震分析、性能评估提供输入。1947年Housner^[1]从工程研究的角度将地震动视作大小一定、随机到达的脉冲的叠加来模拟地震动开始，国内外许多学者对地震动合成方法进行了研究，也提出了多种满足部分地震动统计特征，如功率谱、傅氏谱、反应谱等的人造地震动方法。目前地震动的合成方法主要分为三种：一是以直接运用地震动三要素与震源、路径和场地的关系为基础的经验方法，如比例法、经验格林函数法^[2-3]；二是采用随机振动模型的地震动合成方法，如脉冲法、随机点源法、随机有限断层法^[4-6]；三是采用考虑强度和频率非平稳的地震动合成方法，如各种时频非平稳合成、相位谱、相位差谱、三角级数法等^[7-12]。

在工程中应用最广泛的则是拟合目标反应谱的地震动时程合成方法，以合成时程的反应谱与给定的目标反应谱相等或相近为控制条件。但是，仅以匹配目标谱为控制条件会一定程度上得到不唯一的合成地震动时程，而且这些合成地震动由于缺乏更多的约束条件，基于此进行时程分析可能会增加结构响应的离散性和不确定性。究其原因，反应谱虽然可以反映地震动幅值和频谱的特性，但是对地震动三要素中的持时并没有体现。作为地震动三要素之一的持时对于结构的响应有重要的影响^[13-14]，Bommer等^[15]认为地震发生过程中的强震持续时间对于基础材料和结构响应有着重要的影响，特别是在强度和刚度退化时。因此，地震危险性分析也应该包括对强震持时的估计。Mahin^[16]也指出地震动持续时间对结构的弹性变形和耗能需求有显著影响，尤其是对于相对薄弱的短周期结构。Cabanias等^[17]指出使用同时考虑振幅和持时的响应参数可以更好的估计结构的潜在损伤。

地震动的幅值、频谱及持时等特征主要受到震源、路径和场地条件影响和控制。研究表明不同地震构造区域的地震动的幅值、频谱和持时存在区域性差异^[18]，因此，在合成和模拟场地相关地震动时，应该通过约束幅值、频谱和持时这些工程相关的主要参数，以期使得合成的地震动中包含地震动的区域特征。基于此，本研究的目的是建立同时匹配设定区域地震动目标反应谱和持时参数的地震动合成方法。一方面，这样合成的地震动可以综合的考虑到地震动三要素中的幅值、频谱和持时特征；另一方面，增加地震动持时的合成控制条件可以在实现对目标反应谱匹配的基础上，为研究输入地震动持时对结构响应的影响提供基础。

近年来机器学习方法被多个领域广泛应用^[19-21]。地震学是一个数据密集型领域，目前正在经历着大数据界中体量、种类、速度的快速变化^[22-24]，这些变化为基于数据驱动的机器学习在地震学中的应用提供了条件。Khosravikia等^[25]开发了一种人工神经网络用于美国俄克拉荷马州等地区的自然和诱发地震的地震动预测模型；Kotha等^[26]在对日本活动性浅源地震的地震动预测模型回归分析的基础上，利用无监督学习中的聚类技术对代表各个场地的经验放大函数的估计结果进行了分类，进而提出了一种修正的数据驱动的场地分类方法；Alimoradi等^[27]开发了一种将一组稀疏的被称作特征地震记录的标准正交基向量进行线性组合合成地震动的方法，其组合系数是匹配由高斯过程回归得到设定情景下的反应谱确定。这些研究表明，机器学习具有很好的数据分析处理能力，可以为解决复杂的地震动模拟和合成问题提供新的思路。

为了在地震动合成中考虑地震动幅值、频谱和持时的特征^[28]，将地震动数据与机器学习相结合，运用智能算法寻优求解实现。本文首先应用主成分分析^{[27, 29-31]}从区域地震动数据库中提取具有本区域地震动特征的地震动母波，为了实现合成的地震动能够同时考虑地震动幅值、频谱和持时特征，选择同时匹配目标反应谱和目标持时。采用多目标遗传算法求解地震动母波的线性组合系数并且给出合成地震动的评价标准来控制合成效果，从最优解集中筛选出最合理的母波组合系数用于合成地震动。以2019年中国四川长宁M_s6.0地震的主余震数据作为区域地震动数据库，对本文方法进行了可行性和有效性的验证，最后将本文方法同三角级数法进行了对比。

1   基于PCA和MOP的地震动合成方法

1.1   主成分分析方法

为了在合成地震动时考虑地震动区域特征，需要对大量原始地震动数据的特征进行提取，这实质上就是数据降维。数据降维主要有投影和流行学习两种常用的方法。主成分分析(Principal Component Analysis, PCA)是基于投影的方式来实现线性数据降维和去除相关性的一种方法，其在数量地理学、数学建模和数理分析等学科领域中都有广泛应用。本文将PCA算法引入对地震动数据库的降维处理中，提取表征数据库主要特征成分的特征地震动母波，将这些特征地震动母波用于合成地震动时程。PCA算法的核心是先识别出最接近数据的超平面，然后，将数据投影到该超平面之上，进而实现降低数据向量的维度并去除掉各个分量之间的相关性。目标超平面的识别可以转化为如何得到投影矩阵。该问题是通过优化目标函数算法来实现的，这与其他机器学习算法的思路是一致的。

以一个包含N条地震动记录的数据集为例来说明PCA方法的应用原理。该地震动数据集可视为一个拥有N个维度的数据空间，则可以通过N个两两互相垂直坐标轴上的单位向量的线性组合来表示数据空间中的任意一条地震动记录。将整个数据集投影到第i个的单位向量所在坐标轴上后就称为该投影结果是第i个主成分，即第i条特征地震动母波，且这些用于投影的坐标轴排序是按照数据空间向本轴单独投影后保留数据差异性的大小确定。具体对地震动数据集应用PCA算法的降维流程如下：

1)计算N维实际地震动记录数据集的均值向量，将所有向量减去各自的均值向量，即白化处理；

2)计算数据集的散度矩阵S(或协方差矩阵C)；

3)对散度矩阵进行特征值分解，利用奇异值分解(Singular Value Decomposition, SVD)，得到所有的特征值和特征向量；

4)将所得特征值按从大到小排序，保留最大一部分的K个特征值对应的特征向量，以它们为列，形成降维所需的投影矩阵；

5)将白化后的数据库矩阵右乘投影矩阵，得到降维后全新的正交特征地震动数据，即为特征地震动母波。

对于选择降维后的维度K(即特征地震动母波个数)，可以用SVD分解时产生的散度矩阵S来表示。利用式(1)计算可得满足误差条件最小的K值，若其中t值取0.15，则代表了该PCA算法保留了原数据85%的主要信息，从而可以确定该允许误差下的提取的最少特征地震动母波数。具体t的取值可根据对误差的要求确定，由此选择不同的K值。

1.2   合成系数的控制条件

在得到特征地震动母波的基础上，需要考虑合成地震动时各条母波的对应系数的取值问题。不同于传统的合成地震动时程只考虑目标反应谱的匹配，本文从多参数的角度约束合成地震动时程。工程输入中地震动的特征参数一般包括幅值、频谱和持时。地震动的反应谱可以体现地震动的幅值和频谱两方面的特征，而地震动持时是结构地震反应分析中的重要参数。在结构地震反应的弹性阶段，循环往复的地震作用可能引起构件的低周疲劳，而在结构进入非弹性阶段后，持时的长短可能对结构反应、残余变形和倒塌存在重要影响^[32]。因此，需要引入地震动持时作为约束合成的地震动时程的控制条件。

1.2.1   合成地震动的评价函数

引入合成地震动的评价函数来判定最后的合成结果是否达到目标要求，以得到同地震动三要素匹配良好的合成地震动时程。以$S_a^*( {{T_j}} )$表示由目标地区地震动衰减关系得到的目标反应谱，其中：T_j表示在谱形匹配过程中所选用的自振周期控制点；而${S_a}\left(\displaystyle\sum\nolimits_{i = 1}^{{N}} {k_j}{e_j},T\right)$代表由提取的地震动母波线性组合得到的新地震动时程所计算出的反应谱，其中${\displaystyle\sum\nolimits_{i = 1}^{{N}} {{k_j}e} _j}$为提取的N条地震动母波同各自系数相乘后得到的新地震动时程。对于目标反应谱和合成地震动时程所得反应谱，计算拟合误差绝对值之和。

为了对持时进行约束，首先需要选择适合目标地区的地震动持时预测方程，进而得到目标地震动持时，记为d_5-95，新合成地震动的持时记为D_5-95。持时项误差绝对值为：

1.2.2   合成地震动的评价标准

为了对上述合成地震的两项误差给出一个合理的评价标准。对于地震动预测模型(GMPE)，均采用式(4)和式(5)中的评价标准来筛选出合成的地震动时程。

在评价标准中，以地震动反应谱和持时预测方程中给出的标准差σ为评价条件，合成反应谱在各控制周期点与目标反应谱谱值之差的绝对值小于一倍的反应谱衰减关系中的标准差时，视为合成地震动反应谱与目标反应谱良好匹配；对于地震动持时也以一倍的标准差为控制值来判断合成地震动持时的匹配效果。

1.3   确定合成系数

为了将PCA算法提取出的特征地震动母波线性组合得到合成地震动，需要确定组合系数。由于目标评价函数增加，传统上用于求解评价函数最优化解的单目标优化算法(SOP)已不再适用，而选择多目标优化算法(MOP)中最常用的多目标遗传算法来求解合成系数。其核心是协调各个目标函数之间的关系，从而得到使各子目标函数都尽可能达到最优的最优解集，这不同于单目标优化的解通常唯一，多目标优化问题的解是一组均衡解集。在应用算法得到一系列合成系数解集后，通过设定的评价函数指标筛选出最满意的一组合成系数。

多目标遗传算法中的带精英策略的非支配排序遗传算法(NSGA-II)是2000年Deb等^[33]为改进NSGA而提出的。由于改进后的二代算法引入了快速非支配排序法、精英策略、采用拥挤度和拥挤度比较算子使得算法的计算复杂度大大降低，计算结果在Pareto域中均匀分布，维持了种群的多样性。具体算法流程如下，并以流程框图1说明该优化算法对于合成控制函数的优化过程：

1)随机产生设定规模N的初始种群，经非支配排序后通过遗传算法的选择、交叉、变异后得到第一代子代种群；

2)从第二代开始，将父代与子代种群合并，进行快速非支配排序，同时计算每个非支配层中的个体的拥挤度，根据非支配关系以及个体拥挤度选取合适的个体组成新的父代种群；

3)最后通过遗传算法的基本操作产生新子代种群：循环上述3个步骤，直到达到终止代数，输出最优解集。

图  1  多目标遗传算法流程图

Figure  1.  Flow chart of multi-objective genetic algorithm

下载: 全尺寸图片 幻灯片

1.4   地震动合成整体流程

本文提出的合成方法是应用PCA技术从区域地震动数据库中提取地震动母波，在合成中考虑了对目标反应谱和显著持时进行同时匹配，使用多目标遗传算法确定出母波的线性组合系数得到最终的合成地震动。合成的整体流程见图2。

图  2  地震动合成流程图

Figure  2.  Flow chart of ground motion synthesis method

下载: 全尺寸图片 幻灯片

2   方法验证

2019年6月17日在四川省东南部的宜宾市长宁县发生了M_s6.0级地震。在震后的几日内，震源附近接连发生了多次的强余震，我国强震台网观测到了此次群震的大量地震动^[34]。对于长宁地震这类数据特征明显且记录丰富的数据，应用PCA技术提取的特征母波可以在满足设定的误差条件下，代表整个主余震地震序列中的主要特征信息。通过对长宁地震动的合成和比较可以验证本方法的可行性和有效性。

2.1   地震动数据库

本文收集了此次长宁地震的主震及其4级以上余震的记录作为地震动数据库，具体的地震信息见表1。数据库包含9次地震事件中，来自69个台站的286条水平向地震记录，并采用Boore等^[35]提出的方法对地震动数据进行了滤波和基线调整。

表  1  长宁主余震信息

Table  1.  Changning mainshock-aftershock events

地震名称	发震时间	震中位置		震源深度/km	震级Ms
		北纬/(°)	东经/(°)
长宁主震	2019-06-17	28.340	104.900	16	6.0
长宁余震	2019-06-17	28.430	104.769	16	5.1
长宁余震	2019-06-18	28.389	104.849	10	4.1
长宁余震	2019-06-18	28.389	104.870	10	4.2
长宁余震	2019-06-18	28.379	104.870	14	4.5
长宁余震	2019-06-18	28.370	104.889	17	5.3
长宁余震	2019-06-22	28.430	104.769	10	5.4
长宁余震	2019-06-23	28.389	104.819	14	4.6
长宁余震	2019-06-24	28.440	104.800	10	4.1
注：表中所列地震事件信息来源于中国地震台网中心(CENC)。

下载: 导出CSV 

| 显示表格

2.2   地震动母波的提取

首先按式(1)计算确定提取的地震动母波数目，为了使地震动母波能够保留数据集中95%的主要信息，t取值为0.05，通过计算得到满足误差条件的最小K值为11。在数据库中提取按照保留信息占比大小进行排名后的前11条地震动母波。表2中给出了提取的地震动母波对原始数据库信息保留占比情况，图3给出了排名前八位的地震动母波时程和时频特征，从图3中可以看出提取到的地震动母波与实际地震动记录的时域和频域非平稳特性相似，这一特征使地震动母波适合作为一组基向量用于地震动合成。

表  2  提取的地震动母波保留信息占比

Table  2.  Proportion of information retained by extracted ground motion mother waves

母波编号	1	2	3	4	5	6	7	8	9	10	11
${S_{ii}}$	0.363	0.203	0.104	0.097	0.069	0.039	0.030	0.025	0.010	0.007	0.004
$\displaystyle\sum\nolimits_{i = 1}^K {{S_{ii}}}$	0.363	0.566	0.670	0.767	0.836	0.875	0.905	0.930	0.940	0.947	0.951

下载: 导出CSV 

| 显示表格

图  3  排名前8的地震动母波

Figure  3.  Top 8 mother wave ground motions

下载: 全尺寸图片 幻灯片

2.3   地震动预测模型

尽管从实际地震动数据中提取而来的地震动母波包含了区域地震动的特征信息，但是直接采用提取的母波合成地震动并不能考虑震级、距离和场地等因素的约束，为了考虑这些区域地震相关因素对地震动幅值、频谱和持时的影响，需要在合成地震动时程对地震动的工程相关参数加以约束，使得合成地震动能够与区域实际地震动特征很好的匹配。

四川长宁地震发生于我国西部，因此需要选取适合于此地区的中国西部地区地震动模型^[36]，以及基于我国强震动数据建立的地震动持时预测方程^[37]。通过设定地震情景，利用反应谱和持时地震动模型分别计算得到目标反应谱和目标持时作为对长宁地震进行地震动合成的约束条件。选用的反应谱地震动模型和地震动持时模型的具体形式如下：

图4(a)~图4(d)中分别给出了长宁地震记录在0.01 s、0.20 s、1.00 s三个周期点处的反应谱值和对应地震动模型得到的衰减曲线的对比，以及实际地震动数据的5%~95%显著持时同持时预测方程所得预测值的对比。从图4中可以看出，选取的两个地震动模型可以较好预测本区域的实际地震动特征，也反映出了作为合成约束条件的目标反应谱和持时选取的合理性。

下载: 全尺寸图片 幻灯片

图  4  长宁地震地震动记录反应谱和持时与地震动预测模型对比

Figure  4.  Comparison between response spectra and durations of Changning earthquake and ground motion prediction models

下载: 全尺寸图片 幻灯片

2.4   合成目标地震动时程

为了验证合成方法的可行性，采用本文方法合成一条地震动并将其与实际地震动数据库中选取的一条地震动记录进行匹配，检验匹配结果是否合理。任意选取长宁地震中的一个台站51NXT，选取的地震动信息如下：震级5.4级，震源深度10 km，震中距46.92 km。然后，基于本方法进行合成，表3中给出了合成的地震动加速度幅值、反应谱、持时与实际地震动记录的匹配效果。图5中为合成的地震动与实际地震记录的匹配效果和误差分布图。图5(b)中横纵轴分别表示合成地震动的反应谱和持时与地震动模型所得目标值之间的误差值，具体计算采用式(2)和式(3)，图中分布的每个散点均代表最优解集中的一组母波线性组合系数，正方形散点代表着在评判标准下最终确定的组合系数。图5(a)、图5(c)中的合成地震动即为利用该散点所对应组合系数将母波线性组合得到的结果。结合表3和图5可以看出合成地震动的峰值加速度和反应谱都实现了同实际记录的良好匹配；更重要的是，合成地震动的90%显著持时与实际记录完美匹配，将合成地震的持时控制到了区域持时预测模型的目标水平。实现了考虑区域地震动预测模型约束的地震动合成，可以为工程结构动力时程分析提供体现区域特征的地震动输入。

表  3  合成地震动幅值、频谱和持时的匹配情况

Table  3.  Matching of amplitude, spectrum and duration of synthetic ground motion

匹配情况	震级Ms	震源深度/km	震中距/km	峰值加速度PGA/g	显著持时D5-95/s	反应谱差值/g
实际值	5.4	10	46.92	0.0318	11.045	−
合成值	5.4	10	46.92	0.0333	11.060	0.0784

下载: 导出CSV 

| 显示表格

图  5  合成地震动与实际地震动的比较

Figure  5.  Comparison between synthetic ground motion and actual ground motions

下载: 全尺寸图片 幻灯片

3   与其他方法的比较

三角级数法是工程上广泛应用的合成地震动的方法，为了将本文方法与其结果进行对比，采用三角级数法进行地震动合成。三角级数法合成中采用的地震动时程强度包络函数如式(8)所示。地震动包络函数参数取值根据霍俊荣^[38]建立的回归式(9)确定。

式中：c为衰减关系；t₁和t₂分别是控制平稳段首末时刻；t₃是地震动持时。

式中：M_s为地震震级；D为震中距；t_s为地震平稳段持续时间，${t_2} = {t_1} + {t_{\rm{s}}}$。

合成地震动所需的设定地震信息以及持时和反应谱匹配情况见表4，图6和图7分别给出了两种方法合成的地震动时程、持时曲线、多目标优化最优解集的误差分布、反应谱的谱形匹配误差。在同时满足对目标反应谱和持时的允许误差的解集中，筛选出与反应谱匹配最佳的一组合成系数，并以此系数合成本文方法的地震动。

表  4  设定地震下两种方法持时和反应谱匹配情况

Table  4.  Matching of duration and response spectrum of two synthesis methods under twoearthquake scenarios

震级 Ms	震源 深度/ km	震中距 R/km	场地 条件 VS30/(m·s−1)	目标 持时 D5-95/s	合成持时 D5-95/s	持时 差值/s	反应谱 差值/g
					本文/三角 级数法	本文/三角 级数法	本文/三角 级数法
5	10	10	360	8.469	8.230	0.239	0.0739
					2.215	6.254	0.1155
6	10	10	360	10.897	10.905	0.008	0.3420
					4.360	6.537	0.5022

下载: 导出CSV 

| 显示表格

图  6  M_s=5, R=10 km时两种方法合成结果对比

Figure  6.  Comparison of synthetic results of two methods for M_s=5 and R=10 km

下载: 全尺寸图片 幻灯片

图  7  M_s=6, R=10 km时两种方法合成结果对比

Figure  7.  Comparison of synthetic results of two methods for M_s=6 and R=10 km

下载: 全尺寸图片 幻灯片

对比两种方法合成的地震动，从图6(c)、图7(c)两次设定地震下本文方法的最优解集误差分布上来看，均存在较多组的系数组合可以使得合成地震动满足预先设定的评价标准。在图6(c)、图7(c)中也可以看出单一追求同反应谱的匹配会导致合成地震动的持时与地震动持时预测方程得到的目标持时偏差过大的情况出现，不能在合成地震动中很好地约束工程抗震中所关注的地震动持时，这对于将合成的地震动作为输入研究结构响应会造成影响。这一结果在图6(b)、图6(d)、图7(b)、图7(d)本文方法和三角级数法的合成结果对比中也可以看出。三角级数法考虑通过与反应谱的匹配可以实现对地震动幅值和频谱的约束，但不能对地震动持时参数进行约束。而应用本文的地震动合成方法在约束条件和评价标准的控制下，可以很好的将区域地震动预测方程所体现的区域地震动特征在合成结果中体现。

4   结论与讨论

本文提出一种基于机器学习合成考虑区域地震动幅值、频谱和持时特征的地震动的方法，基于长宁地震的主余震数据对方法进行了验证，并与三角级数方法进行了对比。本方法采用了从区域实际地震动中提取的特征地震动母波，保留了原始地震动的时频特征，合成的地震动考虑了区域地震动的持时特征。总结本文可得如下结论：

(1)将机器学习中的主成分分析应用于从大量实际地震动数据中提取地震动特征成分是可行的、有效的。以降维提取得到的少量特征地震动母波来代表整个数据库的整体特征，反映地震动的震源、路径和场地的综合效应。从母波时程和时频分析中可以看出，特征地震动母波与实际地震动的时频非平稳特征一致，可以将其用于合成包含本区域数据特征的目标地震动。

(2)不同于一般方法仅将目标反应谱作为匹配条件，本文方法可以实现对地震动幅值、频谱和持时进行同时约束，考虑了地震动的时频特征。对地震动三要素的综合考虑也能够更好得到符合实际地震动特征的合成结果。

(3)与三角级数法对比表明，仅匹配反应谱可能会得到与实际地震动持时相差很大的合成地震动，这会对结构动力响应结果引入较大的不确定性。