土石壩滲流性態(tài)分析的IAO-XGBoost集成學(xué)習(xí)模型與預(yù)測(cè)結(jié)果解釋

余紅玲，王曉玲，任炳昱，鄭鳴蔚，吳國(guó)華，朱開渲

(天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

1 研究背景

滲流數(shù)值模擬是分析土石壩滲流性態(tài)的一種有效手段，目前有許多不同的數(shù)值模擬方法應(yīng)用于大壩滲流性態(tài)分析研究中[1-3]。雖然各種滲流數(shù)值模擬方法在大壩滲流場(chǎng)計(jì)算分析中發(fā)揮了重要作用，然而，復(fù)雜地質(zhì)條件下的大壩滲流數(shù)值模型對(duì)計(jì)算資源的要求較高，求解時(shí)間較長(zhǎng)，存在計(jì)算效率較低、難以實(shí)時(shí)分析大壩滲流性態(tài)的問題?；跈C(jī)器學(xué)習(xí)方法的代理模型通過尋求輸入變量與輸出變量間的復(fù)雜響應(yīng)關(guān)系，實(shí)現(xiàn)對(duì)復(fù)雜物理過程的近似描述，可代替真實(shí)系統(tǒng)快速給出所求解[4]。常用的代理模型包括徑向基函數(shù)(Radial Basis Function，RBF)、Kriging模型、支持向量回歸(Support Vector Regression，SVR)、多元自適應(yīng)回歸樣條(Multivariate Adaptive Regression Splines，MARS)等單一代理模型以及由多種機(jī)器學(xué)習(xí)算法構(gòu)成的組合代理模型[5]?；跈C(jī)器學(xué)習(xí)算法建立大壩滲流數(shù)值模擬模型的代理模型，能夠減少計(jì)算量、有效提高大壩滲流性態(tài)分析的效率。然而，現(xiàn)有機(jī)器學(xué)習(xí)算法大多將復(fù)雜的滲流數(shù)理機(jī)制作為“黑箱”處理，難以理解模型決策過程，導(dǎo)致模型可信度較低。因此，構(gòu)建一種高效且具有可解釋性的滲流性態(tài)分析模型具有重要的意義。

極限梯度提升(eXtreme Gradient Boosting，XGBoost)[6]是一種基于梯度提升的集成學(xué)習(xí)算法，其不僅具有較快的運(yùn)算速度和較強(qiáng)的泛化能力，還能夠?qū)斎胩卣鞯闹匾冗M(jìn)行排序，具有一定的可解釋性[7]。目前，XGBoost算法在國(guó)內(nèi)外各領(lǐng)域應(yīng)用較為廣泛。在國(guó)內(nèi)，張福浩等[8]基于XGBoost算法構(gòu)建了滑坡隱患點(diǎn)識(shí)別模型，并通過特征重要性分析獲得了影響滑坡發(fā)生的重要因子；張鈞博等[9]將XGBoost算法應(yīng)用于巖爆烈度分級(jí)預(yù)測(cè)研究中，并對(duì)巖石單軸抗壓強(qiáng)度、單軸抗拉強(qiáng)度、洞室圍巖最大切應(yīng)力、巖石彈性變形指數(shù)和巖體完整性系數(shù)等指標(biāo)進(jìn)行了重要性分析；徐韌等[10]采用基于貝葉斯優(yōu)化算法的XGBoost算法對(duì)大壩變形數(shù)據(jù)進(jìn)行了預(yù)測(cè)，并通過特征重要性分析得出溫度分量對(duì)大壩變形的貢獻(xiàn)較大，而水壓和時(shí)效分量的貢獻(xiàn)較??；張書穎等[11]建立了基于XGBoost算法的纖維增強(qiáng)復(fù)合材料加固鋼筋混凝土梁抗彎承載力預(yù)測(cè)模型，并對(duì)模型進(jìn)行了特征重要性分析，確定了影響加固梁承載力的關(guān)鍵因素；丁陽(yáng)陽(yáng)等[12]將XGBoost算法引入煤體結(jié)構(gòu)識(shí)別中，并結(jié)合特征重要性分析方法輸出了模型構(gòu)建中的敏感參數(shù)。在國(guó)外，Li等[13]采用XGBoost算法預(yù)測(cè)了企業(yè)研發(fā)投資的創(chuàng)新績(jī)效，并通過特征重要性排序發(fā)現(xiàn)了影響企業(yè)創(chuàng)新績(jī)效的關(guān)鍵因素；Su等[14]利用XGBoost算法建立了地表沉降的非線性智能預(yù)測(cè)模型，并通過XGBoost算法的特征重要性評(píng)估發(fā)現(xiàn)了決定地表沉降的最主要因素；Liu等[15]建立了基于XGBoost算法的纖維增強(qiáng)復(fù)合材料殘余抗拉強(qiáng)度和模量預(yù)測(cè)模型，并通過特征重要性分析定量評(píng)價(jià)了各屬性參數(shù)對(duì)預(yù)測(cè)結(jié)果的影響；Yan等[16]提出了基于XGBoost算法的裝配式混凝土建筑投資估算模型，并通過特征指標(biāo)重要性排序得出了影響投資估算結(jié)果的重要指標(biāo)。

盡管傳統(tǒng)的XGBoost算法可以通過特征重要性分析，挖掘影響預(yù)測(cè)結(jié)果的關(guān)鍵特征，但是這只能得出特征的重要程度，難以解釋各特征對(duì)預(yù)測(cè)結(jié)果的影響[17]，模型可解釋性有待增強(qiáng)。Shapley加性解釋(SHapley Additive exPlanation，SHAP)是由Lundberg等[18]提出的一種用于增強(qiáng)機(jī)器學(xué)習(xí)模型可解釋性的統(tǒng)一框架，可以對(duì)每一個(gè)樣本的每一個(gè)特征變量計(jì)算出線性可加的貢獻(xiàn)值，從而達(dá)到解釋的效果[19]。與XGBoost算法的特征重要性分析相比，SHAP可以綜合全局和局部?jī)煞矫孢M(jìn)行模型可解釋性分析。首先，從全局出發(fā)，SHAP不僅可以對(duì)樣本特征的重要性進(jìn)行排序，挖掘影響預(yù)測(cè)結(jié)果的關(guān)鍵特征，還可以定性分析樣本特征與預(yù)測(cè)結(jié)果的正負(fù)相關(guān)性；其次，從局部上看，SHAP可以顯示出單個(gè)樣本中各個(gè)特征對(duì)此樣本的預(yù)測(cè)結(jié)果是如何起作用的，顯著提高預(yù)測(cè)結(jié)果的可信度[20]。因此，本研究將XGBoost算法與可解釋機(jī)器學(xué)習(xí)框架SHAP相結(jié)合，以建立具有較強(qiáng)可解釋性的大壩滲流性態(tài)分析模型。

