基于兩階段分解和可解釋機器學習的中長期徑流預測

薛聯青，周天文，劉遠洪，楊麗娟

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

0 引 言

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

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

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

1 研究方法

1.1 分解方法

1.1.1 基于Loess的季節趨勢分解

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

1.1.2 變分模態分解

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

1.2 機器學習模型

1.2.1 長短期記憶網絡

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

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

1.2.2 卷積神經網絡

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

1.2.3 支持向量回歸

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

1.3 組合模型構建

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

（1）兩階段分解。

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

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

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

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

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

②選用MAE為優化目標函數，求得各粒子對應的適應度值，得到初始的個體最優適應度值、全局最優適應度值及對應的參數；

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

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

（4）結果重構。將趨勢項、季節項和K個分解模態的預測結果線性加和，得到月徑流預測的最終結果。

1.4 評價指標

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

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

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

NSE反映了預測值和實測值的擬合程度：

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

1.5 可解釋機器學習方法

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

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

2 實例應用

2.1 研究區概況

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

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

2.2 數據處理

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

表1 澧水流域數據列表Tab.1 Data list of Lishui Basin

2.3 對照組建模流程

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

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

2.4 結果與分析

2.4.1 模型超參數尋優結果

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

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

2.4.2 徑流預測結果

為定量評估各模型整體預測性能，表3 展示了各模型驗證期和測試期預測值與實測值的MAE、MAPE 和NSE，從表3 可看出：測試期Ⅰ組MAE介于127.44～134.77 m3∕s，MAPE介于28.64%～34.34%，NSE介于0.63～0.64，Ⅱ組和Ⅲ組預測精度相似，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，可見分解預處理可明顯提升模型預測精度，且兩階段分解預處理效果整體優于單一分解；對于LSTM和CNN模型，組內指標變化幅度小于組間變化幅度，可見分解預處理造成模型輸入的改變，進而對模型精度的提升優于模型結構變化的提升；Ⅰ～Ⅳ組中SVR 模型預測精度均不如LSTM 和CNN 模型，這可能由于徑流序列過長，而SVR模型在短時間序列上表現更好有關［9］。

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

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

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

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

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

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

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

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

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

3 結 論

為提高月徑流預測精度，根據月徑流序列強季節性和高度非穩態的特征，構建了基于兩階段分解（STL-VMD）和機器學習的集成模型，研究結論如下。

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

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