智能建造是水利工程的重点领域，对于一些特殊的、复杂的地质环境，重大水利工程智能建造十分需要智能化的强震数据监测。在工程前期设计阶段的抗震安全评价方面，区域的强震数据监测，可以为重大水利工程决策提供依据。由于强震作用时，地震波通过峡谷地形在超深覆盖层中经多次反射传播，加上深覆盖层的阻尼耗散，既有各种放大效应，又有能量耗散。当前，对于超深覆盖层地震波放大效应存在较大分歧，科学合理地确定深埋土体动静力参数以及坝址区强震动荷载输入仍面临技术挑战。这就需要智能化的强震数据监测来得到实测地震记录及超深覆盖层地震放大效应，从而解决设计建造的技术难题。

在地震监测及水利工程设计中，需要大量的工程场地实测地震数据，对场地和结构的响应进行大量的非线性计算，从而研究地震作用下结构响应的规律。然而目前的强震记录观测历史较为短暂，观测受设备限制，缺乏符合设计标准的实测地震波数据，这成为了重大水利工程开发中抗震安全评价的主要技术瓶颈。

另外，在地震研究中，对地震波本身的认识十分匮乏。首先对地震波的认识局限于地震波对于结构的作用，这种对结构响应的研究又局限于以峰值响应、总能量为代表的地震动参数刻画，对于地震波时域、频域特征的认识十分不足。另外，地震作用既有规律性，又有很强的不确定性，因而也存在对其天然不确定性上的认识不足。要在地震波的时域、频域特征上加深认识，需要更加完善、更加真实的人工模拟地震波作为资料支撑，来弥补历史实测强震资料不足的问题。地震动的衰减关系和地震波反应谱的形状与震源机制相关^[1]，而目前对震源机制也缺乏深入的认识。地震波时程是一种复杂的、瞬态的、非平稳的时间序列，其中“非平稳性”(nonstationarity)指地震波时程的幅值和频率随时间而变化，这种幅值和频率的非平稳性能极大地影响结构的非线性地震反应^[2-3]。

然而，地震波这种非常复杂的时序信号，人们现在对它的认识不足，常规的方法对地震波的解释不够。而近年来人工智能方法快速发展，工程问题中逐渐引入具有良好预测能力和泛化能力的机器学习等算法，不仅能够处理复杂的数学、物理计算问题，还能在此基础上发现其中的特性与规律。在地震研究中，Ursino等^[4]将机器学习方法应用于地震事件的鉴别和分类问题；Wiszniowski等^[5]通过机器学习方法来提高地震探测能力；Perol等^[6]进行了波形分类以及对地震的检测和定位工作。XGBoost方法相比于其他机器学习方法，引入了正则化思想，来降低树的复杂度，从而获得更好的模型性能。

本文基于人工智能XGBoost算法，通过对工程场地实测地震波和人工模拟地震波进行判别，进而分析研究地震波的时域和频域特征，从而加深对地震波本身的认识，为进一步研究、模拟地震波提供技术支撑。

1 地震数据库与人工地震波仿真方法 1.1 实测地震波数据库

本文选用的实测地震波数据来源于美国太平洋地震中心下一代强震数据库^[7](PEER-NGA-West2 Database)，NGA-West2项目数据库在一代数据库的基础上进行了扩展。在NGA-West2数据库中，包括了2000年后大量的活跃构造区浅层地壳地震记录的全球地面运动数据，以及1998—2011年加利福尼亚州的部分中小地震。

NGA-West2数据库收录了599次地震动事件的21 336条3分量记录。其中，主要包括中国、美国加州、日本、地中海地区等国家或地区发生的599次地震。数据库覆盖的参数空间为震级3.0~7.9级，震中距离为1~300 km，地表土层30 m的平均剪切波速V_s30范围为94~2 100 m/s(但是震中距离超过400 km、V_s30超过1 200 m/s或小于150 m/s时，数据分布变得稀疏)，图 1为NGA-West2数据库地震动记录的场地、区域和震级-震中距分布情况。

图 1c为数据的震级与震中距关系分布。该数据库收录了在11种阻尼比下，0.01~20 s的111个周期的统一处时间序列和响应谱坐标。震源、路径和场地条件等地震动震源机制和原始数据由开发人员和专业工作小组进行质量检查。