此外，XGBoost集成學(xué)習(xí)算法的模型超參數(shù)較多，超參數(shù)的設(shè)置對(duì)模型預(yù)測(cè)性能具有較大影響?，F(xiàn)有研究大多根據(jù)人工經(jīng)驗(yàn)或網(wǎng)格搜索方法搜尋XGBoost算法的最佳參數(shù)，難以獲得最優(yōu)參數(shù)組合。采用智能優(yōu)化算法對(duì)超參數(shù)進(jìn)行調(diào)整不僅能夠獲得最優(yōu)參數(shù)組合，還可以減少時(shí)間，提升效率[21]。天鷹優(yōu)化器(Aquila Optimizer，AO)是Abualigah等[22]于2021年提出的一種新型智能優(yōu)化算法，具有較強(qiáng)的全局搜索能力和較快的收斂速度。然而，AO算法采用簡(jiǎn)單的隨機(jī)方式對(duì)種群初始化，難以保證初始化種群分布的均勻性和多樣性，并且在開發(fā)階段容易陷入局部最優(yōu)。通過混沌映射產(chǎn)生的混沌序列具備規(guī)律性、遍歷性、隨機(jī)性等特點(diǎn)，可以增加初始種群整體的均勻性和多樣性[23]；天鷹飛行速率是AO算法開發(fā)階段的重要參數(shù)，對(duì)其非線性化對(duì)于增強(qiáng)算法的局部搜索能力具有重要作用。因此，本研究提出基于混沌理論和非線性飛行速率更新策略改進(jìn)的天鷹優(yōu)化(Improved Aquila Optimization，IAO)算法，對(duì)XGBoost集成學(xué)習(xí)算法的超參數(shù)進(jìn)行自適應(yīng)調(diào)優(yōu)，以提高XGBoost算法的預(yù)測(cè)精度。

綜上所述，本研究在可解釋機(jī)器學(xué)習(xí)框架SHAP下，提出一種基于IAO-XGBoost集成學(xué)習(xí)模型的土石壩滲流性態(tài)分析方法，有效解決現(xiàn)有大壩滲流數(shù)值模擬方法計(jì)算效率較低、難以實(shí)時(shí)分析大壩滲流性態(tài)，而基于機(jī)器學(xué)習(xí)方法的代理模型可解釋性較差等不足，從而為大壩運(yùn)行管理人員提供準(zhǔn)確和可靠的大壩滲流性態(tài)分析結(jié)果。

2 研究框架

所提模型的研究框架如圖1所示，主要包括數(shù)據(jù)集生成、基于IAO-XGBoost集成學(xué)習(xí)的預(yù)測(cè)模型建立、基于SHAP理論的預(yù)測(cè)結(jié)果解釋和案例研究四部分。

(1)數(shù)據(jù)集生成。針對(duì)現(xiàn)有復(fù)雜壩基地質(zhì)條件下的土石壩建模方法通常采用先建模后裁剪修型的方式，導(dǎo)致建模過程中涉及大量主觀經(jīng)驗(yàn)判斷，容易影響建模準(zhǔn)確性和建模效率的問題，采用基于有向地層分界面布爾邏輯序列(Boolean Logic Sequence of Oriented Geological Interfaces，BLSOGI)的多地質(zhì)體自動(dòng)建模方法[24]，自動(dòng)建立復(fù)雜壩基地質(zhì)條件下的大壩模型；然后，采用拉丁超立方抽樣(Latin Hypercube Sampling，LHS)方法從上下游水位和各地層滲透系數(shù)的取值范圍中抽取若干樣本點(diǎn)，并將這些樣本點(diǎn)輸入基于CFD技術(shù)的土石壩滲流數(shù)值模型中，通過模擬計(jì)算獲得滲流效應(yīng)量的響應(yīng)值，輸入的樣本點(diǎn)與輸出的滲流效應(yīng)量響應(yīng)值則構(gòu)成大壩滲流性態(tài)指標(biāo)預(yù)測(cè)模型的數(shù)據(jù)集；將數(shù)據(jù)集劃分為訓(xùn)練集和測(cè)試集，其中的訓(xùn)練集用于訓(xùn)練預(yù)測(cè)模型，而測(cè)試集用于檢測(cè)預(yù)測(cè)模型的預(yù)測(cè)性能。

(2)基于IAO-XGBoost集成學(xué)習(xí)的預(yù)測(cè)模型建立。針對(duì)大壩滲流數(shù)值模擬方法存在計(jì)算效率較低、難以實(shí)時(shí)分析大壩滲流性態(tài)的問題，利用XGBoost算法預(yù)測(cè)精度高、運(yùn)算速度快和泛化能力強(qiáng)的優(yōu)勢(shì)，建立大壩滲流性態(tài)指標(biāo)的預(yù)測(cè)模型；進(jìn)一步地，針對(duì)XGBoost算法超參數(shù)的取值問題，在傳統(tǒng)天鷹優(yōu)化算法的基礎(chǔ)上，引入混沌初始化種群和非線性飛行速率更新策略進(jìn)行改進(jìn)，建立改進(jìn)的天鷹優(yōu)化(IAO)算法，并將IAO算法用于自適應(yīng)優(yōu)化XGBoost算法的超參數(shù)n_estimators、learning_rate和max_depth，從而建立基于IAO-XGBoost的大壩滲流性態(tài)指標(biāo)預(yù)測(cè)模型，進(jìn)一步提高大壩滲流性態(tài)指標(biāo)預(yù)測(cè)的精度。

