陳蒙，王華,2*

1 電子科技大學(xué)資源與環(huán)境學(xué)院，成都 611731 2 電子科技大學(xué)(深圳)高等研究院，廣東深圳 518110

0 引言

地震引起的強(qiáng)地面運(yùn)動會造成建筑結(jié)構(gòu)破壞和倒塌，是震后人員傷亡和財產(chǎn)損失的主要原因，同時還會誘發(fā)滑坡等次生災(zāi)害.地震發(fā)生后，快速估算地震動強(qiáng)度參數(shù)(峰值加速度PGA、峰值速度PGV、峰值位移PGD和加速度反應(yīng)譜PSA)得到地震動圖(ShakeMap)(Wald et al., 2005; 李俊等, 2010; 陳鯤等, 2010; 劉軍等, 2016)，可用于評估震后損失，指導(dǎo)應(yīng)急救援和災(zāi)后重建工作.地震發(fā)生前，分析潛在震源區(qū)可能產(chǎn)生地震動強(qiáng)度，可獲得地震動參數(shù)區(qū)劃圖或地震危險性圖, 以作為房屋保險等的重要參考，更有助于國土空間規(guī)劃和重大工程項(xiàng)目(例如核電站、水庫大壩、跨海大橋等)的設(shè)計和建設(shè).

計算(或預(yù)測)地震動強(qiáng)度參數(shù)的方法主要有三種：數(shù)值模擬、地震動預(yù)測方程(Ground Motion Prediction Equation，也稱作地震動衰減關(guān)系)和機(jī)器學(xué)習(xí).常用的數(shù)值模擬方法包括有限差分法(Graves, 1998; Kozdon et al., 2013; Zhang et al., 2014)、有限元法(Bielak et al., 2005; Duan, 2010)、譜元法(Komatitsch et al., 2004; Lee et al., 2009)、有限體積法(Benjemaa et al., 2009)以及這些方法的結(jié)合(Meng and Wang, 2018; Wu et al., 2018).數(shù)值方法可以模擬斷層破裂和波場傳播的全過程，具有明確的物理意義.但是準(zhǔn)確的地震動模擬需要精確的震源信息和地下結(jié)構(gòu)速度模型，并且計算量巨大.地震動預(yù)測方程(Boore et al., 2014; Campbell and Bozorgnia, 2014; Chiou and Youngs, 2014; Idriss, 2014)形式明確、計算速度快，是地震動圖和地震危險性分析中常用的方法.但是現(xiàn)代地震動預(yù)測方程復(fù)雜(Abrahamson et al., 2014; Campbell amd Bozorgnia, 2014)，函數(shù)形式和特征變量的選取沒有統(tǒng)一標(biāo)準(zhǔn)，依賴研究人員的主觀性，且易忽略各特征之間的非線性耦合作用.

隨著人工智能技術(shù)的發(fā)展，用于預(yù)測地震動強(qiáng)度參數(shù)的機(jī)器學(xué)習(xí)方法受到越來越多的重視，相關(guān)方法主要可以分為兩類：

(1)基于純數(shù)據(jù)驅(qū)動演化的預(yù)測模型

(2)基于監(jiān)督學(xué)習(xí)回歸模型的地震動強(qiáng)度參數(shù)預(yù)測

主要是通過多種地震記錄信息，訓(xùn)練得到可用于預(yù)測地震動強(qiáng)度參數(shù)的機(jī)器學(xué)習(xí)模型.例如Arjun和Kumar(2009)利用單隱藏層前饋神經(jīng)網(wǎng)絡(luò)和84456條K-NET記錄到的強(qiáng)震動數(shù)據(jù)，訓(xùn)練出了適合日本地區(qū)的PGA預(yù)測機(jī)器學(xué)習(xí)模型.Alavi和Gandomi(2011)采用模擬退火算法優(yōu)化用于地震動強(qiáng)度參數(shù)預(yù)測的神經(jīng)網(wǎng)絡(luò)，解決了人工神經(jīng)網(wǎng)絡(luò)參數(shù)學(xué)習(xí)過程中的局部收斂問題.Ameur等(2018)采用隨機(jī)自適應(yīng)模糊神經(jīng)網(wǎng)絡(luò)，同時利用模糊邏輯和人工神經(jīng)網(wǎng)絡(luò)的優(yōu)點(diǎn)，提高訓(xùn)練速度，構(gòu)建了PGA、PGV和PGD預(yù)測機(jī)器學(xué)習(xí)模型.Derakhshani和Foruzan(2019)考慮深度學(xué)習(xí)在模擬復(fù)雜非線性關(guān)系和自動特征提取方面的優(yōu)勢，采用深度學(xué)習(xí)模型預(yù)測了PGA、PGV和PGD.此外，Wiszniowski(2019)還提出采用可逐步訓(xùn)練的Fahlman級聯(lián)相關(guān)神經(jīng)網(wǎng)絡(luò)來構(gòu)建地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型.Hamze-Ziabari和Bakhshpoori(2018)提出采用以分類與回歸樹(Classification and Regression Trees, CART)和M5′模型樹為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)進(jìn)行地震動強(qiáng)度參數(shù)預(yù)測.