首先在NGA-West2数据库中选取实测地震波数据，通过观察波形，计算峰值加速度PGA，最终得到6 967条水平方向，3 489条竖直方向，共10 456条NGA-West2数据库实测地震波数据。所选数据的震级范围为3.0~7.9级，震中距范围为1~300 km。由于本文希望从地震波时程波形本身入手，对数据进行分辨，故对地震波数据的震源、场地条件等参数不作过多限制。

1.2 小波包分析

小波包(wavelet packet)分析技术是近年来在小波分析(wavelet)的基础上发展起来的先进的信号处理方法。对比小波分析，小波包对信号的分解、重构具有更高的时频分辨率，是一种更精细的信号分析方法。小波包的母函数是一种在某个时域段和频域段上分布的局部函数。运用Meyer正交小波包为母函数，对地震波时间序列进行多级分解，得到与母函数相对应的小波包系数，而小波包系数的绝对值在时域-频域的分布，被称为“小波包谱”。大量的地震波分解实例表明，小波包系数的正负号可以近似为随机分布的。因此，利用小波包谱可以对地震波进行有效重构。用小波包技术合成人工地震波是近年来由美国斯坦福大学Baker研究团队首先提出的^[8]。

图 2为一个小波包分解实例，数据来源为1994年美国加州北岭(Northridge)地震，通过对加速度时程实测数据进行分解得到，小波包谱的时频分布很清晰地表现出地震波的非平稳特性：该地震波前5 s的小波包系数主要集中在1~10 Hz的频率分布范围；在20 s后，地震波频率主要分布在0.2~3 Hz。而传统的Fourier谱(FAS)只能表征地震波在频域上的强度分布，不能对地震波的频率随时间变化的特性进行刻画。

Baker团队对美国太平洋地震中心下一代强震数据库(PEER-NGA)中的1 408个地震波进行了小波包时频域特征的分解和统计分析，总结出13个小波包参数(表 1)，描述小波包谱在时频域的相关统计特征，并建立了小波包参数的衰减关系模型。其中，E_acc和E(a)_major表征地震波能量，E(t)_major、E(f)_major、E(t)_minor和E(f)_minor表示主要和次要能量组小波包分布在时间轴和频率轴的形心，S(t)_major、S(f)_major、S(t)_minor和S(f)_minor表示主要和次要能量组小波包分布的标准差，ρ(t, f)_major和ρ(t, f)_minor表征时频域分布相关性，S(ξ)为随机参数。清华大学黄杜若团队^[9]首先发展了小波包时频域特征参数之间的空间相关性，并提出了预测模型，首次生成了在空间上分布的人工地震波。进一步的研究结果还表明，小波包参数的空间相关度与区域场地条件相关^[10]。在均匀地质场地，小波包参数的空间相关性可以达到30 km，比不均匀地质场地的相关性大。

1.3 SIMQKE人工地震波模拟方法

由于缺乏地震记录，以及广泛使用地震动时程的动力学分析来获得结构和二次系统的响应，就需要发展人工模拟地震波方法。现代的地震工程时程分析中使用的人工地震波，主要通过地震波选择与改造法、物理模拟法和随机分析法3类方法生成^[11]。其中，地震波选择与改造法选择合适的地震波实测数据，对它的幅值进行缩放，从而使其与目标反应谱相匹配^[12]，是目前最为广泛使用的方法。目前收集的地震记录中的一个显著缺陷是，从未记录过真正的大地震(里氏8级)中发震断层附近的震动^[13]。SIMQKE方法的目的是关注一种常用的数值模拟方法，该方法的基础是任何周期函数都可以扩展为一系列正弦波。

SIMQKE仿真方法主要是从指定的平滑反应谱计算功率谱密度函数；在与指定反应谱相匹配的时程上生成统计独立的人工模拟加速度时程，并通过迭代过程细化谱匹配。它还对生成的运动执行基线校准，以确保最终地面速度为零，以时程作为输入计算反应谱，其流程如图 3所示。

在模拟过程中，如果计算的反应谱在目标上方和下方变化过大，可使用程序中的循环选项来平滑匹配。在每个循环中，将计算出的实际响应与一组控制频率下的目标进行比较。在每个控制频率下获得期望响应与计算响应的比率，并且功率谱密度的相应值与该比率的平方成比例地进行修改。