(3)基于SHAP理論的預(yù)測(cè)結(jié)果解釋。針對(duì)基于IAO-XGBoost的大壩滲流性態(tài)指標(biāo)預(yù)測(cè)模型僅能得出預(yù)測(cè)結(jié)果和輸入特征的重要性排序，而難以深入分析各樣本的特征如何影響預(yù)測(cè)結(jié)果的問題，將可解釋機(jī)器學(xué)習(xí)框架SHAP理論與IAO-XGBoost算法相結(jié)合，在對(duì)特征進(jìn)行全局重要性分析，挖掘關(guān)鍵特征的同時(shí)，分析輸入特征對(duì)預(yù)測(cè)結(jié)果的正負(fù)相關(guān)性影響，并解釋特征間交互作用以及單個(gè)樣本的特征對(duì)大壩滲流性態(tài)預(yù)測(cè)結(jié)果的影響，從而提高預(yù)測(cè)模型的可信度。

(4)案例研究。將所提模型應(yīng)用于中國(guó)西南某土石壩工程，通過對(duì)比分析驗(yàn)證所提方法的有效性和優(yōu)越性。

3 土石壩滲流性態(tài)分析的IAO-XGBoost集成學(xué)習(xí)模型

3.1 XGBoost集成學(xué)習(xí)模型XGBoost是由Chen和Guestrin提出的一種基于梯度提升決策樹(Gradient Boosting Decision Tree，GBDT)模型優(yōu)化的集成學(xué)習(xí)方法[6]。XGBoost的基本思想是通過加入新的弱學(xué)習(xí)器擬合前一次訓(xùn)練的殘差，并在訓(xùn)練結(jié)束時(shí)得到每個(gè)樣本的預(yù)測(cè)分?jǐn)?shù)，最后將所有弱學(xué)習(xí)器中的預(yù)測(cè)分?jǐn)?shù)相加而獲得樣本的預(yù)測(cè)值[25]。目前，XGBoost已經(jīng)成功應(yīng)用于眾多領(lǐng)域，且在訓(xùn)練樣本有限、訓(xùn)練時(shí)間較短等場(chǎng)景下具有獨(dú)特優(yōu)勢(shì)[25]。這些特性正好適合于根據(jù)滲流數(shù)值模擬數(shù)據(jù)建立大壩滲流性態(tài)指標(biāo)預(yù)測(cè)模型。基于XGBoost建立大壩滲流性態(tài)指標(biāo)預(yù)測(cè)模型的主要原理如下。

(1)

式中：fk為一棵回歸決策樹；F為所有可能的回歸決策樹的集合；K為回歸決策樹的總數(shù)；fk(xi)為第k棵回歸決策樹對(duì)第i個(gè)滲流性態(tài)指標(biāo)樣本xi的計(jì)算分?jǐn)?shù)。

滲流性態(tài)指標(biāo)預(yù)測(cè)模型訓(xùn)練過程中的目標(biāo)函數(shù)是算法的核心，如下式所示：

(2)

(3)

式中ft(xi)為第t輪迭代時(shí)新加入的回歸決策樹。

對(duì)上式進(jìn)行二階泰勒展開并去掉常數(shù)項(xiàng)constant后，則目標(biāo)函數(shù)表達(dá)為：

(4)

式中g(shù)i和hi分別為損失函數(shù)l的一階和二階導(dǎo)數(shù)。

進(jìn)一步地，將目標(biāo)函數(shù)轉(zhuǎn)換為一個(gè)關(guān)于wj的一元二次方程求最小值的問題，即

(5)

(6)

(7)

(8)

(9)

3.2 基于IAO算法改進(jìn)的XGBoost集成學(xué)習(xí)模型XGBoost超參數(shù)較多，超參數(shù)設(shè)置不準(zhǔn)確將導(dǎo)致大壩滲流性態(tài)指標(biāo)預(yù)測(cè)效率低下、精度降低。然而，超參數(shù)優(yōu)化過程實(shí)質(zhì)上是一個(gè)黑盒函數(shù)優(yōu)化問題，若優(yōu)化參數(shù)過多，則容易使模型冗余，導(dǎo)致計(jì)算復(fù)雜性增加，并影響系統(tǒng)整體性能[21]。因此，本文僅選取對(duì)XGBoost算法預(yù)測(cè)性能影響較大的關(guān)鍵參數(shù)，如n_estimators、max_depth和learning_rate等超參數(shù)進(jìn)行優(yōu)化。其中，n_estimators 為集成算法中弱評(píng)估器的數(shù)量，此參數(shù)值越大，模型的學(xué)習(xí)能力越強(qiáng)，但模型也越容易過擬合；max_depth控制模型中樹的最大深度，其值越大模型越復(fù)雜，且模型容易過擬合；learning_rate參數(shù)控制迭代速率，可以防止模型過擬合[26]。

現(xiàn)有研究大多根據(jù)人工經(jīng)驗(yàn)或網(wǎng)格搜索方法搜尋XGBoost的最佳超參數(shù)。然而，依靠人工經(jīng)驗(yàn)尋找超參數(shù)需要豐富的專業(yè)背景知識(shí)和大量的實(shí)驗(yàn)[21]，而網(wǎng)格搜索方法容易受到維度約束，搜索范圍有限，難以找到最優(yōu)參數(shù)[27]。采用智能優(yōu)化算法對(duì)超參數(shù)進(jìn)行調(diào)整不僅能夠獲得最優(yōu)參數(shù)組合，還可以減少時(shí)間，提升效率。

3.2.1 改進(jìn)的天鷹優(yōu)化算法 天鷹優(yōu)化器(Aquila Optimizer，AO)是Abualigah等[22]于2021年提出的一種新型智能優(yōu)化算法，具有較強(qiáng)的全局搜索能力和較快的收斂速度。天鷹優(yōu)化算法的靈感來(lái)自于天鷹的狩獵行為。天鷹主要有四種狩獵行為，即高空翱翔和垂直俯沖攻擊、等高線飛行和短滑翔攻擊、低空飛行和慢速下降攻擊以及行走攻擊并捕獲獵物。這四種狩獵行為分別屬于擴(kuò)大搜索、縮小搜索、擴(kuò)大開發(fā)和縮小開發(fā)的四個(gè)階段。然而，AO算法采用簡(jiǎn)單的隨機(jī)方式對(duì)種群初始化，難以保證初始化種群分布的均勻性和多樣性，并且在開發(fā)階段容易陷入局部最優(yōu)。針對(duì)上述問題，本文采用混沌理論中搜素速度較快的Tent混沌映射代替?zhèn)鹘y(tǒng)天鷹優(yōu)化算法中的隨機(jī)初始化，以增強(qiáng)種群的均勻性和多樣性；然后，引入非線性飛行速率更新策略替代原來(lái)的線性飛行速率更新策略，以提高天鷹優(yōu)化算法跳出局部最優(yōu)的能力，從而提出改進(jìn)的天鷹優(yōu)化(IAO)算法。

