基于特征提取和谱匹配优化的区域地震动时程构建方法

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

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

在结构抗震设计中，重要建筑和不规则建筑通常需要通过动力时程分析得到结构的非线性反应特性，选择不同的地震动作为结构非线性动力分析的输入荷载会得到不同的结果，而地震动输入结构前通常需要调整其反应谱与设计反应谱一致[1]，因此，如何获得合理、真实且匹配目标反应谱的地震动时程，是减少结构响应不确定性的重要问题。目前相关地震动模拟方法主要有三角级数法[2−4]、经验格林函数法[5−6]、随机有限断层法[7− 8]和宽频带模拟方法[9]等，这些方法通常根据与震源相关的频域参数合成人工地震动，难以全面考虑区域地震构造背景、地壳介质和场地特征等引起的差异[10−11]。以三角级数法为例，由于受强度包线的控制，该方法生成的地震动时程具有很强的相似性，不能针对不同区域合成考虑区域性差异的地震动。虽然目前考虑场地项合成地震动有传递函数法和四分之一波长法[12]，但这要求详细、准确的区域地质资料。由于目前地震动模拟很难考虑地震发生全过程的所有因素，且合成的地震动缺乏区域性差异，这些方法构造的人工地震动与真实地震动之间有较大的差别，在用于抗震设计或评估时会增大输入地震动的不确定性。

基于实际强震记录合成反应谱匹配的地震动主要有频域和时域方法。频域方法[13− 15]首先利用一组特定形式的函数分解地震动时程，再通过调整这组函数的幅值使其反应谱与目标反应谱一致，最后根据新的幅值重组得到新地震动时程，重复以上步骤最终实现反应谱匹配的目的。与频域方法相比，时域方法的匹配精度更高，且合成地震动相较原始地震动具有改变程度小的优点[16]，因此时域上反应谱匹配方法受到了广泛的研究。KAUL[17]首先提出了时域上反应谱匹配的方法，之后LILHANAND 和TSENG[18]将这种反应谱匹配方法应用于匹配多阻尼反应谱，但是其采用的小波不能保持地震动的非平稳性。ABRAHAMSON[19]借鉴前者的算法提出了新的小波函数，虽然这种调整保留了原始地震动的非平稳性，但未考虑其对速度和位移时程的影响，需要通过基线校正来消除速度和位移时程的漂移。HANCOCK 等[20]通过调整小波函数消除了速度和位移时程的漂移，由于其小波形式中包含了基线校正，导致其无法得到解析解，并影响计算速度和效率。AL ATIK和ABRAHAMSON[16]提出了一种不产生速度和位移时程漂移的调整小波函数，从而使反应谱匹配过程中能够得到小波系数的解析解。张郁山和赵凤新[21]构造了一种增量位移小波函数，通过微分、时移以及线性调幅得到加速度小波函数。ADEKRISTI 等[22]利用Broyden 修正方法来降低反应谱匹配中由于不同小波之间相互影响所产生的非线性问题。HONG 和HUANG[23]每次只匹配单个周期的反应谱值，避免了由多个小波函数引起的非线性问题。ZENGIN 和ABRAHAMSON[24]将反应谱匹配方法应用于匹配瞬时功率谱，并得到能够同时匹配反应谱和瞬时功率谱的地震动，从而为合成近场脉冲型地震动提供了新方法。虽然这些时域上地震动合成方法的原理基本相同，但是用于调整反应谱的小波函数在改进，解决反应谱匹配中非线性问题的方法也在改进。

目前频域和时域合成地震动的方法均有缺陷。频域方法在调整函数幅值时改变了地震动的非平稳特征[16]，且反应谱匹配的收敛性较差。时域方法虽然能保留地震动的频谱特征，但由于引入反应谱峰值时间不变的假设来简化小波幅值的求解过程，将非线性问题用线性方程组近似求解，从而忽略了小波的相互作用以及地震动与其反应谱之间复杂的关系，导致这类方法存在迭代效率低、收敛性较差和匹配精度低等问题。此外，由于地震的产生受震源效应、传播路径和场地效应的影响[25−26]，且地震动存在随机不确定性[26]，频域和时域方法都无法合成出合理考虑区域特征的地震动。

相较于频域合成的地震动，时域方法合成的地震动的非线性动力反应与实际地震动的结果更为相似[27]，因此，本文基于时域方法合成地震动。为了实现考虑区域特征、降低不确定性且匹配反应谱的目标，本文提出了一套完整且高效的地震动合成新方法。该方法首先基于主成分分析算法得到包含目标区域地震动特征的种子地震动，再基于粒子群算法优化时域上反应谱匹配方法，用此方法调整种子地震动的反应谱，最终得到可模拟目标区域地震的合成地震动。为了验证该方法的可行性、合理性和有效性，本文基于中国四川地区的地震动数据，给出了该方法的实际应用案例与分析结果。

