薛聯(lián)青，周天文，劉遠(yuǎn)洪，楊麗娟

（1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2. 皖江工學(xué)院水利工程學(xué)院，安徽 馬鞍山 243031）

0 引 言

準(zhǔn)確可靠的中長期徑流預(yù)測對流域防洪抗旱、水資源規(guī)劃與管理具有重要意義［1］，但徑流在時間和空間尺度上表現(xiàn)出不確定性［2］，對未來的徑流進(jìn)行精準(zhǔn)預(yù)測具有很大難度。近年來，基于數(shù)據(jù)驅(qū)動的機器學(xué)習(xí)模型因其強泛化性和較好的魯棒性而廣泛地被用于徑流預(yù)測［3，4］。然而，由于徑流序列的非線性和復(fù)雜性，在缺少對原始數(shù)據(jù)預(yù)處理的情況下機器學(xué)習(xí)模型難以識別這些特征，從而給準(zhǔn)確的徑流預(yù)測帶來挑戰(zhàn)［5］。通過選擇合適的數(shù)據(jù)預(yù)處理方法如分解技術(shù)，可以濾除混合噪聲和干擾，從而提高預(yù)測精度［6］。如梁浩等［5］基于EMD（Empirical Mode Decomposition）、EEMD（Ensemble Empirical Mode Decomposition）和WD（Wavelet Decomposition）分解方法構(gòu)建多種混合模型，對渭河流域月徑流進(jìn)行預(yù)測，得出混合模型預(yù)測精度均高于單一模型的結(jié)論。

單一分解技術(shù)往往不能完全處理隨機和不規(guī)則數(shù)據(jù)序列的非平穩(wěn)性［7，8］，近年來，有學(xué)者在單一分解技術(shù)的基礎(chǔ)上提出了兩階段分解，如Chen 等［9］構(gòu)建了基于EEMD-VMD 兩階段分解的SVM 模型，對樂昌峽流域平石水文站年徑流進(jìn)行預(yù)測，提高了預(yù)測精度。但月徑流較年徑流系列表現(xiàn)出更強的非線性和季節(jié)性，EEMD 在模態(tài)混疊和分解效率上仍有不足［10，11］，而STL 能夠正確揭示和呈現(xiàn)數(shù)據(jù)中的基本特征［12］，有更高的分解效率，適用于月徑流序列。此外，以往的研究大多未對模型預(yù)測結(jié)果的成因進(jìn)行解釋［8］，而SHAP 是一種基于博弈論的全局靈敏度分析方法［13］，可用于分析模型各輸入特征對預(yù)測結(jié)果的貢獻(xiàn)，從而提高模型的可靠性和應(yīng)用價值。

研究構(gòu)建了基于兩階段分解（STL-VMD）的機器學(xué)習(xí)月徑流預(yù)測模型，同時構(gòu)建了基于未分解、單一分解（STL、VMD）預(yù)處理的對比模型，以期揭示STL-VMD 對不同機器學(xué)習(xí)模型的改善效果，最后基于SHAP 值對最優(yōu)模型進(jìn)行可解釋性分析以提高模型的可信度。

1 研究方法

1.1 分解方法

1.1.1 基于Loess的季節(jié)趨勢分解

STL 是以局部加權(quán)回歸（Loess）作為平滑方法的時間序列分解方法，具有線性最小二乘回歸的簡單性和非線性回歸方法的適應(yīng)性等優(yōu)點［14］。該方法將原始的時間序列分解為趨勢項、季節(jié)項和殘差項，使用穩(wěn)健估計避免了在趨勢項和季節(jié)項中由異常值引起的誤差，且趨勢項平滑程度和季節(jié)項變化速率都可通過相應(yīng)的參數(shù)進(jìn)行控制，較EEMD 的月徑流序列分解效率更高。

1.1.2 變分模態(tài)分解