Tent混沌映射是一個(gè)二維混沌映射，具有均勻的分布函數(shù)和良好的相關(guān)性。Tent混沌映射將正態(tài)分布序列轉(zhuǎn)換為混沌序列，能夠增強(qiáng)種群間不同個(gè)體之間的聯(lián)系，增加初始種群整體的隨機(jī)性，從而提升算法前期的全局搜索能力。采用Tent混沌映射生成混沌序列初始化天鷹位置xi，j的過程如下：

首先，由Tent混沌映射生成混沌序列：

(10)

式中：i=1，2，…，N，N為種群中天鷹的總數(shù)；j=1，2，…，D，D為待優(yōu)化超參數(shù)的維數(shù)；xi，j為[0，1]內(nèi)的隨機(jī)值。

然后，將式(10)中的混沌序列映射到解的搜索空間，得到混沌初始化種群：

(11)

式中LBj和UBj分別為第j維待優(yōu)化參數(shù)取值的下界和上界。

天鷹優(yōu)化算法的全局搜索能力較強(qiáng)，但其局部開發(fā)能力較弱，容易陷入局部最優(yōu)解而導(dǎo)致尋優(yōu)精度降低。雖然天鷹優(yōu)化算法在開發(fā)階段中引入了隨機(jī)性較強(qiáng)的萊維(Levy)飛行策略，增強(qiáng)了算法跳出局部最優(yōu)的能力，但是萊維飛行策略在迭代后期容易出現(xiàn)搜索距離過大而導(dǎo)致難以發(fā)現(xiàn)全局最優(yōu)位置的問題。天鷹飛行速率G2是天鷹優(yōu)化算法在縮小開發(fā)階段中的一個(gè)重要參數(shù)，對(duì)于提高算法優(yōu)化性能具有重要意義。然而，天鷹飛行速率G2隨著迭代次數(shù)的增加而呈現(xiàn)線性下降的變化趨勢(shì)。相關(guān)研究[28-30]表明，智能優(yōu)化算法中的非線性參數(shù)對(duì)于增強(qiáng)算法全局、局部搜索能力或提高收斂速度具有重要作用。因此，采用非線性飛行速率更新策略對(duì)天鷹飛行速率G2進(jìn)行改進(jìn)，以彌補(bǔ)算法后期局部搜索能力較弱而使群體陷入早熟的缺陷。改進(jìn)后的G2可表示為G′2，表達(dá)式如下所示：

(12)

式中：t為當(dāng)前迭代次數(shù)；I為最大迭代次數(shù)。

種群混沌初始化增強(qiáng)了原算法的全局搜索能力，非線性飛行速率更新策略增強(qiáng)了原算法的局部搜索能力，因此將兩者結(jié)合起來(lái)可以增強(qiáng)天鷹優(yōu)化算法的整體尋優(yōu)能力。

3.2.2 IAO-XGBoost集成學(xué)習(xí)模型實(shí)現(xiàn)流程 利用IAO算法的尋優(yōu)能力和XGBoost算法的學(xué)習(xí)能力，提出基于IAO算法改進(jìn)的XGBoost(IAO-XGBoost)集成學(xué)習(xí)模型，該模型充分通過IAO算法找到XGBoost訓(xùn)練過程中使其預(yù)測(cè)結(jié)果和模擬值之間誤差最小時(shí)的一組最優(yōu)參數(shù)值，達(dá)到改進(jìn)XGBoost算法預(yù)測(cè)性能的目的。IAO-XGBoost集成學(xué)習(xí)算法的實(shí)現(xiàn)流程如圖2所示。

4 面向大壩滲流性態(tài)智能分析的可解釋機(jī)器學(xué)習(xí)框架SHAP

對(duì)于大壩滲流性態(tài)指標(biāo)預(yù)測(cè)而言，除了準(zhǔn)確性和預(yù)測(cè)速度之外，了解模型為何會(huì)對(duì)不同的樣本輸入特征得出相應(yīng)的預(yù)測(cè)結(jié)果，關(guān)系到模型預(yù)測(cè)結(jié)果的可信度，同樣是大壩運(yùn)行管理人員關(guān)心的問題。機(jī)器學(xué)習(xí)算法在進(jìn)行預(yù)測(cè)時(shí)，其內(nèi)部運(yùn)作過程類似于“黑箱”，難以直觀地被人類理解，因此有必要對(duì)其預(yù)測(cè)結(jié)果進(jìn)行可解釋性分析[31]。在大壩滲流性態(tài)預(yù)測(cè)中，可解釋性意味著找出影響大壩滲流性態(tài)指標(biāo)預(yù)測(cè)的關(guān)鍵特征，闡明模型的決策過程，從而提供大壩運(yùn)行管理人員可理解的滲流性態(tài)指標(biāo)預(yù)測(cè)結(jié)果。

可解釋機(jī)器學(xué)習(xí)是指使機(jī)器學(xué)習(xí)系統(tǒng)的行為和預(yù)測(cè)對(duì)人類可理解的算法或模型[32]。Shapley加性解釋(SHapley Additive exPlanation，SHAP)是由Lundberg等[18]提出的一種用于增強(qiáng)機(jī)器學(xué)習(xí)模型可解釋性的統(tǒng)一框架，可以對(duì)每一個(gè)樣本的每一個(gè)特征變量計(jì)算出線性可加的貢獻(xiàn)值，從而達(dá)到解釋的效果[19]。盡管傳統(tǒng)的XGBoost算法可以進(jìn)行特征重要性分析，挖掘影響預(yù)測(cè)結(jié)果的關(guān)鍵特征，但是這只能得出特征的重要程度，難以解釋該特征是如何影響預(yù)測(cè)結(jié)果的[17]。與XGBoost算法的特征重要性分析相比，SHAP不僅可以對(duì)樣本特征的重要性進(jìn)行排序，挖掘影響預(yù)測(cè)結(jié)果的關(guān)鍵特征，并定性分析樣本特征與預(yù)測(cè)結(jié)果的正負(fù)相關(guān)性，還可以顯示出單個(gè)樣本中各個(gè)特征對(duì)此樣本的預(yù)測(cè)結(jié)果是如何起作用的，顯著提高預(yù)測(cè)結(jié)果的可信度[20]。