1 方法

1.1 主成分分析提取种子地震动

主 成 分 分 析[28](Principal components analysis,PCA)是一种降低数据维度的算法，该算法通常用于去除高维数据的噪声和提取重要特征[29]。高维数据在某一维度上特征的重要程度用数据投影到该维度后的方差衡量，方差越大则数据的特征越重要。PCA 算法将m维数据映射到n维上(m>n)，不仅要让数据在这n维上对应的方差最大，且这些维度之间相互正交，为此该算法首先找出使数据投影的方差最大的坐标轴，然后从与之前所有坐标轴正交的超平面中找出使数据投影的方差最大的坐标轴，以此类推，最终找出n个坐标轴并完成数据降维。将这一方法应用于地震工程，能从一组地震动数据中得到包含原始数据特征的种子地震动。

主成分分析方法得到种子地震动的具体步骤如下：

1)将m条地震动数据按行组成矩阵X。

2)对矩阵X按行进行零均值化处理。

3)计算矩阵X的协方差矩阵C，并求出矩阵C的特征值与特征向量：

4)将矩阵C的特征值与特征向量按特征值从大到小排列，取前n个特征向量按顺序组成矩阵P，则矩阵P右乘矩阵X即可得到新地震动数据—矩阵Y。

5)矩阵Y中有n条种子地震动，种子地震动的主成分贡献率(Contributing Rate of Principal Component, CRP)越大，其包含的原始数据的区域地震动特征[30](eigenquakes)越明显，主成分贡献率的计算公式如下：

其中，λi为矩阵C的第i个特征值(降序)。

时域上反应谱匹配方法的基本思路是在地震动的加速度时程上添加一组幅值很小的小波，使其反应谱值与目标反应谱值之间的差值减小，通过迭代计算实现反应谱匹配；求解小波幅值是一个非线性问题，为了用线性方程解决该问题，假设在地震动加速度时程上添加小波后，自振周期为Ti的单自由度结构的加速度达到峰值的时间ti不变。

反应谱匹配的方法如下：

1)计算周期为Ti时的反应谱差值，公式如下：

式中：n为周期的个数；Qi为目标反应谱值；Ri为地震动a(t)的反应谱值；Pi由结构加速度峰值确定，当加速度峰值大于零时，Pi=1，否则Pi=−1。

2)在a(t)上叠加小波函数[16]fj(t) 以消除反应谱差值Ej，此时Ri的改变量为cij：

式中：hi(t)为单自由度结构的加速度脉冲响应函数[16]；ti为加速度达到峰值的时间。

3)当周期为Ti时，令添加的所有小波对反应谱值的总改变量等于反应谱差值Ei，则：

1.2 反应谱匹配方法

其中，bj为小波函数fj(t)的幅值。

4)根据式(5)可求解出所有小波的幅值，则新地震动a1(t) 的计算公式如下：

其中，γ 为在0 到1 之间的缩放系数[16]。

用a1(t)代替a(t)，重复以上步骤，直到反应谱之间的差值小于设定差值为止。

1.3 基于粒子群算法的优化原理

目前时域上反应谱匹配方法并没有对式(6)中γ 的取值做出明确的说明，因此现有方法普遍存在迭代效率低、反应谱匹配精度低的问题。本文采用机器学习中的粒子群算法[31](Particle swarm optimization, PSO)来寻找每次迭代过程中γ 的最优解，由于PSO 算法能够共享所有粒子的信息来调整粒子在寻优过程中的速度和位置，采用该算法能快速找到问题的最优解，从而提高反应谱匹配方法的迭代效率与匹配精度。

粒子群算法的原理如下：

1)参数设置与种群初始化。

设置相关参数，例如：粒子的数目、搜索范围和速度范围以及迭代次数等，根据设置的参数随机确定粒子的位置和速度。

2)计算极值，确定最优解。

根据定义的适应度函数计算粒子的适应度值，找出每个粒子的个体极值，再从所有个体极值中找出群体极值作为最优解。

3)更新粒子的速度和位置，公式如下：

式中：i=1, 2, ···,N，N为选取的粒子总数；t为当前迭代步数；vi为粒子的速度；xi为粒子的当前位置；w为惯性因子(非负)，用于调节算法的寻优能力；r1和r2是0 到1 之间的随机数且相互独立；c1、c2是学习因子；pbesti为个体最优解；gbest为全局最优解。

本文的惯性因子采用线性递减权值策略，公式如下：

式中：g为迭代次数；w(g)为第g次迭代的惯性因子；Gk为最大迭代次数；wini为初始惯性权值；wend为最终惯性权值。

4)判断迭代次数与最优解是否符合要求。

