基于可解釋機器學習的黃河源區徑流分期組合預報

關鍵詞：中長期徑流預報；分期組合；機器學習；可解釋性；黃河源區

中圖分類號：Ｐ３３８＋ ．２；ＴＶ８８２．１ 文獻標志碼：Ａ ｄｏｉ：１０．３９６９／ ｊ．ｉｓｓｎ．１０００－１３７９．２０２４．０９．００８

引用格式：黃強，尚嘉楠，方偉，等．基于可解釋機器學習的黃河源區徑流分期組合預報［Ｊ］．人民黃河，２０２４，４６（９）：５０－５９．

０引言

黃河源區年徑流量占全流域的３４．１％，是黃河重要的產流區和水源涵養區，素有“黃河水塔”之稱［１－２］ 。黃河流域屬于資源性缺水地區，源區來水的變化對全流域水資源開發、利用和管理影響巨大，加強黃河源區徑流預報研究至關重要。此外，全球升溫對黃河源區在內的冰凍圈、大氣圈、水圈、生物圈等多個圈層產生了顯著影響［３］ ，綠色低碳發展已經成為實現溫控目標的全球共識。黃河源區跨越我國一、二級階梯，地勢落差大，水能資源蘊藏量達２ ０００ 萬ｋＷ。在水能資源富集區規劃建設水風光一體化清潔能源基地，是黃河等我國主要流域綠色低碳轉型的重要舉措［４］ 。流域水電站是改善風、光出力隨機波動性的調節中樞，徑流則是決定其調節能力的主控因素，準確預報徑流情勢可為流域清潔能源高效利用、綠色低碳轉型提供關鍵技術支撐。因此，開展黃河源區徑流預報研究對于流域水資源科學調配和綠色低碳發展具有重要意義。

近年來，基于機器學習的數據驅動模型在徑流預報領域快速發展，并得到廣泛應用［５］ 。Ｌｕ 等［６］ 將ＸＧ?Ｂｏｏｓｔ（ Ｅｘｔｒｅｍｅ Ｇｒａｄｉｅｎｔ Ｂｏｏｓｔｉｎｇ） 模型、ＲＦ （ ＲａｎｄｏｍＦｏｒｅｓｔ）模型、ＡｄａＢｏｏｓｔ（Ａｄａｐｔｉｖｅ Ｂｏｏｓｔｉｎｇ）模型應用于富春江水庫的日出庫流量預測，結果表明ＸＧＢｏｏｓｔ 模型預測精度最高且預報能力比較穩定。Ｈａｎ 等［７］ 將深度學習方法作為預測模型誤差后處理器，校正了俄羅斯河徑流預報結果，有效提升了預測精度，并有助于揭示多種不確定因素對模型預測結果的影響。

雖然機器學習善于模擬高維非線性復雜系統，在徑流預報研究中表現良好，但是其結構復雜且具有“黑箱”屬性，無法真正理解和刻畫水文過程及其物理規律，導致預測結果難以被充分信任［８］ 。機器學習的可解釋性分析旨在揭示模型輸入因子對預報對象的驅動機制，提高機器學習模型決策過程的透明度。Ｌｕｎｄｂｅｒｇ 等［９］ 將博弈論與局部解釋性關聯起來，提出了機器學習模型的ＳＨＡＰ 可解釋性分析框架，闡明了各輸入變量對機器學習決策的影響機理。Ｌｉｕ 等［１０］ 構建了一個月預見期的ＸＧＢｏｏｓｔ 徑流預測模型，并使用ＳＨＡＰ 可解釋性分析方法計算了ＳＨＡＰ 總效應值、主效應值、交互作用值和損失值，量化了ＸＧＢｏｏｓｔ 徑流預報模型中氣象要素、歷史徑流、全球氣候因子等輸入變量對徑流預報結果的貢獻率。Ｊｉｎｇ 等［１１］ 將ＬＳＴＭ 模型應用于漢江上游日徑流預測中，利用ＩＧ（Ｉｎｔｅｇｒａｔｅｄ Ｇｒａｄｉｅｎｔ）可解釋性方法研究了徑流預測的潛在決策機制。