可解釋機(jī)器學(xué)習(xí)框架SHAP的核心在于諾貝爾獎(jiǎng)得主Lloyd Shapley提出的Shapley值[33]，該值主要用于解決合作博弈論中的分配均衡問題，即解決如何分配多人合作所創(chuàng)造的總價(jià)值的公平分配問題[34]。在機(jī)器學(xué)習(xí)中，可以把輸入特征看作合作人員，把模型預(yù)測(cè)結(jié)果看作所創(chuàng)造的總價(jià)值，則Shapley值可計(jì)算出各特征對(duì)模型預(yù)測(cè)結(jié)果的貢獻(xiàn)。

假設(shè)第i個(gè)滲流樣本為xi，其第j個(gè)特征值xi，j的Shapley值?xi，j的計(jì)算公式如下所示：

(13)

式中：M是滲流數(shù)據(jù)集中所有特征的集合，其維度為m；S是從M中抽取出來(lái)的子集，其大小為|S|；fxi(S)是只使用特征集合S時(shí)，模型對(duì)滲流樣本xi的預(yù)測(cè)值，當(dāng)S是空集時(shí)，fxi(S)的值稱為基礎(chǔ)值，相當(dāng)于模型的預(yù)測(cè)值在所有滲流樣本上的平均值；fxi(S∪{xi，j})是在特征集合S的基礎(chǔ)之上，添加特征值xi，j時(shí)模型對(duì)滲流樣本xi的預(yù)測(cè)值在所有樣本中的平均值；M{xi，j}表示M中剔除第j個(gè)特征之后的全部特征子集。

基于Shapley值，SHAP對(duì)滲流性態(tài)指標(biāo)預(yù)測(cè)模型的解釋可以用一個(gè)加性模型表達(dá)如下：

(14)

式中：?xi，0為整個(gè)滲流性態(tài)指標(biāo)預(yù)測(cè)模型輸出的基線，即所有訓(xùn)練樣本的預(yù)測(cè)均值；?xi，j為Shapley值，表示第i個(gè)滲流樣本的第j個(gè)特征對(duì)預(yù)測(cè)模型輸出f(xi)的貢獻(xiàn)值；zj為簡(jiǎn)化特征指示系數(shù)，取值為0或1，zj=1表示該特征存在于待解釋樣本xi中；zj=0則表示該特征不存在于xi中。

對(duì)于上式中的Shapley值?xi，j，當(dāng)?xi，j>0時(shí)，說(shuō)明第i個(gè)滲流樣本的第j個(gè)特征提升了預(yù)測(cè)值，對(duì)模型輸出有正向作用；而當(dāng)?xi，j<0時(shí)，說(shuō)明該特征降低了預(yù)測(cè)值，對(duì)模型輸出起反向作用。可見，SHAP算法的作用，是與各個(gè)特征的貢獻(xiàn)值共同將模型的預(yù)測(cè)結(jié)果從基線值推演到最終的模型預(yù)測(cè)值。

5 案例研究

以我國(guó)西南某土石壩工程為研究對(duì)象，該土石壩工程為黏土心墻堆石壩，壩基河床覆蓋層深厚，根據(jù)其物質(zhì)組成、結(jié)構(gòu)特征、成因和分布情況等可自下而上分為4大層7個(gè)亞層，即①漂(塊)卵(碎)礫石層、②-1漂(塊)卵(碎)礫石層(圖3中未切到②-1層)、②-2碎(卵)礫石土層、②-3粉細(xì)砂及粉土層、③-1含漂(塊)卵(碎)礫石土層、③-2礫石砂層和④漂卵礫石層。所選工程的典型剖面地質(zhì)圖如圖3所示，由于③-2層分布范圍較小，故將其與③-1層合并為一層(以③層表示)進(jìn)行研究。

采用基于BLSOGI的多地質(zhì)體自動(dòng)建模方法，建立該土石壩工程的壩基地質(zhì)模型。由于大壩樁號(hào)0+80—樁號(hào)0+251附近的區(qū)域?yàn)橹攸c(diǎn)關(guān)注區(qū)域，故選取該區(qū)域作為滲流分析的計(jì)算區(qū)域，如圖4所示。

圖4 研究區(qū)域地質(zhì)模型圖

在本案例研究中，將壩基防滲帷幕(即圖3中的混凝土防滲墻、覆蓋層帷幕灌漿和基巖帷幕灌漿)的最大水力梯度作為計(jì)算分析的滲流性態(tài)指標(biāo)，將上游水位Hu、下游水位Hd以及壩基覆蓋層①層、②-1層、②-2層、②-3層、③層、④層的滲透系數(shù)(分別表示為K1、K2、K3、K4、K5、K6)作為壩基防滲帷幕最大水力梯度的影響變量。

5.1 IAO-XGBoost集成學(xué)習(xí)模型建立與結(jié)果分析通過對(duì)該工程自2013年5月1日到2021年12月25日的上下游水位監(jiān)測(cè)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，設(shè)置上下游水位的取值范圍。將壩基和壩體各區(qū)域設(shè)置為各向同性的均質(zhì)多孔介質(zhì)，其中壩基各覆蓋層滲透系數(shù)的取值范圍如表1所示，壩體各區(qū)域及壩基其余地層的滲透系數(shù)則根據(jù)現(xiàn)場(chǎng)勘探數(shù)據(jù)和室內(nèi)試驗(yàn)數(shù)據(jù)設(shè)置。

表1 上下游水位和壩基各覆蓋層滲透系數(shù)的取值范圍

根據(jù)表1所示的上下游水位和壩基各覆蓋層滲透系數(shù)的取值范圍，采用最優(yōu)拉丁超立方抽樣方法抽取200組樣本點(diǎn)，并將各樣本點(diǎn)逐個(gè)輸入到滲流數(shù)值模型中進(jìn)行模擬計(jì)算，從而獲得對(duì)應(yīng)的壩基防滲帷幕最大水力梯度值。上下游水位和壩基各覆蓋層滲透系數(shù)的樣本點(diǎn)，與對(duì)應(yīng)的壩基防滲帷幕最大水力梯度模擬值構(gòu)成樣本數(shù)據(jù)集，用于對(duì)壩基防滲帷幕最大水力梯度預(yù)測(cè)模型進(jìn)行訓(xùn)練和測(cè)試。其中，訓(xùn)練集含有160個(gè)樣本，測(cè)試集含有40個(gè)樣本。

