基于PCA 及PSO 智能算法的地震动合成方法—以中国西部中强地震为例

胡进军，张 辉，靳超越，王中伟，胡 磊

(中国地震局工程力学研究所地震局地震工程与工程振动重点实验室，黑龙江，哈尔滨 150080)

随着城市化进程的加速，城市的规模越来越大，并且逐渐向地震区延伸，庞大复杂的城市基础设施增加了潜在的地震危险性。地震危险性分析和评估一直是地震工程研究的关键问题[1 − 3]，合理、准确的地震动输入是进行震前的地震风险评估和震后的地震灾害损失评估的前提。然而由于地震受到诸多因素的影响，震源、路径以及场地的未知性和不确定性一直是模拟和合成地震动中难以解决的问题。

确定工程场地地震动的方法有多种，基于概率的地震危险性分析(PSHA)方法是中常用的一种，其从概率的角度定量描述了地震作用[4 − 5]。基于地震动预测模型(GMM)[6 − 7]，利用PSHA 方法可以建立一套具有不同超越概率的危险谱，作为确定设防水准的基础。但是仅反应谱还不能满足结构的抗震设计和性能评估所需，峰值加速度(PGA)以及谱加速度(Sa)也常被用来表征地震动的强度指标[8 − 9]，一些其他的地震动参数如峰值速度(PGV)、峰值位移(PGD)、有效峰值加速度(EPV)、累积绝对速度(CAV)也在逐渐引起人们的重视。但是地震动信号是一种非平稳的随机过程，难以仅通过幅值或者能量参数就能够完全包含其复杂过程对结构的影响。因此地震动时程是结构非线性动力时程分析的必须，这个所需的时程可以通过人造或者对天然的地震动调幅得到[10]。除了基于PSHA 确定输入地震动的方法以外，基于物理过程的地震动模拟也是常用的方法，其可以考虑地震断层的破裂和传播过程，具有物理机制[11]，但是这种确定性的模拟方法依赖于震源模型、速度结构模型以及采用的模拟方法，需要准确的震源、路径和场地的模型，而目前工程上很难快速准确确定这些参数，因此，模拟的结果存在很大的不确定性。

近年来全球强震台站逐渐增多，区域强震数据日益丰富，我国的强震数据记录的数目已经超过了3 万条，而全球的强震数据已超过数十万条[12 − 13]。这为研究基于区域强震数据和数值模拟方法结合构建符合本区域地震构造、路径和场地特征的地震动提供了可能。但是面对数万计的地震动数据，如果采用人工的方法处理和挖掘将非常复杂和繁琐，在以往基于地震动数据的统计研究中，研究者一般先要筛选一部分数据，再针对选取的地震动数据进行分析，以减少计算量。随着计算机技术的发展以及机器学习理论的逐渐成熟，应用机器学习方法对海量地震动数据进行处理和挖掘成为可能。

为了探讨应用机器学习中的智能算法基于实际地震动合成目标地震动，本文采用机器学习方法中的主成分分析(PCA)算法，从目标区域地震动数据库中提取包含区域特征信息的地震动母波，同时基于目标地区的GMM 得到给定场地的加速度反应谱，基于特征母波和设计谱构建包含本地地震动特征的地震动时程。计算过程中，为了改进地震动时程的合成算法并提高计算效率，本文利用粒子群算法(PSO)快速找到母波地震动的权重系数，使得合成的地震动加速度反应谱与GMM 得到的目标谱误差最小，最终通过母波线性组合得到目标地震动时程。为了阐述本文方法的可行性和合理性，本文结合我国西部地区的中强强震数据开展研究。

1 区域强地震动数据

我国西部地区7 级以下的中强震数据相对丰富，但是7 级～8 级的大震事件仍然比较缺乏，在预测模型回归时缺少大震强震数据的约束，为了减少模型的不确定性。本文采用了2007 年−2019 年间四川以及周边省份的震级5.0 级～7.0 级的中强震地震动数据，以及基于此数据建立的中强震地震动预测模型。详细的地震信息如表1 所示，震中和台站位置见图1。数据库中包含了21 次地震中174 个台站的4551 条水平向地震动记录。对原始地震动记录进行了滤波和基线调整[14 − 15]。

表1 选取的西部地区的中强震Table 1 Selected earthquake events in west region of China

图1 台站和地震震中分布图Fig.1 Map of stations and earthquake epicenters

2 区域地震动信息的提取

研究表明：地震动由于受到地震构造、地壳结构和场地条件的影响，不同区域地震动可能具有不同的特征[15 − 17]。在模拟设定区域的地震动时需要考虑本区域的实际地震动的特征信息，因此，需要采用合理的方法从实际地震动中提取区域地震动的特征信息。

在数据挖掘和机器学习中，数据一般被表示为向量，与之类似，也可以把一条地震动记录视为1 列向量，那么n 条地震动记录就可视为n 列向量进而组合得到如下矩阵：

A=[ −→α1−→α2−→α3··· −→αn]m×n(1)

