Prediction and analysis model for ground peak acceleration based on XGBoost and SHAP
-
摘要: 为建立一种不依赖土体本构模型，只依靠地震动和场地主要特征的地表加速度峰值预测方法，以日本KiK-net强震台网搜集到的3104组基岩和地表地震动记录为基础，通过特征选择筛选出6个特征参数，以输入地震动加速度峰值和输入地震动卓越频率表征输入地震动特性，以剪切波速达800 m/s时的土层埋深、场地基本周期、基岩剪切波速和地表剪切波速表征场地特性。采用XGBoost模型，构建基于6个特征参数的地表峰值加速度(PGA)预测模型。通过对比实测记录和一维数值模拟计算结果，表明本文建立的XGBoost模型预测结果稳定，能较好的预测PGA，训练集和测试集的决定系数均大于0.925，平均绝对百分比误差均在20%左右。同时引入SHAP对输入特征与预测结果之间的影响和依赖性进行分析，增强了模型的可解释性，同时也为预测结果提供了可靠性支撑。Abstract: In order to establish a prediction method for the ground peak acceleration (PGA) that does not depend on the soil constitutive model but only on the ground motion and site characteristics, six characteristic parameters are chosen through the feature selection based on 3104 groups of bedrock and surface seismic records collected from the KiK-net strong-motion seismograph network of Japan. Then, the input ground motion characteristics are characterized through the peak bedrock acceleration and predominant frequency, and the site characteristics are characterized by the soil depth at shear wave velocity of 800 m/s, site fundamental period, bedrock shear wave velocity and surface shear wave velocity. The XGBoost model in machine learning is used to establish the prediction models for the PGA based on the above six characteristics. It is shown that the prediction results of the XGBoost prediction model are stable and can be used to predict the PGA better by comparing the records and one-dimensional numerical simulation methods. The coefficients of determination of the training set and the test set are greater than 0.925, and the mean absolute percentage errors are about 20%, which is obviously better than the one-dimensional numerical simulation methods. At the same time, the SHAP is introduced to analyze the influence and dependence between the input characteristics and the predicted results, which enhances the interpretability of the model and provides reliability support for the predicted results.
-
Keywords:
- machine learning /
- XGBoost /
- prediction of PGA /
- SHAP /
- interpretability
-
0. 引言
地表峰值加速度PGA可以直观反映地震动强度，是抗震设计规范中的一个重要参数，同时也是工程结构抗震设计中地震作用特性的直观表征参数，在地震预警、地震烈度速报和地震动区划等领域起到重要作用[1-2]，准确获取和预测PGA具有重要意义。PGA的获取可通过强震观测、地震动预测方程及震害调查等途径获取[3-4]，此外，还可以通过数值模拟方法对PGA进行预测。
目前对PGA的预测中，应用最多的是数值模拟方法，其中二维和三维分析方法更接近于实际情况，但使用三维分析所需的计算费用较高，并且难以确定三维非线性本构关系，相对来说，一维分析方法本构关系简单且能在较低的计算成本的基础上考虑地震动中的相对高频成分，所以现行的土层地震分析方法为一维数值模拟方法[5]。频域等效线性化方法和时域非线性方法是较常用的一维数值分析方法，根据这些方法编制了土层地震反应分析的计算程序，如SHAKE 2000[6]、DEEPSOIL[7]、LSSRLI-1[8]、SOILQUAKE[9]等。虽然SHAKE2000和DEEPSOIL在国际上应用广泛，但其在软土场地上的计算结果并不理想。LSSRLI-1与SHAKE2000原理基本相同，也存在类似问题。SOILQUAKE软件虽然解决了软土场地上计算结果小于实测值的问题，但其在高频部分会出现异常放大的问题[10]。除此之外，还可通过对地震动记录的统计分析进行PGA的预测，如唐川等[11]通过建立场地剪切波速和覆盖层厚度的线性组合，给出了多概率水平PGA预测的概率模型，但该方法中给出的场地特征参数的线性组合并无实际的物理意义，因此该预测模型具有一定的随机性。
随着计算机技术的不断发展，机器学习在地震工程领域也得到了广泛的应用。在地震预警方面，Böse等[12]利用机器学习构建了可持续预测PGA和PGV的神经网络模型；余聪等[2]构建了支持向量机的地震动峰值预测模型，采用P波触发后3 s数据，预测后续地震动强度，改善了常规Pd预测模型中的“小值高估、大值低估”的现象，可以很好的用于现地地震预警。在PGA预测方面，Kerh等[13]使用震中距、震源深度和震级等作为输入特征参数，建立神经网络模型以预测台湾高铁系统沿线10个火车站的PGA，并识别出了潜在的危险站点。此外，还有一些学者建立的机器学习模型[14-15]也能很好地预测PGA。但这些模型多采用震级、震源深度等地震动参数作为模型的输入特征参数，对表征场地特性的参数使用较少。
虽然机器学习有优秀的预测性能，但由于模型的复杂性，无法深入了解输入特征与输出结果之间的关系，从而缺乏对预测结果的解释性，这种模型也被称为黑箱模型。目前，一些方法已经应用到对黑箱模型的解释中，如LIME(local interpretable model agnostic explanation)方法[16]、SHAP(shapley additive explanations)方法[17]等。其中，LIME是一种局部解释方法，只能用来评估特征对单个预测的贡献，不能对模型进行全局解释。SHAP方法则没有这种局限性，它既能从局部对单个样本进行解释，也能对全部样本进行全局解释，SHAP方法由于这种独特优势，已成功应用于多个科学领域，如生物医学工程[18]、生态工程[19]、交通安全[20]、混凝土强度预测[21]、钢筋混凝土失效模型影响因素分析[22]等，但在地震工程中应用较少。
本文通过SHAP方法对表征输入地震动特性和场地特性的不同参数进行重要性排序，筛选出具有较强代表性的6个重要参数，以此作为输入特征参数，利用机器学习中的XGBoost建立PGA预测模型，对从日本KiK-net强震台网中收集到的地震动数据进行训练和测试，同时利用SHAP对预测模型进行解释，并分析各特征参数对预测结果的影响。
1. 数据集的建立
本文以日本KiK-net强震台网[23]的地震动记录为基础，依据以下原则筛选：
(1) 记录时间为1996年3月1日至2020年11月1日；
(2) 强震仪为布设在地表的水平场地台站。本文主要以一维水平场地为研究对象，目标是预测水平场地对地震动峰值的放大效应，地形的放大效应不在本文的研究范围；
(3) 地表峰值加速度大于10 gal。
基于上述原则共筛选出40个KiK-net台站，共计3104条井上井下地震记录，根据中国《建筑抗震设计规范：GB50011—2010》[1]中场地类别的划分标准，将所选场地划分成Ⅱ类、Ⅲ类和Ⅳ类，各场地的地震动记录分布情况如表 1所示，PGA的分布见图 1。
表 1 各类场地台站及地震动记录数量Table 1. Numbers of stations and records of ground motion at various sites场地类别 台站数量/个 地震记录/条 Ⅱ类 6 1524 Ⅲ类 26 1296 Ⅳ类 8 284 合计 40 3104 2. 基本方法和原理
2.1 XGBoost回归模型
XGBoost是由Chen等[24]设计，致力于让提升树突破自身的计算极限，以实现运算快速、性能优越的工程目标的一种集成学习算法。它通过应用梯度下降框架来提升弱学习器，是GB算法的高效实现，与其他GB算法相比，XGBoost能在更短的时间内实现更优的预测性能[25]。
对给定训练集N，从中选取一个样本(xi, yi)，其中xi∈Rn, yi∈Rn, i=1, 2, …, n，用K个可加性函数对样本进行预测，预测值可表示为
ˆyi=n∑i=1fk(xi)。 (1) 式中：xi为第i个样本的n维特征向量；yi为第i个样本的真实值，ŷi为第i个样本的预测值；f(x)为将特征映射到树结构叶子节点权重的函数，可表示为f(x)=ωq(x)，ω为叶子节点的权重，q为一个树结构。
XGBoost回归模型的目标函数由损失函数和正则化项组成，保证了模型的精度和泛化能力。目标函数公式如下：
Obj=n∑i=1l(yi,ˆyi)+K∑k=1Ω(fk)。 (2) 式中：Obj为目标函数；l(yi,ˆyi)为损失函数；Ω(f)为正则化项。模型通过迭代的方式进行训练，迭代的方式为
ˆyi(t)=ˆyi(t−1)+ft(xi)。 (3) 则目标函数的迭代形式为
Obj(t)=n∑i=1l(yi,ˆyi(t−1)+ft(xi))+Ω(ft)。 (4) 对损失函数进行二阶泰勒展开，则目标函数可写成以下形式：
Obj(t)=n∑i=1[l(yi,ˆyi(t−1)+gift(xi))+12hift2(xi)]+Ω(ft) 。 (5) 其中：
gi=∂ˆyt−1il(yi,ˆy(t−1)i)， (6) hi=∂2ˆyt−1il(yi,ˆy(t−1)i)。 (7) 由此可见，XGBoost的目标函数主要依赖于损失函数的一阶和二阶偏导数，这也是该方法的一大优势，它可以不选定具体的损失函数形式，仅依靠输入数据的值进行叶子分裂优化计算，使损失函数的选取与模型算法的优化独立进行，增强了XGBoost的适用性。
2.2 模型评价
本文选择平均绝对误差(mean absolute error，MAE)、均方根误差(root mean squared error，RMSE)、平均绝对百分比误差(mean absolute percentage error，MAPE)和决定系数(R2)做为预测效果的评价指标。其中平均绝对误差是预测值与真实值之间误差绝对值的平均值，反映了预测误差的实际情况；平均绝对百分比误差是预测值与真实值之间的相对误差绝对值的平均值，反映了预测误差与真实值的相对大小；均方根误差是预测值与真实值之间误差平方和与样本数量比值的平方根，反映了误差分布的离散程度；定义预测值与真实值差值的平方和为残差平方和，真实值与其均值之间差值的平方和为总平方和，则R2是1减残差平方和与总平方和比值的差值，反映了模型拟合的优劣程度，R2越接近于1，模型拟合效果越好。其表达式分别为
MAE=1nn∑i=1|yi−ˆyi|， (8) RMSE=√1nn∑i=1(yi−ˆyi)2， (9) MAPE=1nn∑i=1|yi−ˆyiyi|， (10) R2=1−n∑i=1(ˆyi−yi)2n∑i=1(ˉyi−yi)2。 (11) 式中：n为样本量；yi为真实值；ŷi为预测值；ӯ为真实值的平均值。
2.3 模型解释
在模型预测中，了解模型做出某种预测的原因和预测的准确性同样重要。对于简单的模型(如线性模型、单一决策树模型等)，模型本身是可解释性的；对于集成学习算法模型或深度学习模型，该类模型有更好的预测性能，但它们的复杂度也更高，降低了模型的可解释性。为解决XGBoost模型的可解释性问题，本文引入SHAP将预测结果可视化，对模型预测结果进行解释的同时还为模型结果的可靠性提供支撑。
SHAP是由Lundberg等[17]提出的用于解释模型输出的一种解释框架。它基于联盟博弈论[26]的方法，通过计算每个特征对预测结果的贡献值来衡量特征对结果的影响，这个贡献值可能是正值，也可能是负值，其中正值提升了预测结果，负值降低了预测结果。同时，特征在模型中的作用越大，贡献值的绝对值越大，特征重要性越高。
设第i个样本为xi，第i个样本的第j个特征为xij，模型对该样本的预测值为yi，整个模型的基线(模型没有任何输入特征时的预测结果均值)为ybase，f(xij)为第i个样本第j个特征的SHAP值，则SHAP值服从下式：
yi=ybase+f(xi1)+f(xi2)+⋯+f(xij)。 (12) f(xij)的计算公式如下：
f(xij)=∑S⊆F∖{j}|S|!(|F|−|S|−1)!|F|![v(S∪{j})−v(S)]。 (13) 式中：F为样本xi所有影响因素的全集；S为样本xi中任意多个影响因素形成的子集；v(S)为子集S中所包括的影响因素共同作用所产生的贡献；v(S∪{j})-v(S)为影响因素j为该共同作用带来的贡献。
对所有样本的影响因素j累加求均值即为影响因素j的SHAP值f(xj)，公式如下：
f(xj)=M∑i=1f(xij)。 (14) 3. 特征工程
特征工程能从原始数据特征中筛选出更好的特征，以提升模型的预测效果，是机器学习中不可或缺的一部分，具体可以分为特征提取、特征构造和特征选择。
3.1 原始特征
影响场地放大效应的主要因素包含两方面，一方面是输入地震动的特性，如幅值、频率及持时，另一方面则为场地条件，包括场地的剪切波速、覆盖层厚度等。KiK-net台网除提供了台站地表和井下的地震动记录信息外，还提供了对应台站的土层信息，包括土层波速剖面和土层土性概述。为尽量体现场地土层结构的复杂性及输入地震动的特性，本文初步选取以下11个影响PGA的因素来进行分析：
(1) 输入地震动峰值(PBA)，单位gal。
(2) 输入地震动卓越频率(FP)，即输入地震加速度时程傅里叶谱对应的卓越频率，单位Hz。
(3) 场地基本周期(SFP)，单位s，计算公式为
T=√N∑i=1(4hivi)2(2Hihi)。 (15) 式中：T为场地基本周期，hi；Hi和vi分别为第i土层厚度、层中点深度和剪切波速；N为计算土层总数[27]。
(4) 场地覆盖层厚度(OBT)，按《建筑抗震设计规范》场地类别划分中覆盖层厚度计算方法确定，单位为m。
(5) 剖面深度(Depth)，即井下记录仪的埋置深度，单位m。
(6) 剪切波速达800 m/s的土层埋深(D800)，单位m，对于土层剪切波速均小于800 m/s的剖面，取其地震动输入面埋深为D800。
(7) 地表剪切波速(Surface vs)，即场地最上层剪切波速，单位m/s。
(8) 基岩剪切波速(Bedrock vs)，即台站井下记录所在处地层剪切波速，实际工程中为地震动输入界面剪切波速，单位m/s。
(9) 等效剪切波速(vse)，按《建筑抗震设计规范》场地类别划分中等效剪切波速计算公式确定，单位m/s。
(10) 场地地表以下20 m深度范围内的等效剪波速(vs20)，对于钻孔和波速测量深度不足20 m时，可按Boore等[28]建立的速度外推公式进行计算。本文所用场地剖面均大于30 m，可直接计算(m/s)。
(11) 场地地表以下30 m深度范围内的等效剪切波速(vs30)(m/s)。
以上11个影响因素包括反映输入地震动峰值和频谱特性的PBA和FP，另外9个则均是反映场地特性的因素。由于场地结构复杂，可用于描述场地特性的影响因素较多，如对于场地等效剪切波速的描述，既有国内场地类别划分中常用的等效剪切波速vse，也有20 m等效剪切波速，又有国外30 m等效剪切波速，哪一个因素对场地放大效应影响较大，目前并没有较权威的结论，也未达成共识。因而在进行特征选择时，应尽可能包括更多的影响因素，并从中筛选出有较强影响的参数。同时，以上11个影响因素中，各因素也并不相互独立，但利用SHAP方法进行参数重要性筛选时，其算法可有效规避这一限制。
3.2 特征选择
本文在进行特征工程时，仅采用了特征选择。通过特征选择可删除对预测结果没有帮助的特征，以此减少模型特征个数，达到降低模型复杂度、减少运行时间和提高模型预测准确性的目的。将所选特征作为输入参数输入到XGBoost模型中，采用SHAP值对特征进行重要性排序，如图 2(a)所示。由图 2(a)可见，PBA为最重要的影响因素，其SHAP平均值为37.41，D800为次重要影响因素，其SHAP平均值为13.33。对于vs20和vs30，其SHAP值小于1，表明其对模型的贡献较小，可忽略不计。本文选取重要性排序靠前的6个影响因素作为模型的输入特征，分别为PBA、D800、Bedrock vs, FP, Surface vs和SFP。这6个影响因素同时包括了地震动特性和场地条件。
将筛选后的6个特征重新作为输入参数输入到XGBoost模型中，重新计算各个特征的SHAP值，结果如图 2(b)所示。筛选特征后各个特征的重要性排序及SHAP值稍有变化，但幅度不大。且特征筛选前训练集上模型预测结果的决定系数为0.953，特征筛选后训练集上决定系数为0.958，说明模型的性能稍有提高，因此，6参数的模型可替代原11个参数模型，完成本次特征选择。
3.3 模型参数
XGBoost回归模型的重要参数包括：n_estimators参数，用于指定该模型中弱学习器的数量；learning_ rate参数，用于减少每一步的步长，防止因步长太大而跨过了极值点；subsample参数，用于指定从样本中进行采样的比例；booster参数，用于指定要使用的弱学习器类型；max_depth参数，用于指定弱学习器的最大树深度；reg_alpha参数，用于目标函数使用L1正则化并控制正则化强度；reg_lambda参数，用于目标函数使用L2正则化并控制正则化强度。采用十折交叉验证和网格搜索，得到该回归模型的最佳参数如表 2所示。
表 2 XGBoost回归模型最佳超参数Table 2. Best hyperparameters of XGBoost regression model参数 取值 n_estimators 500 learning_rate 0.42 subsample 0.6 booster gbtree max_depth 2 reg_alpha 0 reg_lambda 16 4. 模型预测结果
4.1 整体预测结果
将所选数据的80%作为训练集进行模型训练，剩余的20%作为测试集，进行模型性能测试。模型预测结果如表 3所示。由表 3可见，训练集和测试集的平均绝对误差MAE均较小，在10 gal左右；均方根误差RMSE有所增加，预测结果有一定的离散性；平均相对误差MAPE值，均在20%左右，且R2值均接近于1。从总体上来看，本文所建立的XGBoost模型在训练集和测试集上的预测结果均为理想，可用于一般场地地表峰值加速度的预测。
表 3 XGBoost模型预测结果Table 3. Predicted results by XGBoost model评价指标 MAE RMSE MAPE R2 训练集 10.3 14.55 0.18 0.958 测试集 13.7 20.06 0.20 0.925 为进一步分析预测结果的误差分布情况，将训练集和测试集的预测结果残差(即真实值和预测值的差值)示于图 3。
图 3中黑色实线为平均残差，显示了预测PGA误差的平均变化趋势，黑色虚线为平均残差加减一倍标准差，显示了预测PGA的离散程度。对于训练集和测试集，除极个别工况外，模型预测残差均在正负50 gal内，其平均残差变化趋势一致，说明模型预测结果具有较好的稳定性。当PGA小于40 gal时，平均残差为负值，说明在该区间内模型有轻微的高估；当PGA在40~100 gal时，平均残差基本为0；当PGA在100~500 gal时，平均残差为正值，说明在该区间内模型低估了真实值；当PGA大于500 gal时，训练集平均残差有所降低，测试集的3个样本残差基本为0，说明模型对较大的PGA进行了较为准确的预测。
4.2 与数值模拟方法的对比
将本文的机器学习方法预测结果与数值模拟结果进行对比，通过分析预测结果的准确性，进一步对所给出的算法模型进行评价。数值模拟选择了DEEPSOIL和SHAKE这两个广泛应用的软件，以及笔者给出的改进等效线性化方法SOILRESPONSE，收集了3个程序对测试集所有样本的计算结果[29]。由于台站场地资料中并无土体非线性参数，因而数值模拟时采用了Darendeli[30]的研究成果，按土类和埋深不同共计30组非线性曲线。
将测试集中样本分别按场地类别和PGA大小进行分组，计算并比较不同类别下各种方法的平均绝对百分比误差(MAPE)。对比结果如图 4所示。
由图 4(a)可见，对于3种不同的场地类别，XGBoost模型的MAPE值均在20%左右，平均值为20.4%，预测效果稳定，误差较小。3种数值模拟方法的MAPE均值分别为25.6%，26.2%和28.0%，对比来看，本文所建立的XGBoost预测模型效果更好。
图 4(b)为不同PGA分组下4种方法的预测效果对比，对于PGA小于100 gal的地震动，XGBoost模型预测的MAPE值为20.1%，而3种数值模拟方法的MAPE值均在25%以上；PGA在100~200 gal时，除DEEPSOIL的MAPE值在30%左右，其他3种方法的MAPE值相差不大，在23%左右；PGA大于200 gal时，XGBoost模型的MAPE值仅为11.8%，SOILRESPONSE方法在3种数值模拟方法中预测结果较好，但MAPE值为26.0%，远大于XGBoost模型，DEEPSOIL预测结果最不理想，MAPE值为37.6%。
综上，本文XGBoost预测模型无论在不同场地还是在不同的PGA区间，其预测结果均优于数值模拟方法，尤其对强地震动的预测结果更是明显优于3种数值模拟方法。
5. 预测模型的解释与分析
5.1 整体结构分析
本文所建立的预测模型包含有6个输入特征，图 5是使用SHAP方法计算的特征密度散点图。图 5中，从上至下特征的重要性逐渐减小，散点颜色由蓝至红代表特征的值由小到大，每一个点代表一个样本的SHAP值，该值代表了这个特征对单个预测的贡献，而点的集合表示了特征整体对预测结果影响的方向和大小。
由图 5可见，在这6个特征中，输入地震动峰值PBA是影响模型预测结果的最重要的特征，且PBA对PGA的预测结果呈现出明显的线性影响，PBA的值越高，其SHAP值越高，预测出的PGA值越高，这也和现有认识一致，即输入地震动强度越大，PGA越大。D800、Bedrock vs和输入地震动卓越频率FP是3个次重要的特征，其中D800的值越大，其SHAP值越低，预测的PGA值越低；Bedrock vs的值越小，其SHAP值越低，预测的PGA值越小；FP越大，SHAP值越小，预测的PGA越小。D800越大，说明剪切波速达到800 m/s的土层越深，即场地较近似于深厚软场地，因而放大效应较小。Bedrock vs越小，说明输入界面刚度越小，越接近于土层输入面而非基岩输入面，因而放大效应小。FP越大，即输入地震动卓越频率越大，则其越远离场地的基本周期，难以引起共振效应，因而放大效应也不明显。
对于场地基本周期SFP和地表剪切波速Surface vs，其SHAP值比其他4个特征小，说明这两个特征对PGA的预测结果影响相对较小。总体来说，SFP较小时，SHAP值为正，场地放大效应较明显；Surface vs较小时，对PGA起正向作用。
5.2 特征依赖分析
为深入了解各特征对预测结果影响规律及各特征SHAP值的分布，图 6给出了各个特征的特征依赖图。图 6中，横坐标表示特征值的大小，左轴表示该特征SHAP值的大小，右轴表示与该特征存在最大交互作用的特征值的大小，从蓝到红表示该特征的值由小到大变化。
图 6(a)为PBA的特征依赖图，由图可见，PBA与其SHAP值呈正相关，PBA小于15 gal时，SHAP值为负，即对地表PGA的预测起负向作用，PBA大于15 gal时，SHAP值为正。PBA小于15 gal时，颜色垂直分层不明显；PBA大于15 gal时，红色点主要分布在蓝色下方，说明在该区间时，同一输入PBA下，D800越大，SHAP值越小，预测的PGA值越小，即场地越厚放大效应越不明显。
图 6(b)为D800的特征依赖图，由图可见，D800小于50 m时，SHAP值为正，且PBA越小，SHAP值越小，预测的PGA值越小；D800大于50 m时，SHAP值主要为负值，且总体上呈现PBA越大，SHAP值越小，预测的PGA值越小的趋势，即深厚场地在较大地震动输入时放大效应有减小的趋势。
由图 6(c)所示的Bedrock vs特征依赖图可知，Bedrock vs小于600 m/s时，SHAP值为负，且PBA越大，SHAP值越小，预测的PGA值越小；Bedrock vs大于600 m/s时，SHAP值为正，但除个别较大的PBA样本点，其SHAP值主要集中在0附近，对预测结果的影响不大。
图 6(d)为FP的特征依赖图，由图可见，FP和PBA不存在明显的交互作用，当FP小于7 Hz时，SHAP值主要集中在0附近，说明较小的FP值对PGA的预测结果影响不大；当FP大于7 Hz时，SHAP值为负，且随着FP的增大，其SHAP值呈下降趋势，说明FP越大，预测的PGA值越小。这种变化趋势与现有经验也是一致的，当输入地震动卓越频率较大时，与场地的卓越频率相差较大，避免共振现象的发生，降低了场地放大效应。
图 6(e)为SFP的特征依赖图，当场地周期值小于1 s时，其SHAP值为正，提高了PGA的预测结果；当周期值大于1 s时，SHAP值为负，且PBA越大，SHAP值越小，预测的PGA值越小，即软弱场地对较大地震动的放大起负向作用。
图 6(f)为Surface vs的特征依赖图，由图可见，Surface vs与PBA不存在明显的交互作用，当Surface vs小于110 m/s时，SHAP值主要为正值，提高了PGA的预测结果；当Surface vs大于110 m/s时，SHAP值主要为负值，降低了PGA的预测结果。
5.3 算例分析
为进一步解释各个特征对预测结果的影响，分别绘制不同场地样本的预测瀑布图，如图 7所示，左轴为特征名称和特征值，红色表示对预测结果有正向作用，即提高了预测结果，蓝色表示对预测结果有负向作用，即降低了预测结果。E[f(χ)]为模型基准值，从基准值开始，依次向上加上(红色)或减去(蓝色)各特征的SHAP值，得到该样本的预测结果f(χ)。图 7(a)为Ⅱ类场地FKSH10台站1104111716NS地震动预测结果瀑布图，该样本实际PGA为322.557 gal，本文模型预测PGA为322.674 gal；图 7(b)为Ⅲ类场地TCGH16台站1103111446EW地震动预测结果瀑布图，该样本实际PGA为1195.893 gal，本文模型预测PGA为1130.024 gal；图 7(c)为Ⅳ类场地IKRH02台站1809060308EW地震动预测结果瀑布图，该样本实际PGA为294.558 gal，本文模型预测PGA为306.267 gal。
对比图 7可知，3个输入地震动峰值PBA均大于50 gal，为较强地震动输入，因而都起到了正向作用。图 7(a)场地的剪切波速达800 m/s的土层埋深D800为16 m，地震动输入面的剪切波速Bedrock vs为870 m/s，属于浅硬场地基岩输入，因此这两个特征均起正向作用；图 7(b)场地D800为112 m，虽然属于深厚场地，但由于输入地震动峰值较大，D800仍起到正向作用，但作用不明显；Bedrock vs为680 m/s，属于基岩输入，此特征起正向作用；图 7(c)场地为深厚场地非基岩面输入，D800和Bedrock vs均起负向作用。图 7(a)场地表层较硬，因而起负向作用，反观图 7(b)和图 6(c)，其地表剪切波速分别为80 m/s和70 m/s，对场地放大则起正向作用。图 7(b)场地的基本周期与输入地震动的卓越频率相差不大，因而SFP和FP起到了较明显的正向作用，而图 7(a)和图 7(c)中二者相差较大，因而作用都不是很明显。
由图 7可知，瀑布图可清楚地解释各预测结果与各影响因素之间的关系，较以往的“黑箱”模型有了较大的进步。需要说明的是，由于瀑布图是单个样本的预测结果，每个样本的特征重要性顺序并不相同，且与总体样本的特征排序(图 5)也会有所不同。
6. 结论
本文利用XGBoost建立了基于基岩地震动和场地主要参数的PGA预测模型，同时采用SHAP方法对模型进行了解释和特征影响分析，为地震区划和工程结构抗震设计提供了一种简单实用、无需建模分析的快速预测地表地震动峰值的方法。主要得到以下5点结论。
(1) 通过特征选择，选取了PBA、D800、Bedrock vs、FP、SFP和Surface vs等6个参数作为机器学习的输入特征。
(2) 建立了基于XGBoost的PGA预测模型，决定系数R2大于0.925，模型的拟合效果较好；测试集预测结果的MAPE为20.3%，误差较小。
(3) 对于3种不同场地类别，XGBoost模型预测效果稳定，MAPE均在20%左右，预测结果均优于数值模拟方法。
(4) 对于不同的PGA区间，XGBoost模型预测效果均优于数值模拟方法。特别是在强地震动下，XGBoost模型的MAPE仅为11.8%，明显优于数值模拟方法，为强地震动下的场地反应预测提供了新方法。
(5) 利用SHAP对模型进行解释，结果表明PBA是影响PGA预测结果的最主要特征，且PBA的SHAP值随PBA的增大而增大；D800和Bedrock vs是次重要的特征，特征值小时，SHAP值较大；FP、SFP和Surface vs的重要性相对较小，SHAP值主要分布在正负20区间内。SHAP值对预测结果的解释与现有认识一致，进一步证明了模型预测结果的可靠性。
-
表 1 各类场地台站及地震动记录数量
Table 1 Numbers of stations and records of ground motion at various sites
场地类别 台站数量/个 地震记录/条 Ⅱ类 6 1524 Ⅲ类 26 1296 Ⅳ类 8 284 合计 40 3104 
下载: 导出CSV
表 2 XGBoost回归模型最佳超参数
Table 2 Best hyperparameters of XGBoost regression model
参数 取值 n_estimators 500 learning_rate 0.42 subsample 0.6 booster gbtree max_depth 2 reg_alpha 0 reg_lambda 16
下载: 导出CSV
表 3 XGBoost模型预测结果
Table 3 Predicted results by XGBoost model
评价指标 MAE RMSE MAPE R2 训练集 10.3 14.55 0.18 0.958 测试集 13.7 20.06 0.20 0.925
下载: 导出CSV
-
[1] 建筑抗震设计规范: GB 50011—2010[S]. 北京: 中国建筑工业出版社, 2010. Code for Seismic Design of Buildings: GB 50011—2010[S]. Beijing: China Architecture & Building Press, 2010. (in Chinese)
[2] 余聪, 宋晋东, 李山有. 基于支持向量机的现地地震预警地震动峰值预测[J]. 振动与冲击, 2021, 40(3): 63-72, 80. https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ202103010.htm YU Cong, SONG Jindong, LI Shanyou. Prediction of peak ground motion for on-site earthquake early warning based on SVM[J]. Journal of Vibration and Shock, 2021, 40(3): 63-72, 80. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ202103010.htm
[3] BOORE D M, ATKINSON G M. Ground-motion prediction equations for the average horizontal component of PGA, PGV, and 5%-damped PSA at spectral periods between 0.01 s and 10.0 S[J]. Earthquake Spectra, 2008, 24(1): 99-138. doi: 10.1193/1.2830434
[4] DU K, DING B R, BAI W, et al. Quantifying uncertainties in ground motion-macroseismic intensity conversion equations. A probabilistic relationship for western China[J]. Journal of Earthquake Engineering, 2020: 1-25.
[5] 张震. 场地地震反应一维数值分析方法对比分析[D]. 廊坊: 防灾科技学院, 2020. ZHANG Zhen. Comparison on one Dimension Numerical Methods of Site Seismic Response Analysis[D]. Langfang: Institute of Disaster Prevention, 2020. (in Chinese)
[6] SCHNABEL P B, LYSMER J, SEED H B. SHAKE: A Computer Program for Earthquake Response Analysis of Horizontal Layer Sites[R]. Berkeley: University of California, 1972.
[7] HASHASH Y M, PARK D. Non-linear one-dimensional seismic ground motion propagation in the Mississippi embayment[J]. Engineering Geology, 2001, 62(1): 185-206.
[8] 李小军. 一维土层地震反应线性化计算程序[M]. 北京: 地震出版社, 1989. LI Xiaojun. A Computer Program for Calculating Earthquake Response of Ground Layered Soil[M]. Beijing: Seismological Press, 1989. (in Chinese)
[9] 袁晓铭, 李瑞山, 孙锐. 新一代土层地震反应分析方法[J]. 土木工程学报, 2016, 49(10): 95-102, 122. https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201610015.htm YUAN Xiaoming, LI Ruishan, SUN Rui. A new generation method for earthquake response analysis of soil layers[J]. China Civil Engineering Journal, 2016, 49(10): 95-102, 122. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TMGC201610015.htm
[10] SUN R, YUAN X M. A holistic equivalent linear method for site response analysis[J]. Soil Dynamics and Earthquake Engineering, 2020, 141: 106476.
[11] 唐川, 陈龙伟. 场地校正的地表PGA放大系数概率模型研究[J]. 工程力学, 2020, 37(12): 99-113. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX202012011.htm TANG Chuan, CHEN Longwei. Probability modelling of pga amplification factors corrected by site conditions[J]. Engineering Mechanics, 2020, 37(12): 99-113. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX202012011.htm
[12] BÖSE M B, HEATON T, HAUKSSON E. Rapid estimation of earthquake source and ground-motion parameters for earthquake early warning using data from a single three- component broadband or strong-motion sensor[J]. Bulletin of the Seismological Society of America, 2012, 102(2): 738-750.
[13] KERH T, TING S B. Neural network estimation of ground peak acceleration at stations along Taiwan high-speed rail system[J]. Engineering Applications of Artificial Intelligence, 2005, 18(7): 857-866.
[14] DERRAS B, BARD P Y, COTTOM F, et al. Adapting the neural network approach to pga prediction: an example based on the kik-net data[J]. Bulletin of the Seismological Society of America, 2012, 102(4): 1446-1461.
[15] DHANYA J, RAGHUKANTH S T G. Ground motion prediction model using artificial neural network[J]. Pure and Applied Geophysics, 2018, 175(3): 1035-1064.
[16] RIBEIRO M T, SINGH S, GUESTRIN C. "Why should I trust you?": explaining the predictions of any classifier[C]// Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. San Francisco, CA, USA, 2016: 1135-1144.
[17] LUNDBERG S M, LEE S I. A Unified approach to interpreting model predictions[C]// Conference and Workshop on Neural Information Processing Systems. California: NIPS Press, 2017: 4765-4774.
[18] LUNDBERG S M, NAIR B, VAVILALA M S, et al. Explainable machine-learning predictions for the prevention of hypoxaemia during surgery[J]. Nature Biomedical Engineering, 2018, 2(10): 749-760.
[19] CHA Y, SHIN J, GO B, et al. An interpretable machine learning method for supporting ecosystem management: Application to species distribution models of freshwater macroinvertebrates[J]. Journal of Environmental Management, 2021, 291: 1-13.
[20] PARSA A B, MOVAHEDI A, TAGHIPOUR H, et al. Toward safer highways, application of XGBoost and SHAP for real-time accident detection and feature analysis[J]. Accident Analysis and Prevention, 2020, 136(C): 105405.
[21] LYNGDOH GIDEON A, MOHD Z, ANOOP K N, et al. Prediction of concrete strengths enabled by missing data imputation and interpretable machine learning[J]. Cement and Concrete Composites, 2022, 128: 104414.
[22] MANGALATHU S, HWANG S H, JEON J S. Failure mode and effects analysis of RC members based on machine-learning-based SHapley Additive exPlanations (SHAP) approach[J]. Engineering Structures, 2020: 1-10.
[23] National Research Institute for Earth Science and Disaster Resilience(NIED) Strong-Motion Seismograph Networks (KiK-net)[OL]. https://www.kyoshin.bosai.go.jp/.
[24] CHEN T, HE T. Higgs boson discovery with boosted trees[J]. JMLR: Workshop and Conference Proceedings, 2015, 42: 69-80.
[25] CHEN T, GUESTRIN C. XGBoost: A scalable tree boosting system[C]// Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. San Francisco, CA, USA, 2016: 785-794.
[26] SHAPLEY L S. A value for n-person games[J]. Contributions to the Theory of Games, 1953, 2(28): 307-317.
[27] 齐文浩, 薄景山, 刘红帅. 水平成层场地基本周期的估算公式[J]. 岩土工程学报, 2013, 35(4): 779-784. http://www.cgejournal.com/cn/article/id/15016 QI Wenhao, BO Jingshan, LIU Hongshuai. Fundamental period formula for horizontal layered soil profiles[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(4): 779-784. (in Chinese) http://www.cgejournal.com/cn/article/id/15016
[28] BOORE D M, THOMPSON E M, CADET H. Regional correlations of VS30 and velocities averaged over depths less than and greater than 30 meters[J]. Bulletin of the Seismological Society of America, 2011, 101(6): 3046-3059.
[29] 孙锐, 袁晓铭. 全局等效线性化土层地震反应分析方法[J]. 岩土工程学报, 2021, 43(4): 603-612. doi: 10.11779/CJGE202104002 SUN Rui, YUAN Xiaoming. Holistic equivalent linearization approach for seismic response analysis of soil layers[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(4): 603-612. (in Chinese) doi: 10.11779/CJGE202104002
[30] DARENDELI M B. Development of A New Family of Normalized Modulus Reduction and Material Damping Curves[D]. Austin: The University of Texas at Austin, 2001.
-
期刊类型引用(5)
1. 李威，熊凌，罗钟邱，吴经纬，万诗斐，但斌斌. 基于加权聚类和DNN的KR法脱硫剂加入量预报模型. 炼钢. 2025(01): 12-18+44 . 
百度学术
2. 谌柳谦. 紧邻深基坑的历史保护建筑保护措施关键技术研究. 建筑施工. 2024(01): 81-84 . 百度学术
3. 彭白雪，陈清华，季家东. 基于XGBoost和SHAP的制冷系统故障分析. 低温与超导. 2024(07): 89-96 . 百度学术
4. 曹放，孙徐，张钰. 基于“XGBoost—SHAP”的可解释性崩塌落石风险预测在公路工程中的应用. 工程技术研究. 2024(14): 1-4 . 百度学术
5. 龙潇，孙锐，郑桐. 基于卷积神经网络的液化预测模型及可解释性分析. 岩土力学. 2024(09): 2741-2753 . 百度学术
其他类型引用(6)
-
其他相关附件
计量
- 文章访问数: 364
- HTML全文浏览量: 50
- PDF下载量: 103
- 被引次数: 11