当达到最大迭代次数或者最优解满足设定值时，迭代结束，否则返回步骤2。

在上述原理中，单个粒子的优劣要通过适应度函数来评价。本文构造了四种误差函数,其中任意一种均可作为适应度函数，从而满足不同情况下对反应谱匹配精度的要求。四种误差函数的计算公式如表1 所示。

表1 误差函数的定义Table 1 Definition of error function

式中：N为选取周期点的个数；Ti为第i个自振周期；Sac(Ti)为周期为Ti时新地震动的反应谱值，Sat(Ti)为周期为Ti时目标反应谱值。

本文应用粒子群算法求反应谱匹配过程中γ 最优解的流程图如图1。

图1 粒子群算法求最优解γ 的流程图Fig. 1 Flow chart of particle swarm optimization for solving optimal solution γ

2 地震动合成整体思路

由上述方法可知本文合成地震动的整体思路分为两步：首先从目标区域选取一组区域地震动，用PCA 算法从中得到包含区域地震动特征的种子地震动；再用基于PSO 算法优化的时域内反应谱匹配方法调整种子地震动，通过叠加小幅值的小波调整其反应谱与目标反应谱一致，最终得到合成地震动。整体流程图如图2 所示。

图2 本文的地震动合成流程图Fig. 2 The flow chart of ground motion synthesis

PCA 算法得到种子地震动的优点：1)降低了地震动的不确定性，地震动在时域和频域上的不确定性会导致结构地震响应也具有不确定性[32]，而基于最大方差理论的PCA 算法能从一组地震动中得到不确定性更低的种子地震动；2)具有区域特征，由于PCA 算法能找到一组数据的特征并得到包含此特征的新数据，因此，如果选择的一组地震动具有区域特征，则种子地震动也具有区域特征。此外，从目标区域的地震动数据库中挑选地震动应遵循震级、震中距和场地条件等与设定条件基本相符的原则，则用PCA 算法得到的种子地震动也满足这三个条件的约束。正是由于种子地震动具有上述两个优点且满足一定的约束条件，用它来模拟设定地震在目标区域可能产生的地震动是合理的。

由于种子地震动的反应谱一般与设定区域的目标反应谱不匹配，因此，还需要用优化的反应谱匹配方法调整种子地震动的反应谱与目标反应谱一致。由于该方法不会对种子地震动的时域和频域特征产生较大的改变，因此，最后合成的地震动可用于模拟设定地震在目标区域产生的地震动。

AL ATIK 和ABRAHAMSON[16]提出的反应谱匹配的改进方法为了保证迭代的收敛性，在单次匹配过程中使匹配周期从短周期逐渐扩大到全周期，且迭代过程中未优化缩放系数γ 的取值，因此，该方法的迭代效率和匹配精度相对较低。用PSO 算法优化该方法后，求解γ 的最优值能提高迭代的收敛性，同时，全周期匹配也避免了多次叠加小波降低匹配精度的问题，因此，本文的优化方法提高了迭代的收敛性、效率和匹配精度。

3 方法的验证

为了验证本文提出的基于实际记录的地震动时域合成方法的可行性，以四川地区的地震动为基础进行目标区域的地震动构建。本文的反应谱匹配方法仅考虑单一阻尼，下文算例中所有反应谱对应的阻尼比均为5%。

3.1 基于PCA 方法获取种子地震动

假定需要合成目标区域为四川地区、震级5 级、震中距30 km、场地条件为土层的地震动，为了体现区域地震动特征，选取34 条四川雅安余震的地震动记录，以此合成雅安地区的种子地震动，这些地震动记录的震级为4.5～5.4、震中距为25 km～35 km，详细的地震信息如表2 所示。

表2 区域地震动信息Table 2 Regional ground motion information

通过PCA 算法的分析计算，可依据主成分贡献率从大到小选取对应的种子地震动，本文选取对应CRP最大的4 条地震动作为种子地震动，依次将其记为地震动1～4，其加速度时程图如图3 所示。

图3 种子地震动的加速度时程Fig. 3 The acceleration time history of the seed ground motion

3.2 优化反应谱匹配方法合成地震动

为了使目标反应谱具有四川地区的地震特性，本文选择了基于实际数据得到的川滇地区地震动预测模型[33]，根据设定的震级和震中距即可得到目标反应谱(控制周期点见参考文献[33]，不含PGA)。如图4 所示，灰色曲线是目标反应谱，其他曲线是4 条种子地震动的反应谱。

图4 种子地震动反应谱与目标反应谱Fig. 4 Seed ground motion response spectrum and target response spectrum