變分模態(tài)分解（VMD）是一種自適應(yīng)、完全非遞歸的模態(tài)變分和信號處理方法［15］，在獲取分解分量的過程中通過迭代搜尋變分模型最優(yōu)解來確定每個分量的頻率中心和帶寬，從而能夠自適應(yīng)地實現(xiàn)信號的頻域剖分及各分量的有效分離，具有較強的噪音魯棒性和較弱的端點效應(yīng)，解決了EMD 和EEMD 模態(tài)混疊的問題，有更大的潛力準(zhǔn)確預(yù)測徑流時間序列［11］。

1.2 機器學(xué)習(xí)模型

1.2.1 長短期記憶網(wǎng)絡(luò)

長短期記憶網(wǎng)絡(luò)（LSTM）是一種特殊的遞歸神經(jīng)網(wǎng)絡(luò)（Recurrent Neural Network， RNN）［16］，通過門結(jié)構(gòu)設(shè)計選擇性的遺忘或更新歷史信息，解決了傳統(tǒng)RNN的梯度消失和梯度爆炸問題。模型包含輸入層、隱藏層和輸出層，輸入層需提供數(shù)據(jù)類型為［z，m，n］的三維張量X，其中z為批次大小（batch_size），即輸入數(shù)據(jù)的序列長度；m為時間步長（time_step），即模型依次輸入滑動窗口t-m至t-1時刻的特征序列值（圖1紅框所示），根據(jù)徑流序列的周期性，本文m取值為12，即根據(jù)過去12 個月的特征序列值預(yù)測當(dāng)月徑流值；n為輸入特征種類（feature）。

圖1 LSTM模型輸入數(shù)據(jù)處理過程Fig.1 LSTM model input data processing process

1.2.2 卷積神經(jīng)網(wǎng)絡(luò)

卷積神經(jīng)網(wǎng)絡(luò)（CNN）是一類包含卷積計算且具有深度結(jié)構(gòu)的前饋神經(jīng)網(wǎng)絡(luò)［17］，由多層感知機（Multilayer Perceptron，MLP）演變而來。典型的卷積神經(jīng)網(wǎng)絡(luò)由卷積層、池化層和全連接層構(gòu)成，分別實現(xiàn)對輸入數(shù)據(jù)特征的自動提取、降維和輸出。卷積層通過創(chuàng)建卷積核對輸入進(jìn)行卷積，以生成輸出張量，輸出維度通過參數(shù)filters 的大小確定，當(dāng)卷積層為模型第一層時，需提供與LSTM 模型輸入數(shù)據(jù)類型相同的三維張量X［z，m， n］，m取值及初始輸入特征種類均與LSTM模型相同。

1.2.3 支持向量回歸

支持向量回歸（SVR）是基于統(tǒng)計學(xué)理論的機器學(xué)習(xí)方法［18］，遵循結(jié)構(gòu)風(fēng)險最小化原則，通過將低維樣本數(shù)據(jù)映射到高維特征空間，并對支持向量進(jìn)行回歸擬合，實現(xiàn)預(yù)測。與LSTM 和CNN 模型不同，SVR 模型輸入數(shù)據(jù)X類型為形如（n_samples， n_features）或（n_samples， n_samples）的數(shù)組或者稀疏矩陣，其中n_samples為樣本數(shù)量，n_features為特征數(shù)量。

1.3 組合模型構(gòu)建

基于月徑流序列的數(shù)據(jù)特征及上述各模型的基本原理，本文根據(jù)“分解-預(yù)測-重構(gòu)”的思想，構(gòu)建了基于STL-VMD 的月徑流預(yù)測集成模型，該模型主要包括兩階段分解、輸入特征優(yōu)選、機器學(xué)習(xí)模型預(yù)測和結(jié)果重構(gòu)4 個階段，構(gòu)建流程概述如下：

（1）兩階段分解。

①STL 分解：首先采用STL 方法將原始的徑流序列Q分解為趨勢項QTRE、季節(jié)項QSEA和殘差項QRES，即：

②VMD 分解：采用VMD 方法對隨機性較強的殘差項QRES進(jìn)一步分解以剔除噪聲，得到K個分解模態(tài)VMFi（i=1，2，…，K），即：