以上工作存在三方面問題：(1)震后震害概率評估和地震危險性分析是相關(guān)工作的關(guān)鍵，除了直接給出地震動強(qiáng)度參數(shù)的預(yù)測結(jié)果，相關(guān)結(jié)果的不確定度(或置信程度)也尤為重要，而這是決策樹、神經(jīng)網(wǎng)絡(luò)等單點(diǎn)回歸模型所無法實(shí)現(xiàn)的；(2)演化建模方法中，為保證方程的合理性，其復(fù)雜度受限，導(dǎo)致預(yù)測精度欠佳；(3)基于神經(jīng)網(wǎng)絡(luò)(特別是深度學(xué)習(xí))方法的預(yù)測精度高，然而其可解釋性飽受地球物理領(lǐng)域?qū)W者質(zhì)疑.針對這些問題，本文擬采用自然梯度提升(NGBoost)算法(Duan et al., 2020) 替代決策樹或神經(jīng)網(wǎng)絡(luò)模型，通過給出預(yù)測結(jié)果及置信區(qū)間，更好地完成了地震動強(qiáng)度參數(shù)預(yù)測任務(wù).為解決模型的可解釋性問題，本文結(jié)合SHAP(SHapley Additive exPlanations)值機(jī)器學(xué)習(xí)解釋算法(Lundberg et al., 2020)，分析在PGA和PGV預(yù)測過程中NGBoost模型的各個輸入特征的重要度和它們對模型預(yù)測結(jié)果的影響，以給出機(jī)器學(xué)習(xí)模型的預(yù)測過程.

本文將先介紹訓(xùn)練地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型所使用的強(qiáng)震動數(shù)據(jù)和數(shù)據(jù)篩選方法.接著介紹構(gòu)建概率密度分布預(yù)測機(jī)器學(xué)習(xí)模型所使用NGBoost算法原理.之后介紹地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型的訓(xùn)練過程和結(jié)果，以及基于SHAP值的機(jī)器學(xué)習(xí)模型解釋.最后討論了和其他地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型的比較，以及對各個輸入特征SHAP值變化規(guī)律的物理解釋.

1 NGA-WEST2強(qiáng)震動數(shù)據(jù)庫和數(shù)據(jù)篩選

淺部地殼是具有較大破壞力的地震的多發(fā)區(qū)，為訓(xùn)練出適合活躍構(gòu)造區(qū)淺部地殼地震的地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型，本研究選擇采用NGA-WEST2強(qiáng)震動數(shù)據(jù)庫.NGA-WEST2強(qiáng)震動數(shù)據(jù)庫是太平洋地震工程研究中心(Pacific Earthquake Engineering Research Center)為發(fā)展適合淺部地殼地震的下一代地震動預(yù)測方程而建立的數(shù)據(jù)庫(Bozorgnia et al., 2014)，是目前最完備的強(qiáng)震動數(shù)據(jù)庫.2015年1月新修訂的NGA-WEST2強(qiáng)震動數(shù)據(jù)庫共收錄來自599個地震(包含233個M>5的全球地震和266個3.5

為獲得更好的機(jī)器學(xué)習(xí)模型性能，需要對訓(xùn)練數(shù)據(jù)進(jìn)行篩選與清洗.參考其他地震動預(yù)測方程(Abrahamson et al., 2014; Boore et al., 2014)和地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)研究，本研究采用如下標(biāo)準(zhǔn)從NGA-WEST2強(qiáng)震動數(shù)據(jù)庫中篩選數(shù)據(jù)：(1)只選用包含完整元數(shù)據(jù)信息(震級、震源機(jī)制、地下30 m平均S波速度VS30等)的強(qiáng)震動數(shù)據(jù)；(2)去除S波觸發(fā)的波形不完整的地震動記錄，去除因采樣率或增益異常而導(dǎo)致反應(yīng)譜異常的地震動記錄；(3)去除深震和俯沖帶地震(參見Abrahamson et al.(2013)表2.1)，只選取發(fā)生在活動構(gòu)造區(qū)的地殼地震；(4)按照Abrahamson et al.(2013)公式2.1確定不同年代和地區(qū)數(shù)據(jù)選取的最大震中距(見附圖A1)，去除Joyner-Boore震中距過大的地震記錄；(5)去除地下以及兩層建筑物以上的強(qiáng)震動記錄(參見Boore(2013)表2.1)，只選擇能夠反映自由場地效應(yīng)的地震記錄；(6)由于余震的高頻部分比普通地震系統(tǒng)性地低約20%～40%(Boore and Atkinson,1989,1992； Abrahamson and Silva, 2008)，故去除余震數(shù)據(jù)；(7)去除觀測記錄少于3個的地震記錄.最終，選取了282個地震的12107條記錄用于機(jī)器模型訓(xùn)練.震級分布范圍為3.1～7.9，Joyner-Boore斷層距分布范圍為0.1～348.9 km，VS30分布范圍為89.3～2016.1 m·s-1.具體數(shù)據(jù)分布如圖1所示.

2 地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型訓(xùn)練和結(jié)果分析

2.1 自然梯度提升算法

自然梯度提升算法(Natural Gradient Boosting，NGBoost)(Duan et al., 2020)是一種新的提升算法，它通過引入自然梯度解決了利用提升算法進(jìn)行概率密度分布預(yù)測的難題.自然梯度代表黎曼空間中的最速下降方向(Amari, 1998).相比普通梯度，自然梯度可以使機(jī)器學(xué)習(xí)過程更加高效穩(wěn)定(Duan et al., 2020).

(1)