基于上述研究，筆者以黃河源區唐乃亥和瑪曲水文站為研究對象，考慮季節性融雪變化影響，建立中長期徑流分期組合機器學習預報模型，構建機器學習預報模型的可解釋性分析框架，定量識別不同預報因子的貢獻率及多因子的復雜交互作用。

１研究區概況與數據來源

１．１研究區概況

黃河發源于青藏高原的巴顏喀拉山北麓，全長５ ４６４ ｋｍ，流域面積７９．５ 萬ｋｍ^２。其中，唐乃亥站以上為黃河源區（見圖１）。源區河道總長１ ５５３ ｋｍ，面積達１２．１９ 萬ｋｍ^２，占黃河流域總面積的１６．２％。唐乃亥站多年平均來水量２００ 億ｍ^３，年徑流量占全流域的３４．１％［２］ ，源區多年平均徑流量分割統計顯示，降雨、地下水補給（基流）、冰雪及凍土融水分別占年總徑流量的６３．１５％、２６．１８％和９．１７％［１２］ 。黃河源區湖泊數量眾多，水資源豐沛。河段落差達３ ７４５ ｍ，水能資源豐富［１３－１４］ ，是我國九大清潔能源基地之一。

１．２數據來源

研究數據涵蓋了水文、氣象、積雪覆蓋、大尺度環流因子等多源數據，時段跨度為１９６０—２０２０年。具體的數據類型及其來源見表１。

２研究方法

基于積雪覆蓋變化和彈性系數精細劃分徑流的年內預報時段，構建基于ＸＧＢｏｏｓｔ 的分期組合預報模型，并使用ＳＨＡＰ 可解釋性分析方法研究不同預報因子對徑流的驅動機制。技術路線見圖２。

黃河源區年內不同時期的融雪徑流占比差異較大，為提升中長期徑流預報的精度，本研究以融雪水當量彈性系數和流域積雪覆蓋率變化為判斷依據，將年內的徑流預報時段劃分為融雪影響期和非融雪主導期，分別構建徑流預報模型。融雪影響期的劃分標準為：１）融雪水當量彈性系數大于０；２）時段內流域積雪覆蓋率呈下降趨勢；３）流域內積雪覆蓋率不小于年內最大積雪覆蓋率（其中，唐乃亥以上流域為３２％、瑪曲以上流域為３８％）的１０％；４）融雪徑流占比達到１０％。參考文獻［１７］，采用下式估算融雪徑流比。

ＳＷ ＝ＳＷ ／ （SＷ＋Ｐ）×１００％ （４）

式中：ＳＷ為融雪徑流占比，ＳＷ 為融雪水當量。

２．２中長期徑流分期組合預報模型

采用分期組合預報策略，分別對融雪影響期和非融雪主導期開展徑流預報，再將兩個時期的預報結果組合，得到年內所有月份的預報結果。分期組合預報策略見圖３。

備選的預報因子分為３ 類，包括水文氣象因子（Ⅰ類，歷史月平均流量、月累計降水量等８ 種）、大尺度環流因子（Ⅱ類，８ 種大尺度環流因子與太陽黑子指數）和融雪水當量（Ⅲ類）。針對融雪影響期和非融雪主導期，采用點間互信息法和相關系數法依次篩選預報因子的最佳組合及最優滯時。需要注意的是，非融雪主導期的融雪徑流占比低，融雪水當量可不作為該期的主要預報因子。

從“預報因子類型”和“是否進行預報因子篩選”兩個角度出發，設置多種組合方案，在融雪影響期和非融雪主導期分別優選出性能最好的預報模型，具體的方案設置見表２。

基于優選出的融雪影響期和非融雪主導期預報模型，建立中長期徑流分期組合預報模型。此外，為評估分期組合預報策略在研究區的適用性，對比分期組合模型和傳統不分期模型的預測性能，對比方案設置見表３。