（2）輸入特征優(yōu)選。計算滯時為一個月的徑流序列與前期徑流、氣象數(shù)據(jù)、氣候指數(shù)等14種特征的皮爾遜相關(guān)系數(shù)，選擇皮爾遜相關(guān)系數(shù)絕對值大于0.2的特征作為模型輸入特征［9］。

（3）機器學(xué)習(xí)模型預(yù)測。將優(yōu)選后的特征整理成圖1 機器學(xué)習(xí)模型能夠識別的格式，為進(jìn)一步提高預(yù)測精度，通過粒子群優(yōu)化算法（Particle Swarm Optimization， PSO）［19］對模型超參數(shù)進(jìn)行尋優(yōu)，尋優(yōu)步驟如下：

①初始化PSO 參數(shù)，即粒子數(shù)、最大迭代數(shù)和慣性因子等，并設(shè)置待優(yōu)化超參數(shù)搜索范圍，對于LSTM和CNN模型，以定型周期和批次大小為優(yōu)化對象，對于SVR 模型，以懲罰系數(shù)和核函數(shù)系數(shù)為優(yōu)化對象；

②選用MAE為優(yōu)化目標(biāo)函數(shù)，求得各粒子對應(yīng)的適應(yīng)度值，得到初始的個體最優(yōu)適應(yīng)度值、全局最優(yōu)適應(yīng)度值及對應(yīng)的參數(shù)；

③更新粒子的速度和位置；

④令k=k+1，重復(fù)步驟②～③，直到k為最大迭代數(shù)，輸出最終求解結(jié)果。

（4）結(jié)果重構(gòu)。將趨勢項、季節(jié)項和K個分解模態(tài)的預(yù)測結(jié)果線性加和，得到月徑流預(yù)測的最終結(jié)果。

1.4 評價指標(biāo)

選取平均絕對誤差（MAE）、平均絕對百分比誤差（MAPE）和納什效率系數(shù)（NSE）來綜合評價模型預(yù)測精度。

MAE表示預(yù)測值和實測值之間絕對誤差的平均值：

MAPE表示預(yù)測值和實測值之間相對偏離程度：

NSE反映了預(yù)測值和實測值的擬合程度：

式中：n為樣本數(shù)；yi為第i個實測值；為實測值的平均值；xi為第i個預(yù)測值。

1.5 可解釋機器學(xué)習(xí)方法

可解釋性不足是機器學(xué)習(xí)模型使用的限制之一，SHAP 是一種基于博弈論的可解釋機器學(xué)習(xí)方法，通過構(gòu)建不同輸入變量的組合比較模型輸出的平均變化，進(jìn)而量化各特征對模型結(jié)果的具體貢獻(xiàn)，每個輸入變量的SHAP 值是該變量邊際貢獻(xiàn)的加權(quán)平均值：

式中：?i表示輸入變量i的SHAP 值，該值為正（負(fù)）時表明變量i對預(yù)測結(jié)果起提升（降低）作用；n代表輸入變量個數(shù)；N代表輸入變量的完整集合；S表示除去變量i的集合，是N的子集；F（S）表示基于輸入S的預(yù)測結(jié)果。

2 實例應(yīng)用

2.1 研究區(qū)概況

澧水位于湖南省西北部，屬洞庭湖四大水系之一，于小渡口注入西洞庭湖，全長390 km，落差1 439 m，流域面積18 583 km2，地勢西北高東南低，流域概況如圖2 所示。石門站為澧水流域控制站，多年平均月徑流量458.51 m3∕s，最大月徑流量2 786.80 m3∕s。流域雨量豐沛，洪水發(fā)生頻繁，大洪水集中在6-8月，下游與長江支流松滋河交匯的松澧地區(qū)，歷年來是洞庭湖區(qū)洪澇災(zāi)害最為嚴(yán)重的地區(qū)之一，開展澧水流域徑流預(yù)測的研究對流域防洪減災(zāi)具有重要意義。

圖2 澧水流域概況圖Fig.2 Overview of Lishui Basin

2.2 數(shù)據(jù)處理