圖2 用于地震動強(qiáng)度參數(shù)預(yù)測的NGBoost算法示意圖x表示輸入特征；y表示觀測地震動強(qiáng)度參數(shù)；θ(m)、Pθ(m)(y|x)、S(Pθ(m),y)和S(Pθ(m),y)分別表示第m個提升階段的地震動強(qiáng)度參數(shù)概率密度分布的參數(shù)向量、地震動強(qiáng)度參數(shù)概率密度分布、預(yù)測結(jié)果評分和評分的自然梯度; m=0為預(yù)測參數(shù)的初始值；第m個提升階段的基礎(chǔ)學(xué)習(xí)器f(m)(x)由輸入特征x和前一階段評分的自然梯度(Pθ(m-1),y)訓(xùn)練得到.Fig.2 Block diagram of NGBoost algorithm for ground motion parameters predictionx is input features; y is observed ground motion parameters; θ(m), Pθ(m)(y|x), S(Pθ(m),y) and S(Pθ(m),y) are predicted parameters of probability distributions, probability distributions, scores of probability distribution and natural gradients of scores at the mth upgrading stage, respectively. If m is 0, the predicted parameters are initial values. The base learner f(m)(x) at the mth boosting stage is trained with input features x and natural gradients of scores (Pθ(m-1),y) at the previous boosting stage.

其中f(m)(x)為第m個提升段使用的基礎(chǔ)學(xué)習(xí)器，ρ(m)為學(xué)習(xí)器比例因子，η為加權(quán)系數(shù).初始值(θ(0))利用M個基礎(chǔ)學(xué)習(xí)器經(jīng)過M個提升階段的M次修正后得到要預(yù)測的概率密度分布的參數(shù)向量(θ(M)).訓(xùn)練中，后續(xù)基礎(chǔ)學(xué)習(xí)器擬合之前基礎(chǔ)學(xué)習(xí)器和學(xué)習(xí)目標(biāo)間的殘差.訓(xùn)練完成后，各學(xué)習(xí)器按比例因子ρ(m)縮放.參數(shù)θ的類型取決于概率密度分布(如采用正態(tài)分布表示各樣本預(yù)測結(jié)果時，參數(shù)θ=(μ,logσ)).通過評分標(biāo)準(zhǔn)S比較預(yù)測概率分布和觀測數(shù)據(jù)以表示預(yù)測結(jié)果和真實(shí)數(shù)據(jù)間的差異來構(gòu)建損失函數(shù).NGBoost算法的詳細(xì)描述參見Duan等(2020)算法1.相比于同樣可進(jìn)行概率密度分布預(yù)測的貝葉斯深度學(xué)習(xí)方法(Neal, 2012)，NGBoost算法計算量更小.

2.2 機(jī)器學(xué)習(xí)模型訓(xùn)練

本文利用機(jī)器學(xué)習(xí)工具Scikit-Learn將數(shù)據(jù)隨機(jī)分成訓(xùn)練集(80%)和測試集(20%).分別來訓(xùn)練和驗(yàn)證評價機(jī)器學(xué)習(xí)模型.輸入特征的選取是機(jī)器學(xué)習(xí)模型性能的關(guān)鍵，但對于地震動強(qiáng)度參數(shù)預(yù)測研究目前尚缺乏統(tǒng)一標(biāo)準(zhǔn).本研究測試了不同輸入特征組合.首先選取矩震級MW、Joyner-Boore斷層距Rjb、地下30 m平均S波速度VS30、斷層類型FM和斷層頂部距地表距離ZTOR等最常用的特征(Alavi and Gandomi, 2011; Derras et al., 2014; Thomas et al., 2016; Ameur et al., 2018; Dhanya and Raghukanth, 2018; Hamze-Ziabari and Bakhshpoori, 2018; Akhani et al., 2019; Wiszniowski, 2019)作為輸入特征.考慮到斷層類型FM是依據(jù)滑動角Rake定義(Boore et al., 2014)，同樣測試了滑動角Rake和斷層傾角Dip作為輸入特征的情況.此外本文還測試了Z2.5(即S波速度達(dá)到2.5 km·s-1時的深度)作為輸入特征時的模型表現(xiàn)，該特征曾用于模擬淺層沉積物放大效應(yīng)之外的盆地放大效應(yīng)(Campbell and Bozorgnia, 2014)，但相比于其他特征屬性，Z2.5值的缺失率過高(24%)，故使用平均值代替缺失內(nèi)容.應(yīng)力降、斷層破裂速度和方向等數(shù)據(jù)缺失率過高(約70%～90%)，沒有被選作輸入特征.Joyner-Boore斷層距(Rjb)隱含了上下盤效應(yīng)(Boore et al., 2014)，因而沒有增加其他輸入特征模擬上下盤效應(yīng).最終，本文選擇以下7種特征組合作為模型輸入：

(1) 矩震級MW+Joyner-Boore斷層距Rjb;

(2) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30;

(3) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30+斷層頂部深度ZTOR;

(4) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30+斷層頂部深度ZTOR+斷層類型FM;

(5) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30+斷層頂部深度ZTOR+滑動角Rake;

(6) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30+斷層頂部深度ZTOR+滑動角Rake+斷層傾角Dip;

(7) 矩震級MW+ Joyner-Boore斷層距Rjb+地下30 m平均S波速度VS30+斷層頂部深度ZTOR+滑動角Rake+斷層傾角Dip+VS達(dá)到2.5 km·s-1的深度Z2.5.