采用五折交叉驗(yàn)證方法對(duì)模型進(jìn)行訓(xùn)練，XGBoost算法的超參數(shù)max_depth、learning_rate和n_estimators通過IAO算法進(jìn)行尋優(yōu)，其尋優(yōu)范圍分別為[2，10]、[0.01，1]和[2，100]。在IAO算法中，天鷹種群數(shù)設(shè)置為20，迭代次數(shù)設(shè)置為100次，以均方誤差(Mean Square Error，MSE)作為適應(yīng)度函數(shù)，則超參數(shù)max_depth、learning_rate和n_estimators的尋優(yōu)結(jié)果分別為2、0.2513和54。

在模型訓(xùn)練完成后，以測(cè)試集對(duì)模型預(yù)測(cè)性能進(jìn)行驗(yàn)證。測(cè)試集的預(yù)測(cè)結(jié)果如圖5所示。從圖中可以看出，壩基防滲帷幕最大水力梯度值的預(yù)測(cè)值與模擬值擬合較好。進(jìn)一步計(jì)算預(yù)測(cè)結(jié)果的誤差指標(biāo)，其決定系數(shù)R2、均方根誤差RMSE、平均絕對(duì)誤差MAE和平均絕對(duì)百分比誤差MAPE分別為0.9303、1.3675、1.0890和1.7030%，表明IAO-XGBoost模型的預(yù)測(cè)結(jié)果較為準(zhǔn)確。

圖5 測(cè)試集的預(yù)測(cè)結(jié)果圖

5.2 模型預(yù)測(cè)結(jié)果解釋

5.2.1 樣本特征全局影響分析 在基于IAO-XGBoost模型獲得壩基防滲帷幕最大水力梯度值的預(yù)測(cè)結(jié)果后，可進(jìn)一步對(duì)輸入特征的重要程度進(jìn)行分析，以挖掘影響模型預(yù)測(cè)結(jié)果的關(guān)鍵特征。本文采用特征平均增益值計(jì)算輸入特征的重要性評(píng)分，則各特征的重要性排序如圖6所示。從圖中可以看出，上游水位Hu對(duì)模型預(yù)測(cè)結(jié)果的影響最大，下游水位Hd次之，滲透系數(shù)K4和滲透系數(shù)K6的影響也相對(duì)較高，而滲透系數(shù)K1、K3、K5和K2的影響相對(duì)較小。

雖然IAO-XGBoost算法通過特征重要性分析能夠獲得影響預(yù)測(cè)結(jié)果的關(guān)鍵特征，但是難以判斷各輸入特征對(duì)預(yù)測(cè)結(jié)果的正負(fù)影響。將IAO-XGBoost算法與SHAP理論相結(jié)合可有效解決這一問題，提高模型的可解釋性。

將每個(gè)樣本的每個(gè)特征的Shapley值繪制成散點(diǎn)圖，并用顏色的漸變表示特征值的大小關(guān)系，如圖7所示。在圖7中，每個(gè)點(diǎn)表示各特征對(duì)應(yīng)樣本點(diǎn)的Shapley值，各點(diǎn)的顏色越紅表示特征本身的取值越大，而各點(diǎn)的顏色越藍(lán)，則表示特征本身的取值越小；圖中x軸表示Shapley值，當(dāng)特征的Shapley值大于0時(shí)，表示該特征對(duì)模型預(yù)測(cè)結(jié)果有正影響，當(dāng)特征的Shapley值小于0時(shí)，則表示該特征對(duì)模型預(yù)測(cè)結(jié)果有負(fù)影響；y軸為各特征，各特征根據(jù)其重要性大小從上到下排列。

圖7 特征重要性散點(diǎn)圖

由圖7可知，上游水位Hu、下游水位Hd、滲透系數(shù)K4和滲透系數(shù)K6對(duì)模型預(yù)測(cè)結(jié)果的影響較大。其中，上游水位Hu對(duì)于模型預(yù)測(cè)結(jié)果的影響最大，而且當(dāng)上游水位Hu的值逐步增大時(shí)，其對(duì)應(yīng)的Shapley值大于0并且會(huì)隨著特征值的增大而增大，對(duì)模型預(yù)測(cè)結(jié)果的正影響效果較為顯著；當(dāng)上游水位Hu的值逐漸減小時(shí)，其對(duì)應(yīng)的Shapley值小于0并且也會(huì)隨著特征值的減小而減小，對(duì)模型預(yù)測(cè)結(jié)果的負(fù)影響效果也較強(qiáng)。對(duì)于滲透系數(shù)K6，當(dāng)其特征值增大時(shí)，其對(duì)應(yīng)的Shapley值大于0，但Shapley值增大幅度較小，表示對(duì)模型預(yù)測(cè)結(jié)果的正影響效果不明顯；當(dāng)滲透系數(shù)K6的特征值減小時(shí)，其對(duì)應(yīng)的Shapley值小于0并且對(duì)應(yīng)的Shapley值也在減小，則對(duì)模型預(yù)測(cè)結(jié)果產(chǎn)生負(fù)影響。與其他特征有所不同，下游水位Hd、滲透系數(shù)K4和K1的特征值增大時(shí)，其Shapley值小于0且對(duì)應(yīng)的Shapley值也在減小，對(duì)模型預(yù)測(cè)結(jié)果有負(fù)影響；當(dāng)其特征值減小時(shí)，其Shapley值大于0，則對(duì)模型預(yù)測(cè)結(jié)果有正影響，但該特征對(duì)模型預(yù)測(cè)結(jié)果的正影響沒有負(fù)影響明顯。