水文過程與大時空尺度上的大氣環(huán)流過程、氣象條件和海洋條件等遙相關(guān)［20］，本文共搜集了徑流、降雨、蒸發(fā)、北極濤動等14種特征作為模型初始輸入特征，數(shù)據(jù)類型及相關(guān)信息如表1所示，并基于皮爾遜相關(guān)系數(shù)對初始輸入特征進(jìn)行篩選，系數(shù)絕對值大于0.2 表示兩變量存在相關(guān)關(guān)系，計算滯時（模型預(yù)見期）一個月的徑流序列與14 種特征的皮爾遜相關(guān)系數(shù)，結(jié)果見表1，最終Q、M1、M2、M3、M5、C1、C2、C3、C5、C7 十種特征作為機器學(xué)習(xí)模型輸入特征。

表1 澧水流域數(shù)據(jù)列表Tab.1 Data list of Lishui Basin

2.3 對照組建模流程

基于不同的分解方法構(gòu)建了4組對比模型（Ⅰ～Ⅳ），建模流程如圖3 所示。其中Ⅰ組為未分解預(yù)處理組；Ⅱ組為STL 預(yù)處理組；Ⅲ組為VMD 預(yù)處理組，根據(jù)文獻(xiàn)［4］的模態(tài)個數(shù)確定方法，原始的徑流序列被分解為IMF1～I(xiàn)MF8；Ⅳ組為STL-VMD 預(yù)處理組，徑流殘差項QRES被分解為VMF1～VMF8。對于Ⅱ～Ⅳ組，除根據(jù)皮爾遜相關(guān)系數(shù)選擇的輸入特征外，對各分解分量逐一預(yù)測時，該目標(biāo)分量的前期序列被視為潛在輸入特征。各特征按比例劃分為訓(xùn)練期（70%，1960.01-2002.02）、驗證期（15%，2002.03-2011.01）和測試期（15%，2011.02-2019.12）并進(jìn)行歸一化處理后輸入到各機器學(xué)習(xí)模型中，訓(xùn)練期實現(xiàn)對模型的擬合，驗證期判斷模型是否過擬合或欠擬合，測試期用來評價模型泛化性［21］，通過PSO 對各機器學(xué)習(xí)模型超參數(shù)進(jìn)行尋優(yōu)，預(yù)見期為一個月，輸出為徑流值（Ⅰ組）或各分解分量的值（Ⅱ～Ⅳ組）。

圖3 對照組建模流程Fig.3 Modeling process of the control group

2.4 結(jié)果與分析

2.4.1 模型超參數(shù)尋優(yōu)結(jié)果

本文粒子數(shù)和最大迭代次數(shù)分別設(shè)置為30 和200，慣性因子ω設(shè)為0.8，學(xué)習(xí)因子c1和c2設(shè)為2，對于LSTM 和CNN 模型，定型周期和批次大小搜索范圍分別設(shè)為［100，500］和［16，256］，學(xué)習(xí)率通過觀察測試期和驗證期的loss 曲線手動調(diào)整進(jìn)行尋優(yōu)；對于SVR 模型，搜索范圍統(tǒng)一設(shè)為［1×10-3，8］，核函數(shù)類型通過手動調(diào)整進(jìn)行尋優(yōu)；尋優(yōu)結(jié)果見表2（Ⅰ組）。

表2 模型超參數(shù)（Ⅰ組）Tab.2 Model hyperparameters （group I）

2.4.2 徑流預(yù)測結(jié)果