这样得到的合成地震动具备了原地震动所有的形状特征以及随机性特性。由于地震动数据库的记录数目较大，为了提交计算效率可以对原始的地震动矩阵A 进行降维简化。在线性代数中，一个内积空间的正交基是元素两两正交的基。在二维平面中，任意的二维向量都可以通过一组二维的正交基表示出来：

在三维平面中，任意的三维向量都可以通过一组三维的正交基表示：

当把一条条地震动向量视为一个个列向量时，那么肯定也存在一组正交基能够表示任意一条地震动所构成列向量。主成分分析[18](Principal component analysis，PCA)方法正好可以满足这方面的要求。它可以把数据降维，找出一组符合条件的正交基用于计算任意一条地震动记录。PCA 算法是一种对高维数据降维的方法，并将高维数据中重要的特征保留，去除噪声和不重要的特征。

若有一组如下形式的数据，应用主成分分析的具体步骤为：

首先，对数据进行标准化处理：

其次，计算相关系数矩阵：

接着，用雅克比方法求解出相关系数矩阵R的特征值λ 和特征向量，这里的特征向量就是正交基。

最后，选择重要的主成分，根据方差解释率即：

这里的方差解释率也称主成分贡献率，用于判断单个主成分所包含的原始数据信息的多少，方差解释率越大，所包含的原始信息越多因此本文基于此方法，从目标区域原始地震动数据库中提取含有本地地震动特征信息的母波。

2.1 基于PCA 算法的地震动母波提取

以数据为驱动提取母波的方法的具体流程如图2 所示。首先从区域原始数据库中应用PCA 提取一组标准的正交基向量，并要求这些提取出的正交基向量能够表征地震动时程序列的主要成分。

图2 主成分分析算法提取地震动母波的流程图Fig.2 The flow chart of ground motion mother wave extraction by PCA

应用PCA 算法提取的地震动主要成分的正交基在本文称为地震动母波，提取的地震动母波和原始地震动具有相同的时间采样频率。因此，合成的地震动可以由提取的地震动母波线性组合而成。

式中： ki为系数； ui为提取的地震动母波； n为提取的地震动母波的个数。地震动母波是数据矩阵组成的特征向量，然后根据特征值大小进行排序。

采用上述方法，可以从原始地震动数据库中提取n 条地震动母波，图3 给出了提取的4 条母波，从图3 中给出的地震动时程和傅里叶频谱特性可以看出，提取的母波与实际地震动记录特征非常接近。

2.2 地震动母波的合理性验证

为了验证基于PCA 算法提取的地震动母波合成地震动的合理性，本文以原始地震动数据库中的50 条近场(震中距R<30 km)数据为例提取母波，并进行合成和验证。选取近场数据进行母波提取和验证的原因是由于近场地震动的特征更显著、更加复杂，更具有代表性意义。

图3 从原始地震动数据库中提取的4 条地震动母波及其傅里叶频谱Fig.3 Four ground motion mother waves and their Fourier spectra extracted from the original ground motion database

为了使得提取的母波能够表征原始地震动数据库的特征，首先需要引入主成分累积贡献率的概念，主成分累积贡献率是选择有效主成分的重要依据，它是主成分的方差在所考察的随机变量的总方差中所占的比例；再引入累积方差解释率概念，即多个主成分方差所占的比例之和，它是通过主成分贡献率之和求得。当累积方差解释率比较高时，能够较好的代表数据库的特征。

为了确保能够充分的提取地震动数据库的特征，本文累积方差解释率取值为95%，即当累积方差解释率达到95%时，提取的含有本地地震动特征信息的母波能够很好表征原始地震动数据库的特征。图4 给出了累积方差解释率和母波地震动数量的关系，图中的拐点就是累积方差解释率取值为95%的点。根据图4 分析结果，当满足累积方差解释率为95%时，提取的国内近场地震动数据的母波数目为19 条。

图4 累积方差解释率和母波地震动数量的关系Fig.4 The relationship between the interpretation rate of cumulative variance and the number of ground motions of the mother wave

为了验证提取的母波的合理性，以提取的19 条地震动母波来合成近场数据库中的地震动。首先任意选取近场地震动数据库中的一条记录，计算该条地震动记录的反应谱，对提取的19 条地震动母波进行线性组合，可使得组合的新的地震动反应谱与之前选取的地震动反应谱误差最小，即可得到一条新的合成的地震动。图5 给出了实际地震动与合成地震动的反应谱和地震动时程的比较。从图5 的反应谱和地震动时程的比较中可以发现，PCA 算法提取的地震动母波能够很好的合成原始地震动数据库中的任意一条地震动记录，因此，PCA 提取的地震母波能够很好的表征原始地震动数据库的特征。

图5 实际近场地震动和应用PCA 算法提取的母波合成的地震动比较Fig.5 Comparison between the actual near-field ground motion and the synthetic ground motion by the mother waves extracted by PCA

3 地震动母波权重系数的确定

需要求解出方程的最优解，因此，引入粒子群算法。

3.1 粒子群算法(PSO)