SHAP理論也可對(duì)特征重要性進(jìn)行度量，將每個(gè)特征的Shapley值絕對(duì)值的平均值作為該特征的重要性度量標(biāo)準(zhǔn)，其值越大，則表明該特征越為重要，如圖8所示。由圖8可知，上游水位Hu對(duì)應(yīng)的值最大，則表明該特征在模型中最重要，對(duì)模型預(yù)測(cè)結(jié)果的貢獻(xiàn)最大，然后依次是下游水位Hd、滲透系數(shù)K4和K6對(duì)模型預(yù)測(cè)結(jié)果影響較大。由于特征重要性的度量標(biāo)準(zhǔn)不同，故基于SHAP理論的特征重要性排序結(jié)果與基于IAO-XGBoost算法的特征重要性排序結(jié)果(圖6)并不完全相同。但在圖8的特征重要性排序中，僅滲透系數(shù)K5和滲透系數(shù)K3的排序結(jié)果與IAO-XGBoost算法(圖6)的排序結(jié)果不同，其他特征均與圖6中的排序結(jié)果相同。

圖8 基于SHAP理論的特征重要性排序

5.2.2 特征交互影響分析 兩個(gè)輸入特征對(duì)模型預(yù)測(cè)結(jié)果的交互影響可通過圖9所示的SHAP依賴圖可視化。圖9中橫軸表示某一輸入特征的特征值，左縱軸表示該輸入特征對(duì)應(yīng)的Shapley值，右縱軸為以顏色標(biāo)識(shí)的另一輸入特征的特征值，顏色由藍(lán)色至黃色表示輸入特征的特征值由低到高。

圖9 特征對(duì)模型預(yù)測(cè)結(jié)果的交互影響分析

以滲透系數(shù)K3和K4的交互影響作用為例(即圖9(c)與圖9(d))進(jìn)行說(shuō)明。如圖9(c)，對(duì)于滲透系數(shù)K3而言，滲透系數(shù)K3的值越大使得模型預(yù)測(cè)結(jié)果值減小，當(dāng)滲透系數(shù)K3的值小于0.00225 cm/s左右時(shí)，其絕大部分Shapley值大于0，則對(duì)模型預(yù)測(cè)結(jié)果值有正影響；當(dāng)滲透系數(shù)K3的值大于0.00225 cm/s左右時(shí)，其Shapley值小于0，則對(duì)模型預(yù)測(cè)結(jié)果值有負(fù)影響。如圖9(d)，對(duì)于滲透系數(shù)K4而言，滲透系數(shù)K4的值越大同樣使得模型預(yù)測(cè)結(jié)果值減小，當(dāng)滲透系數(shù)K4的值小于0.0005 cm/s左右時(shí)，其絕大部分Shapley值大于0，則對(duì)模型預(yù)測(cè)結(jié)果值有正影響；當(dāng)滲透系數(shù)K4的值大于0.0005 cm/s左右時(shí)，其Shapley值小于0，則對(duì)模型預(yù)測(cè)結(jié)果值有負(fù)影響。綜合滲透系數(shù)K4和K3分析可發(fā)現(xiàn)，當(dāng)滲透系數(shù)K3的值小于0.00225 cm/s時(shí)，滲透系數(shù)K4越大，滲透系數(shù)K3對(duì)模型預(yù)測(cè)結(jié)果的正影響越大，而當(dāng)滲透系數(shù)K3的值大于0.00225 cm/s時(shí)，滲透系數(shù)K4越大，滲透系數(shù)K3對(duì)模型預(yù)測(cè)結(jié)果的負(fù)影響越大；當(dāng)滲透系數(shù)K4的值小于0.0005 cm/s時(shí)，滲透系數(shù)K3越大，滲透系數(shù)K4對(duì)模型預(yù)測(cè)結(jié)果的正影響越大，而當(dāng)滲透系數(shù)K4的值大于0.0005 cm/s時(shí)，滲透系數(shù)K3越大，滲透系數(shù)K4對(duì)模型預(yù)測(cè)結(jié)果的負(fù)影響越大。

5.2.3 單個(gè)樣本特征影響分析 從測(cè)試集中挑選兩個(gè)預(yù)測(cè)樣本，采用SHAP理論解釋輸入特征對(duì)樣本預(yù)測(cè)結(jié)果的影響，如圖10所示。圖中橫坐標(biāo)為預(yù)測(cè)值的大小，縱坐標(biāo)為輸入特征的數(shù)值，紅色表示輸入特征的Shapley值為正，藍(lán)色表示輸入特征的Shapley值為負(fù)；且紅色和藍(lán)色區(qū)域的面積越大，表示Shapley值越大，則對(duì)樣本預(yù)測(cè)結(jié)果的影響也越大。

圖10 樣本預(yù)測(cè)圖

以圖10(a)為例進(jìn)行說(shuō)明。在圖10(a)中，輸入特征對(duì)樣本6的預(yù)測(cè)結(jié)果影響從大到小分別為上游水位Hu、滲透系數(shù)K6、下游水位Hd、滲透系數(shù)K1、滲透系數(shù)K4、滲透系數(shù)K2、滲透系數(shù)K5和滲透系數(shù)K3。其中，上游水位Hu在樣本6中的特征值為1369.762 m，Shapley值為-3.74，表明上游水位Hu對(duì)此樣本的預(yù)測(cè)值產(chǎn)生負(fù)影響；滲透系數(shù)K6在此樣本中的特征值為0.013 cm/s，它的Shapley值為-2.15，表明滲透系數(shù)K6對(duì)此樣本的預(yù)測(cè)值同樣產(chǎn)生負(fù)影響。同理，下游水位Hd的樣本特征值為1303.556 m，Shapley值為+1.01，則下游水位Hd對(duì)樣本6的預(yù)測(cè)值產(chǎn)生正影響。壩基防滲帷幕最大水力梯度預(yù)測(cè)模型輸出的基線為64.854，在圖10(a)中各個(gè)特征的共同作用下，樣本6的最終預(yù)測(cè)值為59.465。

此外，從圖10中兩個(gè)不同樣本的預(yù)測(cè)圖可以看出，同一輸入特征對(duì)不同樣本的預(yù)測(cè)影響存在差異，故對(duì)輸入特征的解釋分析不能太過絕對(duì)，以免產(chǎn)生誤判[19]。