上述方案的徑流預報模型均基于ＸＧＢｏｏｓｔ模型構建，通過粒子群算法迭代尋優ＸＧＢｏｏｓｔ模型的超參數。模型的率定期設定為１９６０—１９９９年， 驗證期為２０００—２０１４年，測試期為２０１５—２０２０年。

２．３ＳＨＡＰ 機器學習可解釋性分析框架

為了深入挖掘預報模型的決策機制并提升可解釋性，基于ＸＧＢｏｏｓｔ 模型構建徑流分期組合預報模型，并采用ＳＨＡＰ 方法研究其可解釋性分析框架。

通過計算ＳＨＡＰ 值來衡量單變量預報因子、多因子組合對預報結果的貢獻程度，由此可開展預報結果的歸因分析［１８］ 。對于一個已經訓練好的機器學習預報模型，各預報因子的ＳＨＡＰ 值可由下式計算：

根據ＳＨＡＰ 值、ＳＨＡＰ 主效應值和ＳＨＡＰ交互作用值，從各預報因子對預報結果的貢獻程度、預報結果對預報因子的依賴關系、預報因子間的交互作用等角度對建立的徑流預報e206805ff1c1c201b07b8c74d44b5c705ecba92e9ec79d5b53dfcc459f7cd7ef模型進行可解釋性研究。

３結果與分析

３．１融雪影響期與非融雪主導期

基于水量平衡原理，計算了唐乃亥站和瑪曲站以上流域不同月份降水彈性系數和融雪水當量彈性系數（見表４）。結果顯示兩者均為正值，說明降水、融雪水當量與徑流正相關，彈性系數計算結果具有合理性。融雪水當量彈性系數在唐乃亥站以上流域和瑪曲站以上流域的年內變化過程相似，于４—５ 月達最大值，此時徑流對融雪變化最敏感。如圖４ 所示，融雪徑流占比從３ 月起顯著上升，在４—５ 月達到最大（４０％ ～５０％）。由此可見，融雪水是黃河源區地表徑流的重要組成部分［１９－２０］ 。

如圖５ 所示，根據２．１ 節所述的年內預報時段劃分方法，識別出融雪影響期為３—６ 月，非融雪主導期為７ 月—次年２ 月，唐乃亥站和瑪曲站以上流域的預報時段劃分結果一致。

３．２黃河源區中長期徑流分期組合預報

３．２．１融雪影響期和非融雪主導期預報模型優選

根據表２設定的預報模型方案集，分別在融雪影響期和非融雪主導期優選出預報性能最佳的模型方案。

首先，利用點間互信息法初篩各方案（見表２）的徑流預報因子。圖６ 顯示了各預報因子的點間互信息（ＰＭＩ，Ｐｏｉｎｔｗｉｓｅ Ｍｕｔｕａｌ Ｉｎｆｏｒｍａｔｉｏｎ）值。在融雪影響期和非融雪主導期，ＰＭＩ 值較大的前一個月流量（Ｑ）、氣溫、蒸發等氣象要素與唐乃亥站、瑪曲站流量表現出顯著相關性。此外，在融雪影響期，融雪水當量的ＰＭＩ 值均超過０．４（排名第６），也是重要的預報因子。

分別取兩時期內ＰＭＩ 值較大的１０ 種因子作為對應徑流預報時期的預報因子。采用相關系數法確定各方案初篩得到的預報因子的最優滯時，不同滯時預報因子與徑流的相關系數如圖７ 所示。由圖７ 可見，在融雪影響期和非融雪主導期，滯時為一個月的前期流量、前期降水與當月流量之間的相關系數均為１～１２ 個月中的最大值（大于０．６８），關系密切。在融雪影響期，徑流與前一個月融雪水當量的相關系數達０．７５。根據不同滯時預報因子與流量的相關系數大小，確定了不同模型方案預報因子的最優滯時，見表５，其中，Ｑ表示滯時為１ 個月的前期流量，其余因子采用類似的表達方式。