由于NGBoost的基礎(chǔ)學(xué)習(xí)器CART對于特征縮放不敏感，故沒有對特征進(jìn)行線性函數(shù)歸一化和標(biāo)準(zhǔn)化.本研究使用K-折交叉驗(yàn)證和網(wǎng)格搜索確定最優(yōu)的超參數(shù)：基礎(chǔ)學(xué)習(xí)器最大深度(max depth)、學(xué)習(xí)率(learning rate)和基礎(chǔ)學(xué)習(xí)器個數(shù)(number of estimators).基礎(chǔ)學(xué)習(xí)器最大深度取值在3～10；學(xué)習(xí)率取Duan等(2020)的推薦值0.1和0.01.基礎(chǔ)學(xué)習(xí)器個數(shù)統(tǒng)一取為4000次，選取使得驗(yàn)證數(shù)據(jù)集平均評分最小的基礎(chǔ)學(xué)習(xí)器個數(shù).K-折交叉驗(yàn)證實(shí)驗(yàn)中，訓(xùn)練數(shù)據(jù)集被隨機(jī)分成5組，當(dāng)其中1組用于驗(yàn)證時，剩余4組用于訓(xùn)練.完成訓(xùn)練后，取使5次驗(yàn)證誤差均值最小的模型訓(xùn)練整個訓(xùn)練集數(shù)據(jù).表1展示了選取不同特征組合作為輸入時，最佳超參數(shù)下的模型表現(xiàn).可見：(1)選用越多的特征獲得的機(jī)器學(xué)習(xí)模型性能越優(yōu)；(2)使用滑動角作為輸入特征的效果優(yōu)于使用斷層類型；(3)當(dāng)基礎(chǔ)學(xué)習(xí)器最大深度分別取7和6，學(xué)習(xí)率取0.01，基礎(chǔ)學(xué)習(xí)器個數(shù)分別取233和284時，PGA和PGV預(yù)測機(jī)器學(xué)習(xí)模型性能達(dá)到最優(yōu).圖3進(jìn)一步展示了當(dāng)使用所有特征(MW+Rjb+VS30+ZTOR+Rake+Dip+Z2.5)時，驗(yàn)證數(shù)據(jù)集的平均評分隨基礎(chǔ)學(xué)習(xí)器最大深度和學(xué)習(xí)率的變化.最終我們將所有特征屬性作為輸入，利用驗(yàn)證得到的最佳超參數(shù)訓(xùn)練PGA和PGV預(yù)測模型.

2.3 機(jī)器學(xué)習(xí)模型性能分析

圖4顯示了模型預(yù)測結(jié)果和實(shí)際結(jié)果比較，可以看出：(1)模型表現(xiàn)良好，預(yù)測結(jié)果在訓(xùn)練集和測試集上都與真實(shí)值接近；(2)基于NGBoost算法的地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型具有良好的泛化能力(訓(xùn)練集和測試集的預(yù)測表現(xiàn)接近).為進(jìn)一步分析機(jī)器學(xué)習(xí)模型的性能，計算了機(jī)器學(xué)習(xí)模型預(yù)測值與實(shí)際觀測值的比值、相關(guān)系數(shù)R、平均絕對誤差MAE和均方根誤差RMSE，并與其他地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)研究進(jìn)行了比較.測試集結(jié)果顯示：(1)PGA和PGV的預(yù)測值與觀測值間誤差在±30%以內(nèi)的比例分別約為70%和80%(見圖5)，優(yōu)于Derakhshani和Foruzan(2019)利用深度學(xué)習(xí)和NGA-WEST2強(qiáng)震動數(shù)據(jù)庫訓(xùn)練得到的模型(原文分別為65%和60%)；(2)PGA和PGV預(yù)測值與實(shí)測值的相關(guān)系數(shù)分別為0.960和0.976，平均絕對誤差分別為0.511和0.451，均方根誤差分別為0.655和0.585，預(yù)測結(jié)果的相關(guān)系數(shù)均優(yōu)于Derakhshani和Foruzan(2019)深度學(xué)習(xí)的結(jié)果(表2為本模型的性能，表3為Derakhshani和Foruzan(2019)論文的結(jié)果).

表1 對于PGA和PGV預(yù)測，不同輸入特征組合對應(yīng)的最佳超參數(shù)和驗(yàn)證數(shù)據(jù)集平均評分Table 1 For prediction of PGA and PGV, the best hyperparameters and corresponding average scores of validation dataset for different set of input features

圖3 對于PGA和PGV預(yù)測，當(dāng)學(xué)習(xí)率(0.1和0.01)和基礎(chǔ)學(xué)習(xí)器最大深度(3～10)取不同值時，驗(yàn)證數(shù)據(jù)集平均評分隨基礎(chǔ)學(xué)習(xí)器個數(shù)的變化.訓(xùn)練機(jī)器學(xué)習(xí)模型預(yù)測PGA和PGV時使用了所有特征屬性(MW+Rjb+VS30+ZTOR+Rake+Dip+Z2.5)圖中黑點(diǎn)為基礎(chǔ)學(xué)習(xí)器最大深度取不同值的驗(yàn)證數(shù)據(jù)集的最小平均評分.Fig.3 For prediction of PGA and PGV, the variations of average scores of validation dataset with the number of estimators for different learning rate (0.1 and 0.01) and max depth of estimator (3～10). All the input features (MW+Rjb+VS30+ZTOR+Rake+Dip+Z2.5) are used when training these models for prediction of PGA and PGVThe black dots represent minimum average scores of validation dataset for different max depth of estimator.

圖4 機(jī)器學(xué)習(xí)模型預(yù)測的PGA和PGV與實(shí)際觀測值的比較. 模型的訓(xùn)練集(80%)和測試集(20%)從NGA-WEST2強(qiáng)震動數(shù)據(jù)庫中隨機(jī)選取Fig.4 Comparisons between the machine learning model for predicted and observed PGA and PGV. The training (80%) and testing (20%) dataset are randomly selected from NGA-WEST2 strong motion dataset

圖5 機(jī)器學(xué)習(xí)模型預(yù)測PGA和PGV與觀測值的比值的統(tǒng)計直方圖. 模型的訓(xùn)練集(80%)和測試集(20%)從NGA-WEST2強(qiáng)震動數(shù)據(jù)庫中隨機(jī)選取Fig.5 Histograms of ratios of machine learning model for predicted and observed PGA and PGV. The training (80%) and testing (20%) dataset are randomly selected from NGA-WEST2 strong motion dataset