粒子群算法(Particle Swarm Optimization, PSO)是由Kennedy 等[19]和Stefan 等[20]可以用于求解最优化问题，能够有效地实现计算机智能搜索和优化。该方法所求出的解是全局最优解而不是局部最优，它能够找出满足条件的一组 ki使得 S最小，具体的要点如下，流程见图6 所示。

1)参数的初始化。设置初始化参数，如：自变量 ki初始值，最大迭代次数，粒子的最大速度，粒子群的规模以及整个搜索空间。

2)个体极值以及全局最优解。个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。

3)更新速度和位置公式，即式(15)。

式中： ω为惯性因子，当取值较大时寻优能力强；C1和 C2为加速度常数； Pid为个体极值；Pgd为群体极值； Xid为粒子当前的位置； Vid粒子的速度；Maxgen 是迭代的次数。

4)设置迭代次数或者最小误差。

图6 粒子群算法求解权值ki 流程图Fig.6 Flow chart of PSO algorithm to solve weight ki

3.2 权重系数求解

为了使得地震动母波线性组合得到的新的地震动的反应谱与目标反应谱误差最小，图7 给出了基于PSO 算法[21 − 22]求解权重系数，以及基于地震动预测模型合成目标地震动的流程图。首先，选取本地震动数据库区域合适的地震动预测模型，应用PCA 算法提取地震动母波，通过地震动预测模型[23]得到的反应谱与组合地震动母波得到的新的地震动的反应谱匹配，再用PSO 算法快速求解权重系数。PSO 算法的具体参数参考了文献[24]，如表2 所示。

4 目标地震动的合成

图7 应用PCA 和PSO 算法合成地震动时程的流程Fig.7 Flow chart of simulation ground motion time history by PCA and PSO

为了验证本文提出的方法的可行性，分别对中国西部地区的四个设定地震场景下的不同场点进行地震动合成。设定震级、断层距以及场地条件如表3 所示。将设定震级、距离以及场地参数输入到本区域的地震动预测模型中，本文采用了文献[25]基于四川地区的中强震数据建立的地震动预测模型，与本文的研究区域一致。然后基于此模型对设定场点的地震动反应谱进行估计，通过组合母波得到的新的地震动时程并计算其反应谱，当计算的反应谱与地震动预测模型反应谱一致时，则得到最终的地震动时程，这是一个迭代过程。

通过粒子群算法寻优计算出的权重系数值如表4 所示。

图8 中给出了迭代次数和误差S 之间的关系，从图8 中可以看出在迭代到50 次时误差都已收敛，因此，针对本次地震动的计算模拟，可以取迭代次数为50。

图9 给出了合成的地震动时程的反应谱与预测模型得到的目标反应谱的比较，图中给出的分别是不同场点(R=10，30 km)和不同震级(M=5.5，6)的比较。从图9 中能够看出，通过PSO 智能算法求解出的地震动的反应谱能够较好的匹配地震动预测方程得到的目标反应谱。图10 给出了最终合成的地震动时程，从图中可以看出合成的地震动与实际地震动非常接近，具有随机性和非平稳性，包含了区域地震动的特征。因此，地震动母波的线性组合能够得到地震动数据库中的任意地震动数据，合成的目标地震动既匹配了目标谱，有能够很好地代表本区域实际地震动的特征。

表2 PSO 算法参数Table 2 Parameter of PSO

表3 设定地震信息和计算信息Table 3 Scenario earthquake and calculation information

表4 地震动母波的权重系数Table 4 Weight coefficient of the mother wave of the ground motion

图8 迭代次数与误差之间的关系Fig.8 The relation between the number of iterations and the error

图9 合成的地震动的反应谱与目标谱的比较Fig.9 Comparison between the response spectra of the synthesized ground motion and the object spectra

图10 机器学习方法合成的地震动时程Fig.10 Time history of ground motion synthesized by machine learning method

5 结论

为了研究考虑区域地震动特征信息的地震动合成方法，本文引入了机器学习中PCA 算法，从地震动数据库中提取有效的地震动母波信息，结合目标区域的地震动预测模型给出的特定场点的地震动反应谱，通过PSO 算法求解组合地震动母波的权重系数，使得合成反应谱与目标谱误差的最小，最终由地震动母波线性叠加得到目标地震动时程。通过上述研究可以得到如下结论：

(1)应用PCA 算法能够从地震动数据库中能够提取出代表性地震动母波，地震动母波能够合理表征地震动数据库的特性。

(2)应用PSO 算法能够快速高效求解地震动母波的组合权重，PSO 智能算法避免了应用穷举法求解权值，提升了计算速度。

(3)通过PCA 和PSO 智能算法，结合本区域实际地震动和预测模型来合成新的地震动时程，能够合理的包含区域实际地震动的特性，能够匹配目标地震动的频谱特征。

本文提出的方法考虑了区域实际地震动的特征，使得合成的地震动时程既包含了时程上的区域特征，又匹配了目标谱，满足了谱型上的一致性。采用PCA 和PSO 智能算法，提高了计算效率，满足了地震动合成时效性的需求，因此，可为未来面向工程的抗震性能评估提供合理的地震动时空分布场。