为了保证收敛性，种子地震动需要适当的调幅和补零[16]。由于迭代效率的提高，本文设置迭代次数为20 次即可使反应谱差值降到很低的水平，实际应用中也可以预设合适的匹配精度。为了找出γ 的最优解，PSO 算法的粒子数和迭代次数不宜过低，本文设定粒子数为20 个、迭代次数为5 次。为了使合成地震动的反应谱与目标谱之间的差值更小，本文选择反应谱匹配差值的总体平均值(ERR1mean)作为粒子群算法的适应度函数。由种子地震动调整得到的4 条合成地震动的加速度时程如图5 所示。

图5 合成地震动的加速度时程Fig. 5 The acceleration history of the synthetic ground motion

图7 是4 条合成地震动在控制周期点上ERR1和ERR2的曲线图，可以看出在所有周期上ERR1不超过0.045 cm/s2，ERR2不超过0.3%，因此优化的反应谱匹配方法能得到匹配精度很高的合成地震动，至此合成地震动方法的可行性得到了验证。

图6 匹配误差随迭代次数增加的变化曲线Fig. 6 The variation curve of matching error with the increase of iteration number

3.3 合成地震动的结果分析

从图7 可以看出，合成地震动的反应谱与目标谱基本一致，但在匹配反应谱时种子地震动的时程需要调幅和叠加小波函数，因此要从时域和频域上分析此过程是否改变了地震动的区域特征。由于种子地震动调幅为初始地震动的过程不会改变区域特征，且小波函数均叠加在初始地震动上，因此分析比较初始地震动与合成地震动即可确定后者是否仍具有区域特征，为此选择种子地震动1 对应的初始地震动和合成地震动为本节的分析对象。

图7 合成地震动的匹配误差Fig. 7 Matching error of the synthetic ground motion

如图8 所示，初始地震动和合成地震动的加速度和速度时程非常相似，两者的位移时程基本相同，且最大位移相差不大；从图9 可以看出初始地震动和合成地震动的归一化Arias 强度曲线基本一致，且两者的显著持时D70和D90相差不超过0.024 s，故两者的能量分布是一致的，因此这两条地震动在时域上的特征基本没有差别。

图8 初始地震动和合成地震动的时程Fig. 8 Time history of initial and synthetic ground motion

图9 初始地震动和合成地震动的归一化Arias 强度Fig. 9 Normalized Arias intensities for initial and synthetic ground motion

如图10 所示，两条地震动的傅里叶幅值谱基本一致。由于地震动是非平稳信号，因此还需要用短时傅里叶变换分析地震动的时频特征。图11分别是两条地震动的时频图，与图11(a)相比，图11(b)在2 s～5 s 处存在2 Hz～8 Hz 的低频成分减少、在3 s 附近存在大于10 Hz 的高频成分增加的现象，这是由于种子地震动在不同周期下达到加速度峰值的时间ti基本在2 s～5 s，使得叠加小波后这一段时程的改变较大，但从整体上可以看出两条地震动的时频特征基本相同，因此两者在频域上的特征基本没有差别。

图10 初始地震动和合成地震动的傅里叶幅值谱Fig. 10 Fourier amplitude spectra of initial and synthetic ground motion

图11 初始地震动和合成地震动的时频图Fig. 11 Time-frequency diagrams of initial and synthetic ground motion

本节从时域和频域上对合成地震动和初始地震动进行比较，得出两者的相关特征基本相同的结论，而初始地震动保留了种子地震动中的区域特征，因此合成地震动也具有设定区域的地震动特征，同时其反应谱与目标谱一致，所以用本文合成方法得到的地震动模拟目标区域可能产生的地震是合理的。

4 结论与讨论

本文基于实际地震动记录并结合机器学习算法得到具有区域特征且与目标谱一致的合成地震动，本文的主要结论如下：

(1)提出了一种时域上合成地震动的方法，该方法先由PCA 算法获得包含目标区域地震动特征的种子地震动，再根据PSO 算法优化的反应谱匹配方法调整种子地震动，使其反应谱与目标谱一致，最终得到可用于模拟目标区域地震的地震动。

(2)与现有的方法相比，用PCA 算法降维处理地震动数据时不仅降低了地震动的不确定性，同时保留了地震动的区域特征，因此相对于选择实际地震动记录或是人工地震动，用种子地震动模拟目标区域地震更为合理。

(3)用PSO 算法对目前已有的时域上反应谱匹配方法优化，从而在每次反应谱匹配过程中获取用于约束小波函数的最优缩放系数γ，此方法大幅提高了反应谱匹配过程中的迭代效率和匹配精度，并使地震动高精度匹配目标谱成为可能。

由于获取实际地震动特征需要目标区域实际地震动记录，对于强震数据不丰富的地区可从地震地质背景相似地区的地震动记录中提取种子地震动。本文考虑区域实际地震动特征的构建方法为解决面向区域的地震动场构建提供了参考。