表2 NGBoost算法訓(xùn)練出的PGA和PGV預(yù)測機(jī)器學(xué)習(xí)模型的性能指標(biāo). 模型的訓(xùn)練集(80%)和測試集(20%)從NGA-WEST2強(qiáng)震動數(shù)據(jù)庫中隨機(jī)選取Table 2 Performance of the trained NGBoost algorithm for PGA and PGV prediction. The training (80%) and testing (20%) dataset are randomly selected from NGA-WEST2 strong motion dataset

表3 基于深度學(xué)習(xí)的PGA和PGV預(yù)測結(jié)果指標(biāo)(Derakhshani and Foruzan，2019)Table 3 Performance of PGA and PGV prediction by a Deep-Learning method (Derakhshani and Foruzan, 2019)

相比其他地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型，基于NGBoost算法的機(jī)器學(xué)習(xí)模型不僅可以給出準(zhǔn)確預(yù)測值而且可以給出預(yù)測結(jié)果的概率密度分布.圖6展示了對于2004年美國加州ParkfieldMW6.0地震(EQID=1018)、2009年美國加州San BernardinoMW4.45地震(EQID=1018)和2007年日本Chuetsu-okiMW6.8地震(EQID=278)，機(jī)器學(xué)習(xí)模型預(yù)測的PGA和PGV的概率密度分布.約84%的實(shí)際觀測值位于機(jī)器學(xué)習(xí)模型預(yù)測的85%置信區(qū)間內(nèi)，說明機(jī)器學(xué)習(xí)模型可以很好地預(yù)測地震動的概率密度分布，且對于近場和遠(yuǎn)場都有很好預(yù)測.

圖6 基于NGBoost機(jī)器學(xué)習(xí)模型的PGA和PGV的概率密度分布預(yù)測. 數(shù)據(jù)為2004年美國加州Parkfield MW6.0地震(EQID=1018)、2009年美國加州San Bernardino MW4.45地震(EQID=1018)和2007年日本Chuetsu-oki MW6.8地震(EQID=278)圖中黑點(diǎn)和青點(diǎn)分別為訓(xùn)練集和測試集中的觀測值. 藍(lán)線表示機(jī)器學(xué)習(xí)模型預(yù)測值，其中綠線表示85%置信區(qū)間. 色標(biāo)表示預(yù)測的概率密度分布.Fig.6 The predicted probability density distributions by the proposed NGBoost model. Dataset are the 2004 Parkfield MW6.0 earthquake, 2009 San Bernardino MW4.45 earthquake and 2007 Chuetsu-oki MW6.8 earthquakeThe black and cyan dots are the observations in the training and testing datasets, respectively. The blue lines are the predicted values by machine learning models where the green lines indicate the 85% confidence intervals. The predicted probability densities are color-coded shown in the color bar.

3 基于SHAP值的地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型解釋

復(fù)雜模型(如深度學(xué)習(xí))在地震動參數(shù)預(yù)測任務(wù)中給出了相當(dāng)準(zhǔn)確的預(yù)測結(jié)果(Derakhshani and Foruzan, 2019)，但地震專家卻難以分析這些復(fù)雜預(yù)測過程與客觀事實(shí)的相符程度.事實(shí)上，以CART作為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)模型在訓(xùn)練完成后，可通過基尼指數(shù)(Gini importance)或者Saabas 值等方式來計算特征對模型輸出的全局重要性(由于篇幅的限制本文不展開討論).但仍難以給出局部樣本對模型輸出的影響，且無法量化不同特征之間的相關(guān)性(相互耦合)，不足以為預(yù)測過程給出確切的解釋.

為理解模型對地震動強(qiáng)度參數(shù)預(yù)測的機(jī)理，我們引入了一種統(tǒng)一的解釋框架：SHAP(SHapley Additive exPlanations)(Lundberg and Lee, 2017; Lundberg et al., 2019, 2020).對于NGBoost這類決策樹集成模型，SHAP框架結(jié)合了局部解釋方法LIME(Local Interpretable Model-agnostic Explanations)(Ribeiro et al., 2016)和經(jīng)典Shapley值估計方法(Shapley, 1953).基于LIME方法，SHAP框架利用若干個線性模型對待解釋復(fù)雜模型進(jìn)行逼近，并逐個計算這些線性模型中各個特征的Shapley重要性估計值，以此定量表示預(yù)測過程中，每個特征對預(yù)測的貢獻(xiàn).對于帶有M個特征的簡化輸入x′，貢獻(xiàn)值g(x′)定義如式(2)，其中，φ0為平均貢獻(xiàn)值，φj為特征j的Shapley值.當(dāng)令φ0為模型輸出的期望E(f(x))時，SHAP輸出將近似等于模型的真實(shí)輸出.SHAP框架的詳細(xì)介紹見Lundberg和Lee(2017)、Lundberg et al.(2019,2020)，限于篇幅，不再贅述.

(2)

以2007年日本Chuetsu-okiMW6.8地震3條記錄(Joyner-Boore斷層距分別為20.4、107.1和279.3 km)為例，本文利用SHAP框架對PGA預(yù)測結(jié)果進(jìn)行初步解釋.3條記錄中，各個輸入特征的SHAP解釋結(jié)果如圖7所示.可以看出：(1)矩震級MW和Rjb距離對PGA預(yù)測貢獻(xiàn)遠(yuǎn)大于其他特征(SHAP值的絕對值大)；(2)Rjb距離對PGA的SHAP值(即貢獻(xiàn)值)隨距離增大而逐漸減少，當(dāng)Rjb距離較大時，其對于PGA預(yù)測的貢獻(xiàn)為負(fù).這些與地震學(xué)理論相符，證實(shí)了SHAP值機(jī)器學(xué)習(xí)解釋算法的有效性.