5.3 對(duì)比分析與討論基于測(cè)試集數(shù)據(jù)，采用決定系數(shù)R2、均方根誤差RMSE、平均絕對(duì)誤差MAE和平均絕對(duì)百分比誤差MAPE四種誤差指標(biāo)對(duì)IAO-XGBoost與梯度提升決策樹(Gradient Boosting Decision Tree，GBDT)、隨機(jī)森林(Random Forest，RF)、決策樹(Decision Tree，DT)和支持向量回歸(Support Vector Regression，SVR)四種算法的預(yù)測(cè)性能進(jìn)行對(duì)比分析，如圖11所示。其中，GBDT、RF、DT和SVR的超參數(shù)都通過IAO算法尋優(yōu)而得，故在圖11中分別表示為IAO-GBDT、IAO-RF、IAO-DT和IAO-SVR。

圖11 誤差指標(biāo)計(jì)算結(jié)果

從圖11可以看出，相比于其余四種算法，IAO-XGBoost集成學(xué)習(xí)算法具有最高的R2值，而其RMSE、MAE和MAPE值均較小，表明IAO-XGBoost算法的預(yù)測(cè)精度較高。通過進(jìn)一步計(jì)算可知，與IAO-GBDT、IAO-RF、IAO-DT和IAO-SVR算法相比，IAO-XGBoost算法的預(yù)測(cè)精度分別提高了0.52%、11.64%、37.21%和25.07%。相比于IAO-XGBoost算法，IAO-GBDT和IAO-RF的預(yù)測(cè)精度較低，這是由于XGBoost模型通過引入損失函數(shù)的二階泰勒展開和正則項(xiàng)，優(yōu)化了模型復(fù)雜度，避免了模型過擬合，從而提高了模型預(yù)測(cè)精度。此外，IAO-DT和IAO-SVR的預(yù)測(cè)精度均低于IAO-GBDT和IAO-RF的預(yù)測(cè)精度，這是由于IAO-DT和IAO-SVR屬于單一的機(jī)器學(xué)習(xí)算法，而IAO-GBDT和IAO-RF均屬于集成學(xué)習(xí)算法，具有更好的預(yù)測(cè)性能和泛化能力。

IAO-GBDT和IAO-RF這兩種算法也能夠通過特征重要性分析，得出影響預(yù)測(cè)結(jié)果的關(guān)鍵輸入特征，因而也具有一定的模型可解釋性。IAO-GBDT和IAO-RF的特征重要性排序結(jié)果如圖12所示。將圖12與IAO-XGBoost和SHAP理論的特征重要性排序結(jié)果(如圖6和圖8)進(jìn)行對(duì)比分析，可以看出，這四種算法所得的特征重要性排序并不完全相同，但是對(duì)壩基防滲帷幕最大水力梯度預(yù)測(cè)結(jié)果影響最大的輸入特征均為上游水位Hu，且下游水位Hd、滲透系數(shù)K4和滲透系數(shù)K6均對(duì)模型預(yù)測(cè)結(jié)果具有較大影響。此外，IAO-XGBoost、IAO-GBDT和IAO-RF這三種算法只能用于樣本輸入特征的全局解釋，而SHAP理論還能在全局解釋中發(fā)現(xiàn)輸入特征與壩基防滲帷幕最大水力梯度預(yù)測(cè)值的正負(fù)相關(guān)趨勢(shì)，并且能夠?qū)⑼粋€(gè)輸入特征在不同樣本之間進(jìn)行比較，提供模型預(yù)測(cè)結(jié)果的局部解釋，因而SHAP理論具有更強(qiáng)的模型可解釋性。

圖12 IAO-GBDT和IAO-RF的特征重要性排序

6 結(jié)論

本文建立了土石壩滲流性態(tài)分析的IAO-XGBoost集成學(xué)習(xí)模型，并基于SHAP理論對(duì)預(yù)測(cè)結(jié)果進(jìn)行解釋，最后通過案例研究和對(duì)比分析，驗(yàn)證了所提方法的有效性。所得主要結(jié)論如下：

(1)在采用多地質(zhì)體自動(dòng)建模方法和CFD技術(shù)對(duì)土石壩滲流場(chǎng)進(jìn)行數(shù)值模擬的基礎(chǔ)上，提出了基于IAO-XGBoost集成學(xué)習(xí)的滲流性態(tài)指標(biāo)預(yù)測(cè)方法，采用基于混沌初始化和非線性飛行速率更新策略改進(jìn)的天鷹優(yōu)化(IAO)算法對(duì) XGBoost 集成學(xué)習(xí)算法的超參數(shù)進(jìn)行優(yōu)化，實(shí)現(xiàn)了大壩滲流性態(tài)指標(biāo)的高精度預(yù)測(cè)，并顯著提高了大壩滲流性態(tài)分析的效率。

(2)將可解釋機(jī)器學(xué)習(xí)框架SHAP理論與IAO-XGBoost集成學(xué)習(xí)算法相結(jié)合，挖掘了影響大壩滲流性態(tài)指標(biāo)預(yù)測(cè)結(jié)果的關(guān)鍵特征，并分析了單個(gè)樣本的特征對(duì)預(yù)測(cè)結(jié)果的影響，增強(qiáng)了IAO-XGBoost集成學(xué)習(xí)算法的可解釋性，從而提高了大壩滲流性態(tài)指標(biāo)預(yù)測(cè)結(jié)果的可信度。

(3)案例分析表明，IAO-XGBoost集成學(xué)習(xí)算法具有較高的預(yù)測(cè)精度，相比于IAO-GBDT、IAO-RF、IAO-DT和IAO-SVR算法，其預(yù)測(cè)精度分別提高了0.52%、11.64%、37.21%和25.07%；通過可解釋機(jī)器學(xué)習(xí)框架SHAP理論，得出上游水位Hu和下游水位Hd是影響壩基防滲帷幕最大水力梯度預(yù)測(cè)結(jié)果的關(guān)鍵因素；并且相比于IAO-XGBoost、IAO-GBDT和IAO-RF算法的特征重要性分析方法，SHAP理論不僅能夠挖掘影響壩基帷幕最大水力梯度預(yù)測(cè)結(jié)果的關(guān)鍵特征，并在全局解釋中發(fā)現(xiàn)輸入特征與預(yù)測(cè)結(jié)果的正負(fù)相關(guān)趨勢(shì)以及特征間的交互影響，而且能夠?qū)⑼粋€(gè)輸入特征在不同樣本之間進(jìn)行比較分析，提供模型預(yù)測(cè)結(jié)果的局部解釋，具有更強(qiáng)的模型可解釋性。