地震動強度參數(shù)估計的可解釋性與不確定度機器學習模型

陳蒙，王華,2*

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

0 引言

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

計算(或預測)地震動強度參數(shù)的方法主要有三種：數(shù)值模擬、地震動預測方程(Ground Motion Prediction Equation，也稱作地震動衰減關(guān)系)和機器學習.常用的數(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ù)值方法可以模擬斷層破裂和波場傳播的全過程，具有明確的物理意義.但是準確的地震動模擬需要精確的震源信息和地下結(jié)構(gòu)速度模型，并且計算量巨大.地震動預測方程(Boore et al., 2014; Campbell and Bozorgnia, 2014; Chiou and Youngs, 2014; Idriss, 2014)形式明確、計算速度快，是地震動圖和地震危險性分析中常用的方法.但是現(xiàn)代地震動預測方程復雜(Abrahamson et al., 2014; Campbell amd Bozorgnia, 2014)，函數(shù)形式和特征變量的選取沒有統(tǒng)一標準，依賴研究人員的主觀性，且易忽略各特征之間的非線性耦合作用.

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

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

(2)基于監(jiān)督學習回歸模型的地震動強度參數(shù)預測

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

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

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

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

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

為獲得更好的機器學習模型性能，需要對訓練數(shù)據(jù)進行篩選與清洗.參考其他地震動預測方程(Abrahamson et al., 2014; Boore et al., 2014)和地震動強度參數(shù)預測機器學習研究，本研究采用如下標準從NGA-WEST2強震動數(shù)據(jù)庫中篩選數(shù)據(jù)：(1)只選用包含完整元數(shù)據(jù)信息(震級、震源機制、地下30 m平均S波速度VS30等)的強震動數(shù)據(jù)；(2)去除S波觸發(fā)的波形不完整的地震動記錄，去除因采樣率或增益異常而導致反應譜異常的地震動記錄；(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)去除地下以及兩層建筑物以上的強震動記錄(參見Boore(2013)表2.1)，只選擇能夠反映自由場地效應的地震記錄；(6)由于余震的高頻部分比普通地震系統(tǒng)性地低約20%～40%(Boore and Atkinson,1989,1992； Abrahamson and Silva, 2008)，故去除余震數(shù)據(jù)；(7)去除觀測記錄少于3個的地震記錄.最終，選取了282個地震的12107條記錄用于機器模型訓練.震級分布范圍為3.1～7.9，Joyner-Boore斷層距分布范圍為0.1～348.9 km，VS30分布范圍為89.3～2016.1 m·s-1.具體數(shù)據(jù)分布如圖1所示.

2 地震動強度參數(shù)預測機器學習模型訓練和結(jié)果分析

2.1 自然梯度提升算法

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

(1)

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

2.2 機器學習模型訓練

本文利用機器學習工具Scikit-Learn將數(shù)據(jù)隨機分成訓練集(80%)和測試集(20%).分別來訓練和驗證評價機器學習模型.輸入特征的選取是機器學習模型性能的關(guān)鍵，但對于地震動強度參數(shù)預測研究目前尚缺乏統(tǒng)一標準.本研究測試了不同輸入特征組合.首先選取矩震級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波速度達到2.5 km·s-1時的深度)作為輸入特征時的模型表現(xiàn)，該特征曾用于模擬淺層沉積物放大效應之外的盆地放大效應(Campbell and Bozorgnia, 2014)，但相比于其他特征屬性，Z2.5值的缺失率過高(24%)，故使用平均值代替缺失內(nèi)容.應力降、斷層破裂速度和方向等數(shù)據(jù)缺失率過高(約70%～90%)，沒有被選作輸入特征.Joyner-Boore斷層距(Rjb)隱含了上下盤效應(Boore et al., 2014)，因而沒有增加其他輸入特征模擬上下盤效應.最終，本文選擇以下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達到2.5 km·s-1的深度Z2.5.

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

2.3 機器學習模型性能分析

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