进一步，在筛选NGA实测地震波数据后，选用SIMQKE人工模拟地震波方法，将筛选出的实测NGA数据作为输入数据，经过3次迭代，模拟生成对应NGA实测地震波的人工地震波，并对生成的地震波进一步筛选。最终得到2 956条水平方向，1 475条竖直方向，共4 431条SIMQKE人工模拟地震波(图 4)。对人工模拟生成的地震波进行抽查，保证数据质量，如图 5所示为SIMQKE人工模拟地震波方法生成的地震波时程、反应谱及功率密度谱对比，数据来源为1984年意大利5.8级地震，可见SIMQKE方法生成的地震波数据质量良好。

2 XGBoost算法判别实测和人工地震波

XGBoost算法是梯度增强回归树的一种改进算法，它引入了正则化思想，来降低树的复杂度，从而获得更好的模型性能。回归树(regression tree，RT)是梯度增强回归树(gradient boosted regression tree，GBRT)的基本模型。本研究使用小波包参数X₁和X₂作为预测变量的单个回归树模型，用前述的小波包谱分解参数替代，预测器变量空间被划分为多个区域，每个区域由从根部分割节点到对应的叶子节点的路径表示。

如图 6所示，回归树的每个叶子节点都有一个特定的叶得分(w_i)，它代表该区域的预测值结果，并将其拟合为该区域目标变量样本的平均值。回归树的训练是通过遍历搜索最优分割节点，并继续进行分割，直到达到停止条件为止的过程，也就是树的生长过程。

集成建模的思想是利用一组较弱的基础模型来构造一个强大的模型，并行集成技术被广泛应用于降低模型的不确定性。梯度增强回归树(GBRT)作为一种序列集成技术，随着时间的推移发展出一系列基本的回归树(RT)，以扩大模型容量。提升模型的训练过程以正向逐步迭代的方式，可以表示为