為定量評估各模型整體預(yù)測性能，表3 展示了各模型驗證期和測試期預(yù)測值與實測值的MAE、MAPE 和NSE，從表3 可看出：測試期Ⅰ組MAE介于127.44～134.77 m3∕s，MAPE介于28.64%～34.34%，NSE介于0.63～0.64，Ⅱ組和Ⅲ組預(yù)測精度相似，MAE介于92.40～119.57 m3∕s，MAPE介于20.33%～27.97%，NSE介于0.71～0.82，Ⅳ組MAE介于55.65～111.57 m3∕s，MAPE介于9.61%～25.59%，NSE介于0.72～0.93，可見分解預(yù)處理可明顯提升模型預(yù)測精度，且兩階段分解預(yù)處理效果整體優(yōu)于單一分解；對于LSTM和CNN模型，組內(nèi)指標(biāo)變化幅度小于組間變化幅度，可見分解預(yù)處理造成模型輸入的改變，進(jìn)而對模型精度的提升優(yōu)于模型結(jié)構(gòu)變化的提升；Ⅰ～Ⅳ組中SVR 模型預(yù)測精度均不如LSTM 和CNN 模型，這可能由于徑流序列過長，而SVR模型在短時間序列上表現(xiàn)更好有關(guān)［9］。

表3 各組模型徑流預(yù)測精度評價Tab.3 Evaluation of model runoff prediction accuracy in each group

散點圖可更直觀的展示出預(yù)測值與實測值的擬合情況（如圖4 所示），從圖4（a）中可看出Ⅰ組各模型預(yù)測結(jié)果整體偏小，且這一現(xiàn)象在圖4（a）～（d）中高流量事件的預(yù)測尤為明顯；對比圖4（a）和4（b）可看出STL 預(yù)處理對SVR 模型預(yù)測精度的提升較為顯著，決定系數(shù)R2由0.57 提升至0.78；對比圖4（a）和4（c）可看出VMD 預(yù)處理對CNN 模型預(yù)測精度提升較為顯著，R2由0.65 提升至0.88；對比圖4（a）和4（d）可看出STL-VMD 預(yù)處理對LSTM 和CNN 模型預(yù)測精度均有顯著提升，R2由0.64 和0.65提升至0.94 和0.93，預(yù)測值更收斂于實測值。結(jié)合表3 數(shù)據(jù)可看出不同指標(biāo)得到的模型性能優(yōu)劣程度并不完全一致，如對于SVR 模型，根據(jù)MAE、MAPE和R2計算結(jié)果，最優(yōu)模型為STLSVR 模型，而根據(jù)NSE計算結(jié)果，最優(yōu)模型為STL-VMD-SVR模型，這與不同指標(biāo)的評價角度不同有關(guān)。MAE和MAPE分別反映了預(yù)測值可能的誤差范圍、偏離程度，R2側(cè)重于評估預(yù)測值與觀測值的相關(guān)性，當(dāng)模型預(yù)測結(jié)果整體偏高或偏低時，R2往往不能正確的評價預(yù)測效果，而NSE使用平方值計算預(yù)測值與觀測值之間的差異，對峰值十分敏感卻容易忽略低值的預(yù)測誤差［22］，對比圖4（b）和4（d），可看出STL-VMD-SVR 模型的峰值誤差更小，使得該模型NSE優(yōu)于STL-SVR模型。

圖4 各組模型測試期預(yù)測值和觀測值散點圖Fig.4 Scatter plots of predicted and observed values in each group of models during the testing period

高流量事件的準(zhǔn)確預(yù)測對指導(dǎo)流域防洪具有重要意義。從表4中可看出對于測試期三次高流量事件（2015.06、2016.06、2016.07）的預(yù)測，Ⅰ組的相對誤差介于30%～49%，Ⅱ組介于15%～43%，Ⅲ組介于6%～33%，Ⅳ組介于1%～38%，可見未分解預(yù)處理時，各模型對高流量的預(yù)測結(jié)果普遍較差，而對徑流分解預(yù)處理有助于提高模型對高流量事件的預(yù)測精度，但預(yù)測值偏小這一問題依然存在；STL-VMD-LSTM 模型的預(yù)測結(jié)果整體最優(yōu)，3 次高流量事件的相對誤差分別為12.98%、1.34%和6.61%（2015.06、2016.06和2016.07）。

表4 3次高流量事件預(yù)測值及相對誤差Tab.4 Predicted values and relative errors of three high-flow events

2.4.3 STL-VMD-LSTM 模型的可解釋性分析