為進(jìn)一步分析各個特征(MW、Rjb、VS30、ZTOR、Rake、Dip和Z2.5)在整體上如何影響PGA和PGV預(yù)測，我們計算了所有地震動記錄的各特征對PGA和PGV的SHAP值，并繪制了SHAP值摘要圖(圖8).在SHAP值摘要圖中，各個特征按重要度(SHAP值絕對值的平均值)由上至下排列.對于圖中的每一個特征，每個點(diǎn)表示每一次預(yù)測該特征的SHAP值.所有預(yù)測按SHAP值大小沿X軸排列，并在出現(xiàn)相同值時沿Y軸堆疊以顯示數(shù)據(jù)分布.數(shù)據(jù)點(diǎn)顏色表示輸入特征值的大小.如果顏色隨SHAP值連續(xù)變化說明輸入特征對預(yù)測的貢獻(xiàn)隨輸入特征值的大小連續(xù)變化.

圖7 用于PGA預(yù)測的不同輸入特征的SHAP值統(tǒng)計結(jié)果(以2007年日本Chuetsu-oki MW6.8地震的三條Joyner-Boore震中距Rjb為20.4、107.1和279.3 km記錄為例)圖中虛線標(biāo)注了對于不同記錄的機(jī)器學(xué)習(xí)模型預(yù)測值，柱狀圖中各部分的顏色長度表示SHAP值的大小，顏色表示SHAP值的正負(fù)(紅色為正，藍(lán)色為負(fù)).Fig.7 The computed SHAP values of different input features for three strong motion records (the 2007 Chuetsu-oki MW6.8 earthquake with distances at 20.4, 107.1 and 279.3 km)The dashed lines are the predicted PGAs by trained machine learning model where the height of each bar indicates the magnitudes of SHAP values. The color represents the signs of SHAP values (red color indicates positive value, and blue color indicates negative value).

圖8 對于PGA和PGV預(yù)測，所有輸入特征的SHAP值摘要圖圖中各個特征后的數(shù)值表示特征的重要度，每個點(diǎn)表示每次預(yù)測不同特征的SHAP值，顏色代表特征的值.這些點(diǎn)按SHAP值大小沿X軸排列，并在出現(xiàn)相同值時沿Y軸堆疊以顯示數(shù)據(jù)分布.Fig.8 For prediction of PGA and PGV, the SHAP summary plot of all input featuresThe number behind input features are the feature importance. The dots are SHAP values of different features for different predictions. The color indicates the value of input feature. These points are plotted horizontally, and stacking vertically when they have the same SHAP value to show the distribution of data.

SHAP值摘要圖展示了一些工程地震學(xué)中常見可以說明機(jī)器學(xué)習(xí)模型合理性的現(xiàn)象：矩震級MW的特征重要性隨震級增大而增大，較高的矩震級MW普遍帶來較大的正向貢獻(xiàn)，而當(dāng)矩震級較低時，其對地震動參數(shù)的影響則不如高值那樣穩(wěn)定(藍(lán)色范圍明顯寬于紅色范圍)；Rjb距離、VS30的特征重要性則隨特征值增大而減小，且較低的Rjb、VS30對預(yù)測結(jié)果的貢獻(xiàn)值不穩(wěn)定(藍(lán)色范圍明顯寬于其他顏色).此外，圖中還揭示出一些可能被忽視的現(xiàn)象：(1)Rjb對PGA預(yù)測的重要度低于其對PGV的重要度；(2)斷層頂部深度ZTOR的SHAP值在PGA和PGV的預(yù)測結(jié)果解釋中存在“長尾”現(xiàn)象(見ZTOR紅色點(diǎn))，這意味著模型認(rèn)為較大地震會使得PGA和PGV的預(yù)測結(jié)果更大.

4 討論

基于NGA-WEST2強(qiáng)震動數(shù)據(jù)庫，本文利用NGBoost算法訓(xùn)練出了可用于PGA和PGV概率密度分布預(yù)測的機(jī)器學(xué)習(xí)模型.考慮到平均絕對誤差和均方根誤差受單位影響較大，本研究采用預(yù)測結(jié)果與實(shí)際觀測值(包括測試集和訓(xùn)練集)的相關(guān)系數(shù)作為評價指標(biāo)，以便于同已有研究進(jìn)行比較.

與深度學(xué)習(xí)(Derakhshani and Foruzan, 2019)、以CART為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)(Hamze-Ziabari and Bakhshpoori, 2018)、M5′模型樹(Kaveh et al., 2016)、以M5′模型樹為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)(Hamze-Ziabari and Bakhshpoori, 2018)、單隱藏層神經(jīng)網(wǎng)絡(luò)(Alavi and Gandomi, 2011)、最小二乘和遺傳編程混合算法(Gandomi et al., 2011)、多表達(dá)式編程(Alavi et al., 2011)、模擬退火和遺傳編程混合算法(Mohammadnejad et al., 2012)、隨機(jī)自適應(yīng)神經(jīng)網(wǎng)絡(luò)模糊推理系統(tǒng)(Thomas et al., 2016)等諸多方法相比，NGBoost算法的預(yù)測值與觀測值具有更高的相關(guān)系數(shù)(表4).表中大部分研究采用了和本研究類似的數(shù)據(jù)篩選方法，并且NGA-WEST1強(qiáng)震動數(shù)據(jù)庫包含在NGA-WEST2強(qiáng)震動數(shù)據(jù)庫之中.因此，數(shù)據(jù)對于比較結(jié)果的影響可以忽略不計.