根據表５中的不同方案預報因子集合，采用粒子群優化算法優化ＸＧＢｏｏｓｔ 的超參數，分別率定融雪影響期和非融雪主導期的徑流預報模型。并基于測試期的納什效率系數（ＮＳＥ）、確定系數（Ｒ^２）、均方根誤差與標準偏差之比（ＲＳＲ）和均方根誤差（ＲＭＳＥ），分別篩選融雪影響期和非融雪主導期最優預測模型，結果見表６。

由表６ 可以發現：在融雪影響期，以水文氣象要素和融雪水當量作為預報因子時（ＴＮＨ－Ｓｗ－３－Ｓ 和ＭＱＳｗ－３－Ｓ 方案），模型預報能力均達到最優；在非融雪主導期，將水文氣象要素作為預報因子時（ＴＮＨ－Ｐｒ－１－Ｓ和ＭＱ－Ｐｒ－１－Ｓ 方案），其模型預報能力在兩站均為最佳。

通過篩選預報因子并以最優滯時作為輸入（對應表６ 中加星號的方案），模型預報精度可得到提升。例如，在唐乃亥站，經預報因子篩選后的模型方案ＴＮＨ－Ｓｗ－３－Ｓ 和ＴＮＨ－Ｐｒ－１－Ｓ 與未篩選方案ＴＮＨ－Ｓｗ－３ 和ＴＮＨ－Ｐｒ－１ 相比，ＮＳＥ 分別提高約１７％和６％、ＲＳＲ 分別降低約４４％和１７％、ＲＭＳＥ 分別減少約４４％和１７％。

由表６ 還可知：ＴＮＨ－Ｓｗ－３－Ｓ、ＴＮＨ－Ｐｒ－１－Ｓ 和ＭＱ－Ｓｗ－３－Ｓ模型方案表現較佳；在瑪曲站非融雪主導期內，雖然ＭＱ－Ｐｒ－１－Ｓ 的Ｒ２略低于ＭＱ－Ｐｒ－２，但其余３ 個指標均更優，可視為最優模型。據此，得到研究區中長期分期組合預報模型方案如下：由方案ＴＮＨ－Ｓｗ－３－Ｓ 和ＴＮＨ－Ｐｒ－１－Ｓ 可得到唐乃亥站性能最優的分期組合模型（記為ＴＮＨ－Ｆ－０－Ｓ）；由方案ＭＱ－Ｓｗ－３－Ｓ和ＭＱ－Ｐｒ－１－Ｓ 得到瑪曲站性能最佳的分期組合模型（記為ＭＱ－Ｆ－０－Ｓ）。

３．２．２分期組合預報模型性能分析

將得到的最優分期組合模型與傳統不分期模型（建模方案見表３）進行對比，評估分期組合預報策略在研究區的適用性及其預報性能提升效果。

在構建不分期模型時，運用點間互信息和相關系數方法篩選各模型方案（表３ 中Ｚ－１－Ｓ、Ｚ－２－Ｓ、Ｚ－３－Ｓ 和Ｚ－４－Ｓ 方案）中預報因子類型及其最優滯時，進而采用ＰＳＯ 優化的ＸＧＢｏｏｓｔ 機器學習算法構建預測模型，４ 種方案不分期模型性能見表７。由表７ 可見，與不分期模型相比，分期預報模型在唐乃亥站和瑪曲站的ＮＳＥ 分別提高３％和１０％、ＲＳＲ 分別降低１０％和１７％、ＲＭＳＥ 分別減少１０％和１７％，兩站徑流預報的準確率得到提高。綜上可知，本研究構建的分期組合預報模型（即表７ 中ＴＮＨ－Ｆ－０－Ｓ 和ＭＱ－Ｆ－０－Ｓ）的預報精度高，整體性能較傳統不分期模型明顯提升，分期組合預報策略在研究區具有良好適用性。