綜合各評價指標(biāo)計算結(jié)果及高流量事件的預(yù)測精度，STLVMD-LSTM 模型為最優(yōu)模型。為增強該模型的可解釋性，采用SHAP 可視化模型解釋工具探究各輸入特征對預(yù)測結(jié)果的貢獻(xiàn)程度，由1.3 節(jié)可知該模型徑流序列被分解為趨勢項（TRE）、季節(jié)項（SEA）和8 個模態(tài)分量（V1～V8）并逐一預(yù)測，不同目標(biāo)分量對應(yīng)輸入特征的SHAP 值計算結(jié)果如圖5 所示（V3～V8 分量振幅較低且計算結(jié)果與V2 相似，未對其進(jìn)行展示）?？v坐標(biāo)為各輸入特征，由上到下其貢獻(xiàn)度依次遞減，橫坐標(biāo)為SHAP 值，正（負(fù)）值表示對模型預(yù)測的貢獻(xiàn)為正（負(fù)），絕對值越大，貢獻(xiàn)度越大，點的顏色代表輸入特征的值，紅色代表最大值，藍(lán)色代表最小值。

圖5 STL-VMD-LSTM 模型輸入特征SHAP值Fig.5 Input features SHAP value of STL-VMD-LSTM model

從圖5可看出分解后的各分量對預(yù)測結(jié)果的貢獻(xiàn)均優(yōu)于其他輸入特征，進(jìn)一步驗證了兩階段分解對模型精度提升的促進(jìn)作用；除去目標(biāo)特征歷史序列值外，前期降雨（M1）、NINO1+2（C1）和蒸發(fā)（M2）對徑流趨勢項預(yù)測結(jié)果貢獻(xiàn)較大，其中降雨與SHAP值正相關(guān)，表明高（低）降雨值會增大（減小）預(yù)測結(jié)果，蒸發(fā)對預(yù)測結(jié)果整體貢獻(xiàn)為負(fù)［圖5（a）］，符合客觀物理規(guī)律，表明該解釋結(jié)果具有一定的可信度；前期C1和氣溫（M3）、徑流（Q）對徑流季節(jié)項預(yù)測結(jié)果有較大貢獻(xiàn)［圖5（b）］，前期氣溫、NINO1+2、蒸發(fā)對徑流殘差項的低頻分量貢獻(xiàn)較大［圖5（c）～（d）］；同一因子對徑流不同組成的貢獻(xiàn)作用可能不同，如對于徑流分量V1 的預(yù)測，前期氣溫與SHAP 值正相關(guān)［圖5（c）］，而對季節(jié)項、徑流分量V2 的預(yù)測，前期氣溫與SHAP 值呈負(fù)相關(guān)［圖5（b）、5（d）］。

3 結(jié) 論

為提高月徑流預(yù)測精度，根據(jù)月徑流序列強季節(jié)性和高度非穩(wěn)態(tài)的特征，構(gòu)建了基于兩階段分解（STL-VMD）和機器學(xué)習(xí)的集成模型，研究結(jié)論如下。

（1）兩階段分解預(yù)處理可有效地呈現(xiàn)月徑流序列的趨勢、季節(jié)等本質(zhì)特征，對提升徑流預(yù)測精度起促進(jìn)作用，其中對高流量事件預(yù)測精度的提升尤為明顯，STL-VMD-LSTM 模型為最優(yōu)模型，測試期MAE、MAPE和NSE分別為55.65 m3∕s、9.61%和0.93，較單一的LSTM 模型分別改善了56.33%、66.44% 和49.18%；

（2）基于SHAP 值對STL-VMD-LSTM 模型進(jìn)行了可解釋性分析，分析結(jié)果表明歷史降雨、NINO1+2、蒸發(fā)對徑流趨勢項預(yù)測貢獻(xiàn)較大，歷史NINO1+2、氣溫、Q 對徑流季節(jié)項預(yù)測貢獻(xiàn)較大，歷史氣溫、NINO1+2、蒸發(fā)對徑流殘差項低頻分量預(yù)測貢獻(xiàn)較大。