表4 不同研究機(jī)器學(xué)習(xí)模型預(yù)測PGA和PGV與觀測值的相關(guān)系數(shù)(R)的比較

此外，Hamze-Ziabari和Bakhshpoori(2018)指出，使用M5′模型樹做基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí)算法表現(xiàn)優(yōu)于以CART為基礎(chǔ)學(xué)習(xí)器的集成學(xué)習(xí).相比于CART，M5′模型樹的葉節(jié)點(diǎn)使用線性模型，并且經(jīng)過剪枝和平滑處理后，單棵M5′模型樹的性能更優(yōu).本研究對原有NGBoost算法進(jìn)行改進(jìn)，使用M5′模型樹作為基礎(chǔ)學(xué)習(xí)器進(jìn)行了相同的實(shí)驗(yàn).圖A2展示了使用M5′模型樹和CART作為基礎(chǔ)學(xué)習(xí)器構(gòu)建NGBoost模型，對2004年美國加州ParkfieldMW6.0地震、2009年美國加州San BernardinoMW4.45地震和2007年日本Chuetsu-okiMW6.8地震預(yù)測結(jié)果的比較.結(jié)合表4所示，兩者預(yù)測結(jié)果接近，表明M5′模型樹并不能顯著提升NGBoost算法性能.

基于NGBoost算法的地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型不但具有高的預(yù)測精度，而且通過計算SHAP值還可以從數(shù)據(jù)角度分析輸入特征如何影響地震動強(qiáng)度參數(shù)預(yù)測.圖9和圖10分別給出所建立的機(jī)器學(xué)習(xí)模型對PGA和PGV預(yù)測時各輸入特征與其對應(yīng)的SHAP值之間的關(guān)系.

(1)圖9a和圖10a中，隨著矩震級MW增加其對應(yīng)的SHAP值呈近似分段線性增加(虛線為分段處)，且增加的幅度逐漸減少.對近場(<2 km)PGA預(yù)測(圖A3a)，MW對地震動強(qiáng)度預(yù)測的貢獻(xiàn)值不隨MW值增大而增大，表明PGA隨MW的變化是非線性的，存在“近震飽和”現(xiàn)象.“近震飽和”現(xiàn)象是由于對于近場記錄，斷層破裂面上每個子斷層產(chǎn)生的地震波的持續(xù)時間都較短，并且到達(dá)觀測點(diǎn)的時間也不同，很難線性疊加產(chǎn)生較大振幅(Yenier and Atkinson, 2014).這些分段現(xiàn)象也說明了在某些地震動預(yù)測方程中震級項(xiàng)使用分段線性函數(shù)的合理性(Campbell and Bozorgnia, 2014).

圖9 PGA預(yù)測過程中各個輸入特征的SHAP值依賴圖.SHAP值依賴圖展示了輸入特征值變化對模型預(yù)測帶來的影響橫坐標(biāo)表示各個特征的值，縱坐標(biāo)表示各個特征的SHAP值.顏色表示對所考察特征依賴程度最大的其他特征的取值.(g)中的N、S和R分別代表正斷層(-150°

圖10 PGV預(yù)測過程中各個輸入特征的SHAP值依賴圖.SHAP值依賴圖展示了輸入特征變化如何影響模型預(yù)測橫坐標(biāo)表示各個特征的值，縱坐標(biāo)表示各個特征的SHAP值.顏色表示對探究特征影響最大的另一個特征的值(g)中的N、S和R分別代表正斷層(-150°

(2)隨著Rjb距離逐漸增大，Rjb的貢獻(xiàn)值不斷減小，并在約45和150 km后兩次出現(xiàn)減小速率放緩(圖9b和圖10b).45 km后的減小速率放緩可能與面波或Lg波出現(xiàn)相關(guān).150 km后出現(xiàn)的減小速率放緩可能與S波的莫霍面反射波(SmS波)出現(xiàn)相關(guān)(Zhu et al., 2019).

(3)VS30的SHAP值并非一直隨VS30值的減小而增加(圖9c和圖10c).當(dāng)VS30值低于240 m·s-1(圖中虛線對應(yīng)的值)時，其對應(yīng)的SHAP值會出現(xiàn)降低，且SHAP值降低的點(diǎn)MW值較高.這與Aboye等(2011)利用1D等效線性模型計算出的場地放大系數(shù)隨VS30值的變化規(guī)律是一致的，可能是由于地震動強(qiáng)度過高時出現(xiàn)的場地非線性效應(yīng)所致.

除了分析主要輸入特征影響外，ZTOR、Z2.5、Dip和Rake等次重要特征的SHAP值依賴圖也展示出了一些有趣的現(xiàn)象.例如對于活躍構(gòu)造區(qū)地殼地震，斷層較淺時(ZTOR<～5 km)的ZTOR的SHAP值整體上低于斷層較深時的值(圖9d和圖10d).這說明淺層的軟弱帶破裂受速度強(qiáng)化控制，滑動弱化距離更大，需要更多的破裂能，通常會導(dǎo)致較低的破裂速度和滑動速率，進(jìn)而導(dǎo)致較低的PGA和PGV(Somerville and Pitarka, 2006).當(dāng)斷層較淺時，ZTOR的SHAP值隨ZTOR值增大而減小(圖9d和圖10d)，可能是由于發(fā)震深度增加導(dǎo)致的幾何衰減的增大所致.當(dāng)斷層較深時，ZTOR的SHAP值隨ZTOR值增大而增大(圖9d和圖10d)，可能是由于發(fā)震深度增加導(dǎo)致的地震應(yīng)力降的增加或品質(zhì)因子Q的增大所致(Abercrombie et al., 2021).當(dāng)Z2.5>～3 km時，對于PGA和PGV預(yù)測，Z2.5的SHAP值均隨Z2.5值的增大而增大(圖9e和圖10e)，此時Z2.5主要反映3D盆地效應(yīng)，對PGA和PGV都是放大作用(Campbell and Bozorgnia, 2013).當(dāng)Z2.5<～1 km時，對于PGA預(yù)測，Z2.5的SHAP值隨Z2.5值的增大而減小(圖9e)，而對于PGV預(yù)測，Z2.5的SHAP值隨Z2.5值的增大而增大(圖10e).此時，Z2.5主要反映淺層沉積物效應(yīng)(Campbell and Bozorgnia, 2013).淺層厚沉積物的自振頻率低，離頻率較低的速度的頻率更近，離頻率較高的加速度的頻率更遠(yuǎn)，從而導(dǎo)致PGV的增加和PGA的減小.SHAP值摘要圖中觀測到的Rjb對PGA預(yù)測的重要度低于PGV也可能與位移和速度的頻率不同相關(guān)(圖8).頻率越高受幾何衰減影響越明顯.對于PGA和PGV預(yù)測，Dip的SHAP值均隨Dip值緩慢增加(圖9f和圖10f)，表明Dip越大預(yù)測PGA和PGV越大.對于PGA和PGV預(yù)測，Rake的SHAP值隨Rake值近似呈周期性變化，并且當(dāng)?shù)卣馂樽呋瑪鄬?-180°

5 結(jié)論

基于NGA-WEST2強(qiáng)震動數(shù)據(jù)庫，本文利用NGBoost算法訓(xùn)練出了適合預(yù)測活躍構(gòu)造區(qū)地殼地震PGA和PGV的機(jī)器學(xué)習(xí)模型.在測試集數(shù)據(jù)上，PGA和PGV預(yù)測結(jié)果和實(shí)際觀測結(jié)果的相關(guān)系數(shù)可達(dá)0.972和0.984，優(yōu)于目前所有基于NGA-WEST2或NGA-WEST1強(qiáng)震動數(shù)據(jù)庫的地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)模型.得益于NGBoost的預(yù)測機(jī)制，模型除了直接給出實(shí)值結(jié)果，還能為每個預(yù)測結(jié)果提供合理概率密度分布，進(jìn)而為地震危險性概率分析提供定量分析依據(jù).這是目前其他地震動強(qiáng)度參數(shù)預(yù)測機(jī)器學(xué)習(xí)研究所忽視的.