表1 對于PGA和PGV預測，不同輸入特征組合對應的最佳超參數(shù)和驗證數(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預測，當學習率(0.1和0.01)和基礎(chǔ)學習器最大深度(3～10)取不同值時，驗證數(shù)據(jù)集平均評分隨基礎(chǔ)學習器個數(shù)的變化.訓練機器學習模型預測PGA和PGV時使用了所有特征屬性(MW+Rjb+VS30+ZTOR+Rake+Dip+Z2.5)圖中黑點為基礎(chǔ)學習器最大深度取不同值的驗證數(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 機器學習模型預測的PGA和PGV與實際觀測值的比較. 模型的訓練集(80%)和測試集(20%)從NGA-WEST2強震動數(shù)據(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 機器學習模型預測PGA和PGV與觀測值的比值的統(tǒng)計直方圖. 模型的訓練集(80%)和測試集(20%)從NGA-WEST2強震動數(shù)據(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算法訓練出的PGA和PGV預測機器學習模型的性能指標. 模型的訓練集(80%)和測試集(20%)從NGA-WEST2強震動數(shù)據(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 基于深度學習的PGA和PGV預測結(jié)果指標(Derakhshani and Foruzan，2019)Table 3 Performance of PGA and PGV prediction by a Deep-Learning method (Derakhshani and Foruzan, 2019)

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

圖6 基于NGBoost機器學習模型的PGA和PGV的概率密度分布預測. 數(shù)據(jù)為2004年美國加州Parkfield MW6.0地震(EQID=1018)、2009年美國加州San Bernardino MW4.45地震(EQID=1018)和2007年日本Chuetsu-oki MW6.8地震(EQID=278)圖中黑點和青點分別為訓練集和測試集中的觀測值. 藍線表示機器學習模型預測值，其中綠線表示85%置信區(qū)間. 色標表示預測的概率密度分布.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值的地震動強度參數(shù)預測機器學習模型解釋

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

為理解模型對地震動強度參數(shù)預測的機理，我們引入了一種統(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框架利用若干個線性模型對待解釋復雜模型進行逼近，并逐個計算這些線性模型中各個特征的Shapley重要性估計值，以此定量表示預測過程中，每個特征對預測的貢獻.對于帶有M個特征的簡化輸入x′，貢獻值g(x′)定義如式(2)，其中，φ0為平均貢獻值，φj為特征j的Shapley值.當令φ0為模型輸出的期望E(f(x))時，SHAP輸出將近似等于模型的真實輸出.SHAP框架的詳細介紹見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預測結(jié)果進行初步解釋.3條記錄中，各個輸入特征的SHAP解釋結(jié)果如圖7所示.可以看出：(1)矩震級MW和Rjb距離對PGA預測貢獻遠大于其他特征(SHAP值的絕對值大)；(2)Rjb距離對PGA的SHAP值(即貢獻值)隨距離增大而逐漸減少，當Rjb距離較大時，其對于PGA預測的貢獻為負.這些與地震學理論相符，證實了SHAP值機器學習解釋算法的有效性.

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

圖7 用于PGA預測的不同輸入特征的SHAP值統(tǒng)計結(jié)果(以2007年日本Chuetsu-oki MW6.8地震的三條Joyner-Boore震中距Rjb為20.4、107.1和279.3 km記錄為例)圖中虛線標注了對于不同記錄的機器學習模型預測值，柱狀圖中各部分的顏色長度表示SHAP值的大小，顏色表示SHAP值的正負(紅色為正，藍色為負).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預測，所有輸入特征的SHAP值摘要圖圖中各個特征后的數(shù)值表示特征的重要度，每個點表示每次預測不同特征的SHAP值，顏色代表特征的值.這些點按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值摘要圖展示了一些工程地震學中常見可以說明機器學習模型合理性的現(xiàn)象：矩震級MW的特征重要性隨震級增大而增大，較高的矩震級MW普遍帶來較大的正向貢獻，而當矩震級較低時，其對地震動參數(shù)的影響則不如高值那樣穩(wěn)定(藍色范圍明顯寬于紅色范圍)；Rjb距離、VS30的特征重要性則隨特征值增大而減小，且較低的Rjb、VS30對預測結(jié)果的貢獻值不穩(wěn)定(藍色范圍明顯寬于其他顏色).此外，圖中還揭示出一些可能被忽視的現(xiàn)象：(1)Rjb對PGA預測的重要度低于其對PGV的重要度；(2)斷層頂部深度ZTOR的SHAP值在PGA和PGV的預測結(jié)果解釋中存在“長尾”現(xiàn)象(見ZTOR紅色點)，這意味著模型認為較大地震會使得PGA和PGV的預測結(jié)果更大.

4 討論

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

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

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

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

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

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

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

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

(2)隨著Rjb距離逐漸增大，Rjb的貢獻值不斷減小，并在約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).當VS30值低于240 m·s-1(圖中虛線對應的值)時，其對應的SHAP值會出現(xiàn)降低，且SHAP值降低的點MW值較高.這與Aboye等(2011)利用1D等效線性模型計算出的場地放大系數(shù)隨VS30值的變化規(guī)律是一致的，可能是由于地震動強度過高時出現(xiàn)的場地非線性效應所致.

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

5 結(jié)論

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

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

致謝謹此祝賀陳颙先生從事地球物理教學科研工作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ǔ)學習器和使用分類回歸樹作為基礎(chǔ)學習器的預測結(jié)果比較圖中黑點表示測試集中的觀測值，虛線表示使用M 5′模型樹作為基礎(chǔ)學習器的預測結(jié)果，實線表示使用分類回歸樹作為基礎(chǔ)學習器的預測結(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.