為進一步減小徑流預報誤差，提升分期組合預報模型的精度，采用分位數映射方法校正徑流預報誤差，結果見表８。根據《水文情報預報規范》（ＧＢ／ Ｔ ２２４８２—２００８）對預報精度的分級標準，在校正后的ＴＮＨ－Ｆ－０－Ｓ分期組合模型的確定系數Ｒ^２ 為０．９２６，達到甲等精度；ＭＱ－Ｆ－０－Ｓ 分期組合模型的確定系數Ｒ２為０．８５０，接近甲等精度。

３．３分期組合預報模型的可解釋性分析

分期組合模型預報功能的實現， 需要利用ＸＧＢｏｏｓｔ 機器學習方法擬合預報因子與徑流的復雜響應關系。為了深入理解ＸＧＢｏｏｓｔ 在徑流預報中的決策邏輯和工作機制，提升分期組合預報模型的可理解性和透明度，基于上文優選的最優組合預報模型及其預測結果，計算各預報因子的ＳＨＡＰ 值、ＳＨＡＰ 主效應值和ＳＨＡＰ交互效應值，依次解析組合預報模型的單因子貢獻度、因子依賴關系和多因子間交互作用。

３．３．１單因子貢獻度

ＳＨＡＰ 值可以衡量某一模型預報因子對預報結果的貢獻程度，值越大表示其對預報結果的貢獻越大。計算各預報因子的ＳＨＡＰ 均值，并按其絕對值大小進行排序，如圖８ 所示。由圖８可知，在不同水文站和不同預報時段，盡管預報因子對預報結果的貢獻度略有差異，但排名前列的預報因子類型相同。綜合ＳＨＡＰ 均值來看，預測因子按貢獻由大到小依次為降水、前一個月流量、蒸發、氣溫、相對濕度、融雪水當量、氣壓和風速。

采用ＳＨＡＰ 方法進行預報模型可解釋性分析的優勢還在于能夠可視化模型的決策路徑。圖９ 中每條折線代表預報模型對樣本的決策路徑，顏色表示預報流量的大小（藍色表示小流量預報值，紅色表示大流量預報值）。線條方向的變化指示預報因子的影響效果，向左（ＳＨＡＰ 值＜０） 表示會減小預報值， 向右（ＳＨＡＰ 值＞０）則增大預報值。結果表明，在小流量模型中，預報因子的影響較一致，歸納其典型決策路徑為：Ｑ（減小）、ＰＲＥ（減小）、ＥＶＰ（減小）、ＴＥＭ（增大）、ＲＨＵ（增大）、ＷＩＮ（增大）、ＰＲＳ（減小）、ＳＷ（減小）；但隨著預報流量的增大，其預報因子的決策路徑變得更加復雜。

３．３．２預報值對因子的依賴關系

繪制ＳＨＡＰ部分依賴圖，可以識別預報結果與預報因子間線性或非線性、單調或非單調的復雜依賴關系，如圖１０所示。圖１０中，橫、縱軸分別代表預報因子數值與ＳＨＡＰ 主效應值。可以發現：流量預報值與前期降水整體成單調遞增的依賴關系，與蒸發則成單調遞減的依賴關系，且在較大范圍內以線性依賴關系為主。除單調線性關系外，預報徑流與風速還存在非單調、非線性的復雜依賴關系。此外，預報流量值與前期流量、相對濕度整體同樣成單調遞增的關系，與氣溫、氣壓成非單調、非線性的復雜依賴關系。受篇幅限制，此處不再贅述。

對于成單調依賴關系的預報因子，其依賴關系在各徑流預報時段基本相似，但變化閾值存在一定差異。以降水為例，構建的預報模型在不同時期均顯示出了單調遞增的線性依賴關系，在同一預報時期內預報因子閾值接近，推測與研究區高寒山區下墊面對徑流調節機制的影響有關［２１］ 。例如，在融雪影響期，氣溫較低且冰雪、凍土未完全消融，土壤下滲能力弱，少量降水即可產生徑流，在局部依賴圖中表現為ＳＨＡＰ 主效應值開始增長的閾值較小但之后的響應斜率較大；而在非融雪主導期，氣溫和蒸發強度較大，土壤下滲能力增強，需要更加充足的降水才能產生地表徑流，因此相應的閾值也較大。