$\begin{array}{c} {{\hat y}^{(0)}} = 0, \\ {{\hat y}^{(1)}} = \nu {f_1}\left( {x;{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_1}} \right), \\ {{\hat y}^{(2)}} = \nu \sum\limits_{j = 1}^2 {{f_j}} \left( {x;{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right) = {{\hat y}^{(1)}} + \nu {f_2}\left( {x;{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_2}} \right), \\ \cdots \\ {{\hat y}^{(T)}} = \nu \sum\limits_{j = 1}^T {{f_j}} \left( {x;{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right) = {{\hat y}^{(T - 1)}} + \nu {f_T}\left( {x;{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_T}} \right). \end{array}$

(1)

其中，T为用于迭代提升的回归树数量；Θ_j为第j个回归树的结构(包括树的所有分割节点和叶子节点)；ν为收缩因子(也称为满足0＜ν＜1的学习率，用于收缩单个回归树的贡献)；$ {{\hat y}^{(j)}}$为使用第j个回归树预测的目标变量；f_j()是没有收缩的第j个回归树的输出，它使用预测变量x来近似$ y - {{\hat y}^{(j - 1)}}$(即残差)，具有树的结构Θ_j。因此，残差通常会随着回归树数量的增加而减少，梯度增强回归树的整体结构示意图如图 7所示。

为了控制模型的复杂度，发展出了一种可伸缩的树提升系统XGBoost，它在GBRT的传统损失函数中加入了正则化项，遵循正则化学习的一般思想^[15]。XGBoost建模的核心任务是在第j步找到最优的Θ_j并建立f_j (x; Θ_j)，这是通过使目标函数最小化来实现的。

$\begin{array}{l} {\mathit{\boldsymbol{ \boldsymbol{\widehat \varTheta} }}_j} = \arg {\min\limits_{{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}}}\left\{ {\sum\limits_{i = 1}^T L \left[ {{y_i}, \hat y_i^{(j - 1)} + \nu {f_j}\left( {{x_i};{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right)} \right] + } \right.\\ \left. {\mathit{{ {\varOmega} }}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right)} \right\} = \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \{ \sum\limits_{i = 1}^T {\left[ {\hat y_i^{(j - 1)} - } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{{\left. {{y_i} + \nu {f_j}\left( {{x_i};{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right)} \right]}^2} + \mathit{{ {\varOmega} }}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right)} \right\} \end{array}$

(2)

$\mathit{\Omega }\left( {{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right) = \gamma {M_j} + \frac{1}{2}\lambda ||{w_k}|| = \gamma {M_j} + \frac{1}{2}\lambda \sum\limits_{k = 1}^{{M_j}} {w_k^{(j)2}} .$

(3)

式(2)中，$ L(y, \hat y) = {(\hat y - y)^2}$为常用的平方损失函数；Ω(Θ_j)为第j个回归树上的正则化项，正则化项可表示为式(3)，其中，w_k^(j)为第j个回归树中第k个叶节点的叶得分(图 7)；M_j是第j个回归树中叶子节点的数量；γ为回归树进一步划分节点时，所需要的最小损失函数减少值；λ是回归树中叶得分的L2正则化项。通过二阶泰勒展开，式(2)可近似表示为：

$\begin{array}{l} {\mathit{\boldsymbol{ \boldsymbol{\widehat \varTheta} }}_j} \simeq \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \left\{ {\sum\limits_{i = 1}^T {\left[ {L\left( {{y_i}, \hat y_i^{(j - 1)}} \right) + \upsilon g_i^{(j)}{f_j}\left( {{x_i};{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right) + } \right.} } \right.\\ \left. {\left. {\;\;\;\;\;\;\frac{1}{2}{\upsilon ^2}h_i^{(j)}f_j^2\left( {{x_i};{\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_j}} \right)} \right] + \gamma {M_j} + \frac{1}{2}\lambda \sum\limits_{k = 1}^{{M_j}} {w_k^{(j)2}} } \right\} \end{array}$

(4)

$g_i^{(j)} = \frac{{{\rm{ \mathsf{ ∂} }} L\left( {{y_i}, \hat y_i^{(j - 1)}} \right)}}{{{\rm{ \mathsf{ ∂} }} \hat y_i^{(j - 1)}}} = 2\left( {\hat y_i^{(j - 1)} - {y_i}} \right), $

(5)

$h_i^{(j)} = \frac{{{{\rm{ \mathsf{ ∂} }} ^2}L\left( {{y_i}, \hat y_i^{(j - 1)}} \right)}}{{{{\left( {{\rm{ \mathsf{ ∂} }} \hat y_i^{(j - 1)}} \right)}^2}}} = 2.$

(6)

很明显，叶片更多的树(较大的M_j)会受到较大的惩罚，而较大的λ将产生更规则的叶得分分布(w_k^(j)较小)。由于y_i是给定的样品，并且$ {\hat y_i^{(j - 1)}}$已经在第j-1步测定，因此L(y_i, $ {\hat y_i^{(j - 1)}}$)可以被视为一个常数项，并将从式(4)中删除。XGBoost模型在确定所有Θ_j之后完成训练，因此公式(1)可用于执行预测。

前面描述的T、ν、γ和λ是XGBoost的输入参数，与在训练过程中确定的模型参数Θ不同，输入参数应该在训练前指定。除此之外还需要考虑，单个回归树的最大深度d_max(根部分割节点与叶子节点之间距离最远的边数)以及回归树中进一步划分所需的最小实例权重之和w_mc。

3 实测与人工地震波时频特征分析和判别研究

如前所述，目前的人工地震波方法，尽管能够在反应谱、总能量等方面拟合得很好，但是在时域、频域较为全面的特征上仍有很大不足，所以考虑利用前述的小波包技术，求出NGA-West2实测数据和SIMQKE模拟数据的小波包参数，通过研究二者的时频域特性，发现人工模拟地震波究竟在哪些参数上拟合不足。所以对已有的实测和模拟地震波求出其小波包参数，并标注实测或模拟数据，NGA实测真实地震数据标注1，人工地震数据标注0，分组整合成数据文件。图 8给出了本模型的工作流程图。

在建立判别人工和实测地震波的XGBoost模型时，将已有数据打乱并分成两组，其中一组作为训练数据，对XGBoost模型进行训练，另外一组保留进行XGBoost模型的结果验证。在训练过程中，设定一定的训练次数，最终取训练效果最好的作为输出模型。

在对模型进行评价时，将保留的一组数据输入到模型中，模型判断给出输入数据是实测地震波还是模拟地震波，将此结果与输入数据本来的标注对比，得到模型的预测准确率。另外，模型在判断过程中用到的输入参数重要性也是一项重要的评价指标，输入参数的重要性代表了人工模拟地震波和NGA实测地震波在哪些参数上差距比较大。最终，根据小波包参数重要性，得到在判断中较为重要的小波包参数。进一步绘制NGA实测地震波和人工模拟地震波的某个小波包参数随震中距变化的分布图，对比其参数分布进行模型评估。

在求出10 456条NGA-West2数据库实测地震波数据和4 431条SIMQKE方法人工模拟地震波数据的小波包参数后，标注实测或模拟地震数据，并打乱随机分为两组，其中一组作为训练数据，另一组作为验证数据。将训练数据输入XGBoost模型进行若干次训练，取出训练效果最好的一次模型。用此模型分辨剩余的另一组数据，并输出其为实测地震波还是模拟地震波，最后将输出结果与标注的标签对比，得到预测准确率达到了91%。

由XGBoost模型的预测准确率结果可得，对NGA实测地震波与SIMQKE模拟地震波的判断达到了91%，可见SIMQKE人工地震波模拟方法确实与真实地震存在很大差异，进一步研究是哪些参数存在差异。

通过统计XGBoost模型在回归树结构中的小波包参数作为节点的次数，得到模型在判断时不同输入参数次数对判断结果的影响程度，即输入的小波包参数的重要性，如图 9所示。重要性最高的是ρ(t, f)_minor这个参数，它是刻画地震波的时域和频域相关性的一个参数，另外还有E(t)_minor、E(f)_minor和E_acc 3个参数的重要性较高，由此可知，SIMQKE方法虽然在反应谱、能量等特征上和真实地震波匹配得较好，但是在时频域相关性上存在很大差异。

将已有的NGA实测地震波和SIMQKE模拟地震波数据求得的上述4个参数与震中距关系绘图，得到图 10。图中黑色圆圈代表NGA实测地震波数据，红色圆圈代表SIMQKE人工模拟地震波数据，可见二者在ρ(t, f)_minor这个参数上区别十分明显，SIMQKE模拟地震波此小波包参数集中于-0.15~0.05，而NGA实测地震波此参数则分布较为分散。

综上所述，SIMQKE人工地震波模拟方法在反应谱和能量上能够匹配实测数据，而XGBoost方法从时频域特征的小波包参数入手，能够较好地判别SIMQKE人工地震波和实测地震波，可见综合使用小波包技术分析和XGBoost方法进行判断，能够发现人工地震波模拟方法与实测地震数据的区别。

4 结语

本文基于人工智能XGBoost机器学习算法，选用NGA-West2数据库中的实测地震波数据和SIMQKE方法模拟的人工地震波数据，对模拟地震波与实测地震波的时频特征进行分析。利用小波包参数作为回归树的分割节点，迭代提升回归树形成XGBoost的基本结构，建立算法的判别模型，对模拟与实测地震波进行判别，从而对二者的时频特征进行研究。

在筛选NGA-West2数据库并进行SIMQKE模拟后，基于小波包分析方法对已有数据进行时频特征分析，得到小波包参数。使用XGBoost方法，通过对SIMQKE人工模拟方法生成的地震波与NGA-West2数据库中的实测地震波进行判别，对NGA实测地震波与SIMQKE模拟地震波的判断准确率达到91%，说明从地震波小波包分析时频域特征的角度入手，能够较好地分辨出SIMQKE方法模拟的地震波与实测地震波。

进一步对判别实测和人工地震波中小波包参数的重要性进行研究，发现二者差异最大的是用于描述地震波时域、频域相关性的参数。绘制该小波包参数与震中距关系图，发现SIMQKE人工地震波的该参数较为集中，而真实地震波的该参数分布较为均匀。SIMQKE方法在反应谱、能量等特征上和真实地震波匹配得很好，而本研究的XGBoost方法通过分析地震波的时频域特征，能够发现SIMQKE人工地震波与实测地震波的区别。本文提出的智能方法基于XGBoost机器学习算法，从地震波的小波包参数入手，对实测和人工地震波进行分析，对地震波的时域、频域特性的研究分析具有参考意义。在接下来的研究中，将丰富输入参数，强化模型判别能力，如断层类型等震源机制的影响以及场地条件的影响等，从而进一步对地震波的特征进行深入分析与研究。