為解釋機(jī)器學(xué)習(xí)模型，分析機(jī)器學(xué)習(xí)模型的預(yù)測過程是否符合內(nèi)在物理規(guī)律和實(shí)際觀測，研究各個特征對預(yù)測結(jié)果的影響和重要程度，計算了各個特征的SHAP值.SHAP值顯示的各個輸入特征對PGA和PGV預(yù)測的影響基本與地震動預(yù)測方程研究一致，但同時也展示出一些以往研究所忽視的現(xiàn)象.例如，對于活躍構(gòu)造區(qū)地殼地震，淺部破裂(ZTOR<～5 km)ZTOR的SHAP值整體上低于深部，且多為負(fù)值，表明相比于深部破裂淺部破裂有著更低的PGA和PGV.這可能是由于淺部破裂受速度強(qiáng)化控制所致.并且淺部ZTOR的SHAP值隨ZTOR值增大而減小，深部ZTOR的SHAP值隨ZTOR值增大而增大.這可能與隨發(fā)震深度變化的幾何衰減、應(yīng)力降和品質(zhì)因子Q的變化相關(guān).當(dāng)Z2.5較小(Z2.5<～1 km)時，Z2.5對PGA和PGV的影響是相反的.可能是因?yàn)楫?dāng)Z2.5較小時，Z2.5主要反映淺層沉積物效應(yīng)而不再是3D盆地效應(yīng).淺層沉積物厚度的變化會影響共振頻率的變化.加速度和速度頻率不同，受到的影響不同.Rjb對PGA和PGV預(yù)測的重要度不同，也可能與加速度和速度的頻率不同相關(guān).

致謝謹(jǐn)此祝賀陳颙先生從事地球物理教學(xué)科研工作60周年．

附錄

圖A1 對于不同年代(1933—2000，2001—2005和2006—2011年)和不同震級的地震，選取地震記錄時使用的最大Joyner-Boore震中距Fig.A1 For different times (1933—2000, 2001—2005 and 2006—2011) and moment magnitudes, the maximum Joyner-Boore epicentral distances used for selecting seismic data

圖A2 對于2004年美國加州Parkfield MW6.0地震(EQID=1018)、2009年美國加州San Bernardino MW4.45地震(EQID=1018)和2007年日本Chuetsu-oki MW6.8地震(EQID=278),使用M 5′模型樹作為基礎(chǔ)學(xué)習(xí)器和使用分類回歸樹作為基礎(chǔ)學(xué)習(xí)器的預(yù)測結(jié)果比較圖中黑點(diǎn)表示測試集中的觀測值，虛線表示使用M 5′模型樹作為基礎(chǔ)學(xué)習(xí)器的預(yù)測結(jié)果，實(shí)線表示使用分類回歸樹作為基礎(chǔ)學(xué)習(xí)器的預(yù)測結(jié)果.Fig.A2 For the 2004 Parkfield MW6.0 earthquake, 2009 SanBernardino MW4.45 earthquake and 2007 Chuetsu-oki MW6.8 earthquake, the comparisons of prediction results using M 5′ model trees and classification and regression trees as base learnersThe black dots indicate observations in testing dataset, the dashed lines indicate predicted results by trained models using M 5′ model tree as base learner, the solid lines indicate predicted results by trained models using classification and regression tree as base learner.