３．３．３兩因子的交互作用

計算預報因子兩兩組合的ＳＨＡＰ 交互作用值，揭示多因子交互作用對預報流量的影響。圖１１ 為降水與其他預報因t9/YGQpi2FhawAuHBCgpMQ==子的ＳＨＡＰ 交互作用散點圖（以唐乃亥站非融雪主導期ＴＮＨ－Ｐｒ－１－Ｓ 模型、瑪曲站非融雪主導期ＭＱ－Ｐｒ－１－Ｓ 模型為例）。圖中，縱軸表示ＳＨＡＰ交互作用值，散點顏色表示預報因子（降水）值，橫軸是與降水組合的預報因子。ＳＨＡＰ 交互作用值的絕對值越大，交互作用越強；交互作用值為正表示有增大徑流量的效果，負值則相反。

以ＳＨＡＰ 交互作用絕對值判斷交互作用強度，發現降水與氣溫、降水與蒸發、降水與前期流量的交互作用較強。其中，在唐乃亥站非融雪主導期ＴＮＨ－Ｐｒ－１－Ｓ 模型中，前期降水與氣溫交互作用值達－１００，表明在低氣溫下的降水以積雪為主，導致徑流預報值容易出現減小的情況。在瑪曲站非融雪主導期ＭＱ－Ｐｒ－１－Ｓ模型中，前期降水與蒸發交互作用值約為７０，表明當降水量較大、蒸發較少時，會導致徑流量增大。

同時還發現ＳＨＡＰ 交互作用散點分布具有拖尾式或階躍式的特征。拖尾式的特點為在某一區間內交互作用變化不顯著（ＳＨＡＰ 交互作用值基本保持不變），區間外連續變化顯著（ 散點連續分布）， 如ＴＮＨ－Ｐｒ－１－Ｓ模型中降水與蒸發、降水與風速的組合。階躍式的特點為在某一區間內散點分層（分區）現象明顯，反映出當給定某一降水時交互作用變化對橫坐標預報因子的變化不敏感，此時交互作用強度主要由降水量決定，如ＭＱ－Ｐｒ－１－Ｓ 模型中降水與流量、降水與氣壓組合等。

４結論

以黃河源區唐乃亥站和瑪曲站為研究對象，提出了分期組合預報策略，構建了中長期徑流分期組合預報機器學習模型及其可解釋性分析框架，主要結論如下。

１）黃河源區存在顯著的季節性積融雪過程，基于水量平衡原理并結合積雪覆蓋率、融雪水當量變化，定量劃分年內的徑流預報時段，其中融雪影響期為３ 月至６ 月，非融雪主導期為７ 月至次年２月。

２）與傳統模型相比，采用分期組合預報模型能提高黃河上游唐乃亥站和瑪曲站徑流預報精度。與不分期模型相比，分期組合預報模型在唐乃亥站和瑪曲站的ＮＳＥ分別提高３％和１０％、ＲＳＲ 分別降低１０％和１７％、ＲＭＳＥ 分別降低１０％和１７％，經誤差校正后唐乃亥站預報模型Ｒ^２為０.９２６，瑪曲站預報模型Ｒ^２ 為０.８５０。

３）基于ＳＨＡＰ 機器學習可解釋性分析框架，識別出預報因子對徑流預報結果的貢獻程度從大到小依次為降水、前一個月流量、蒸發、氣溫、相對濕度、融雪水當量等。徑流對不同預報因子的依賴關系復雜，其中降水、前期流量和相對濕度以單調遞增的依賴關系為主，蒸發在一定范圍內成單調遞減的線性依賴關系，不同預報因子之間交互作用具有拖尾式或階躍式的特征。