基于水沙組合分類的黃河中下游水沙變化特點研究

蔡蓉蓉，張紅武，卜海磊，張 宇

（1. 清華大學水沙科學與水利水電工程國家重點實驗室，北京 100084；2. 鄭州清大水利工程技術咨詢有限公司，河南鄭州 450003）

1 研究背景

黃河流域的徑流主要來自內蒙古自治區托克托縣河口鎮以上，泥沙主要來源于河口鎮以下至陜西潼關之間，水沙異源［1］，且黃河歷來“水少沙多，水沙關系不協調”，加上下游洪災頻發［2-3］，下游河段淤積嚴重，成為舉世聞名的“地上懸河”。1960 年代以來，黃河下游洪水漫灘機遇減少，泥沙主要淤積在主槽及嫩灘，形成“二級懸河”［4］。下游水沙條件大為改變后，一些河段“二級懸河”發展迅速［5］，“二級懸河”形勢依然嚴峻［6］，河道治理仍面臨挑戰。河床的沖淤變形與水沙組合（徑流量與輸沙量組合）類型密切相關，研究黃河的水沙組合類型有助于深入了解黃河水沙特點、進一步揭示河床沖淤變形規律，對判斷未來的水沙情勢及制定黃河未來治理方略大有裨益。

黃河水沙變化研究，一直是制定黃河治理開發方案的基礎和難點，也是技術界關注與爭論的焦點，尤其自1980 年代以來，有關部門先后圍繞該問題開展相關項目的科學研究，取得了一系列成果［7］。特別地，有學者利用經驗或理論方法，利用實測水沙資料對黃河的水沙關系展開了有關研究。尹學良認為大水期黃河河槽發生沖刷、小水期發生淤積，提出山東段黃河河槽的沖淤分界流量在1800 m3/s左右［8］。邰淑彩借鑒并改進了劉善均研究水文站典型水沙年的方法，與劉善均［9］認為河流的水量和沙量一般呈正相關關系不同，邰淑彩［10］認為黃河的水量和沙量可呈現出非正相關關系，她仍將水（沙）豐隸屬度及其均值作為分類依據，利用閾值β（水沙大年下限）、a（水沙中年界限）、α（水沙小年上限），對渭河華縣站的水文泥沙年資料進行分類并選取了包括小水大沙年在內的典型水沙年。目前，黃河水沙情勢已發生改變［7］，沙量銳減。實測資料表明，黃河中游潼關站年輸沙量已由年均16億t（1919—1959年）銳減至不足3億t（2000—2018年）。信忠保等學者［11］認為退耕還林還草工程使植被恢復可能是導致相應區域產沙量顯著下降的重要原因。馬麗梅、趙躍中、蔣觀滔等學者［12-14］的研究也給出了類似結論，都表明植被覆蓋面積不斷恢復是黃河沙量銳減的重要原因。由于黃河問題的復雜性與研究的重要性，“十三五”期間，國家重點研發計劃首批設立“黃河流域水沙變化機理與趨勢預測”項目（2016YFC0402400），對黃河水沙關系及水沙變化機理繼續展開深入研究。

近年來，以人工神經網絡為代表的機器學習方法在水科學研究中得到了廣泛應用［15］，有學者已將神經網絡技術應用于降雨、徑流預報及模擬等研究中［16］，這為研究黃河水沙關系提供了新的方法。自組織映射（Self-Organizing Map，SOM）模型作為一種典型的神經網絡模型，具有聚類、模式識別等眾多功能，如伊璇等［17］利用該模型進行了滇池流域的分類和無資料區的徑流模擬。鑒于黃河中游潼關站來沙量基本能反映黃土高原的侵蝕量及產沙強度［18］，同時考慮到由于黃河支流伊洛、沁河的匯入及三門峽水庫與小浪底水庫的攔截作用，潼關站與花園口站相比，年徑流量一般偏小，年輸沙量一般偏大，可以說潼關站水沙條件更為不利。本文選取潼關站為代表站，利用該站的實測水沙資料，基于自組織映射-K均值聚類耦合方法（以下簡稱SOM-K法），對該站的水沙組合進行分類，揭示出不同階段各水沙組合類型的變化規律，為深入了解黃河中下游潼關以下河段水沙變化特點提供參考。

2 研究方法

2.1 SOM模型Kohonen參考人腦神經元的特性，于1982年提出SOM模型［19］。SOM是一種無監督學習模型，可將復雜的高維輸入數據通過非線性映射反映到低維（一般是二維）網格，并保留數據的拓撲關系，從而對輸入數據進行有效聚類［19］。SOM 模型由輸入層和可代表網格拓撲結構的輸出層組成，輸入（出）層由存放輸入（出）向量的神經元組成，通常情況下，輸出神經元以二維網格的形式排列，每一個輸入神經元與每一個輸出神經元間通過權值向量連接（圖1）［20］。

圖1 SOM模型結構

在利用SOM 模型進行聚類分析前，需明確距離度量方法，數據初始化、數據訓練所采用的算法類型，確定SOM 模型的網絡尺寸、網格及映射類型、鄰域類型、訓練次數等參數。實踐中一般使用歐氏距離作為距離度量方法，采用隨機初始化方法（random initialization）進行數據初始化［21］。批量訓練（batch training）算法中每次將所有數據向量輸入模型，順序訓練（sequential training）算法中每次僅輸入一個數據向量［22］，前者可節省計算時間，更適用于實際［21］。輸出層神經元數量確定方法沒有統一的標準：Melssen 等［23］認為，輸出層神經元數量應大于所期望聚類組數的2 倍并遠小于輸入數據的樣本數；Lin等［24］認為輸出層神經元數量應盡可能大以確保得到最大的分類；Abrahart等［25］及Parasuraman 等［26］采用反復試驗（trial and error）的方法選取輸出層神經元數量；亦可根據樣本數量n采用公式m=5 n 確定輸出層神經元數量m［22］，并根據輸入數據協方差矩陣最大兩個特征值的比例確定網格的邊長比［27］。為便于可視化，網格類型（grid shape）一般采用六邊形（hexagonal），映射類型（map shape）一般定為片狀（sheet）［27］。在水資源領域使用SOM 模型時，一般選擇高斯（Gaussian）鄰域函數［15］。使用順序訓練算法時，為使結果收斂，訓練次數應至少是輸出層神經元數量的500倍［27］。本文利用mean quantization error（QE）與topographic error（TE）度量SOM模型的聚類效果，QE代表輸入向量與獲勝神經元的平均距離，TE代表數據中第一獲勝神經元與第二獲勝神經元不相鄰的比例。QE和TE越小，說明SOM模型的運算結果越好，但對于小尺寸網格TE并非決定性因素［28］。

2.2 K均值聚類法K均值聚類法是一種尋找輸入數據集中未知子類的無監督學習方法，旨在將輸入數據劃分到指定數量的類中，使類內差異盡可能小而類間差異盡可能大［29-30］。

K均值聚類法的算法步驟如下：（1）確定想要得到的輸入數據集的子類數K；（2）為每個輸入向量隨機分配一個1至K的編號作為該輸入向量的初始子類；（3）將各子類中輸入向量的均值向量作為相應子類的類中心；（4）計算每個輸入向量與各子類類中心的距離，將每個輸入向量分配到距離其最近的類中心所在的子類中；（5）計算調整后新類的類中心，重復步驟（4），直至迭代結束或者子類類中心沒有任何變化［30］。K均值聚類法得到的可能是局部最優解，所得結果依賴于步驟（2）中輸入向量的類別初始化情況，因此，必須多次運行步驟（2）—（5），選擇一個最優的情況［30］。

2.3 SOM-K法SOM模型具有自組織性，算法簡單，但SOM模型網絡收斂時間過長［31］。K均值聚類法效率較高，但聚類結果依賴于所選取的初始聚類中心，容易陷入局部最優情況造成聚類效果不佳［31-32］。K均值聚類法處理噪聲數據的能力遜于SOM模型［33］。黃河徑流量與輸沙量數據變幅較大，變化規律十分復雜，為了取得較好的聚類效果，本文使用自組織映射-K均值聚類耦合方法，即SOM-K法進行分析：將原始數據（矩陣尺寸n×b，n代表樣本數量，b代表特征個數，考慮年徑流量及年輸沙量時b=2）輸入SOM模型得到輸出結果，該輸出結果（矩陣尺寸a×b，a為SOM模型輸出層神經元數量）表示SOM模型對于原始數據矩陣的學習結果，將其作為K均值聚類法的輸入數據進行訓練，最終得到聚類結果。本文利用MATLAB SOM Toolbox 工具箱［34］實現SOM 模型，使用MATLAB R2017a 中的kmeans函數實現K均值聚類。

3 研究區域與數據

黃河發源于青藏高原巴顏喀拉山北麓，在山東省東營市墾利區流入渤海，全長約5464 km，流域面積約75萬km2［2］（圖2）。以內蒙古自治區托克托縣河口鎮與河南省滎陽市廣武鎮桃花峪為分界點將黃河分為上中下游三段。黃河上游長約3472 km，流域面積約占總流域面積的51.3%；中游長約1224 km，流域面積約占總流域面積的45.7%；下游長約768 km，流域面積約占總流域面積的3.0%［2］。

圖2 黃河流域示意圖

本文利用1919—2018 年及1987—2018 年黃河中游潼關站的年徑流量及年輸沙量數據（均為日歷年數據，其中潼關2018年輸沙量包括1億多噸萬家寨等水庫排沙量，取為3.7億t），分別構造水沙數據矩陣，矩陣維度分別為100行（即100年）2列（年徑流量，年輸沙量）及32行2列，基于SOM-K法進行年尺度水沙組合類型分類。由于潼關站早期實測數據匱乏，1919—1951年采用陜縣站數據。年徑流量及年輸沙量數據單位不同，數量級上存在較大差異，故本文使用標準化數據進行分析，數據標準化公式為：

式中：Xst為標準化后的數據；為數據的平均值；σ為數據的標準差。

4 研究結果

4.1 1919—2018年潼關站水沙組合分類將歐氏距離作為距離度量方法，選取隨機初始化方法、批量訓練算法、六邊形網格形狀、片狀映射結構、高斯鄰域函數配置SOM模型，訓練次數選為輸出層神經元數量的500倍。綜合多方面考慮，結合Melssen等提出的方法與反復試驗法確定輸出層神經元數量。根據前人經驗，黃河的水沙組合可分為沙多水多、沙多水中、沙多水少、沙中水多、沙中水中、沙中水少、沙少水多、沙少水中及沙少水少9類，故輸出層神經元數量宜大于18；由于輸入數據樣本數較少，且“遠小于”的概念較為模糊，本文認為輸出層神經元數量上限取為樣本數的1/2 即可。將輸出層最大邊長定為10，結合上述兩條原則，訓練1919—2018年百年數據時共選取39組輸出層神經元網格尺寸（表1）。每一種尺寸的SOM模型均進行500次訓練，選擇QE最小的一次結果作為各尺寸的最終結果。通過比較不同尺寸SOM模型的結果，將輸出層神經元網格尺寸選為3×10，輸出層神經元數量為30，該尺寸下500次訓練的部分QE結果如表2所示（QE的最小值加粗表示）。

表1 輸出層神經元網格尺寸（39組）

表2 部分試驗QE結果（輸出層神經元網格尺寸：3×10）

SOM模型的可視化聚類結果如圖3所示，W代表徑流量，Ws代表輸沙量，TG代表潼關站，兩個3×10六邊形拓撲結構代表兩個變量的標準化值分布，標準化值越大，六邊形顏色越深。采用K均值聚類法將SOM模型的結果分為4大類，使用輪廓值（silhouette value）［35］評價聚類效果。輪廓值在-1到1之間，如果大多數輸入向量擁有較大的輪廓值，可以認為聚類結果是合理的，如果大多數輸入向量擁有較小或者為負數的輪廓值，表明聚類不合理，分得的子類數目可能過多或過少［35］。本次聚類效果如圖4所示，大多數輸入向量擁有較大的輪廓值，將輸入數據集分為4個子類是合理的。

圖3 SOM模型可視化聚類結果（1919—2018年數據）

圖4 聚類效果圖（1919—2018年數據）

潼關站4個子類的年徑流量及年輸沙量分布如圖5所示，依據各子類的年徑流量均值及年輸沙量均值將4個子類依次命名為沙少水中類型、沙多水多類型、沙多水中類型及沙中水中類型，各水沙組合的年徑流量均值及年輸沙量均值如表3 所示。值得注意的是，水沙組合類型中的“多”、“中”，“少”并非絕對的概念，而是與研究時期有關的相對值。

圖5 潼關站4個子類的年徑流量及年輸沙量分布（1919—2018年數據）

表3 潼關站4種類型的水沙均值及出現頻率（1919—2018年數據）

4種類型的出現頻率如表3所示，沙少水中類型出現頻率最大，為42%，沙多水多類型出現頻率仍有20%。周恩來總理曾指出“水土保持是根治開發黃河的基礎”［36］，據《黃河水土保持志》記載，黃河流域的水土保持群眾運動在1958—1960年進入高潮［37］，張含英［38］認為1960年是水土保持工作中具有決定意義的一年，劉瑞龍［39］在這一年提出必須大抓黃河中上游的水土保持工作，可以說1960年是水土保持工作達到高潮并具有重要意義的一年。劉家峽及龍羊峽水庫分別建成于1968年10月及1986年10月［40］，龍劉水庫的修建使黃河干流的水沙條件發生不小的改變［41］。1998 年特大洪水后朱镕基總理提出要“封山植樹、退耕還林”［42］，國務院發出緊急通知，2000 年前要對毀林開墾的林地實現全部還林［43］。劉家峽水庫為不完全年調節水庫，龍羊峽水庫雖為多年調節水庫改變了黃河上游的年際水量［44］，但龍劉水庫聯合運用主要影響徑流年內分配，且黃河沙量主要來源于中游，潼關站年水沙量受龍劉水庫影響較小，故僅將1960 年及2000 年作為時間節點將1919—2018 年分為3 個階段：第1 階段為水土保持運動達到高潮前（1919—1960 年）；第2 階段為水土保持運動達到高潮后至實施退耕還林前（1961—1999年）；第3階段為實施退耕還林后（2000—2018年）。如表4所示，沙少水中類型主要在實施退耕還林后的時期出現，在3 個階段中的出現頻率相對增加，這與現有認識相符，說明本文方法可將水沙類型正確分類；沙多水多及沙多水中類型出現在實施退耕還林前的時期，在3 個階段中的出現頻率相對減少；沙中水中類型在3 個階段中的出現頻率先增后減。1986 年以來黃河年徑流量及年輸沙量大幅度減少［45］，這一點亦可通過水沙組合分類結果看出：1986 年后潼關站有利于下游河道輸沙的大水類型的出現幾率明顯減少（表5）。考慮到黃河水沙條件在1986年后的顯著變化，需將1986 年后潼關站的年徑流量及輸沙量數據進行單獨分類，以明確黃河水沙條件發生顯著變化后潼關站水沙組合類型的變化。

表4 4種類型在3個階段中的出現頻數（頻率）（1919—2018年數據）

表5 1986年后4種類型分別出現的頻率

4.2 1987—2018年潼關站水沙組合分類將輸出層最大邊長定為10，訓練1987—2018年的數據時共選取14組輸出層神經元網格尺寸（表6）。每一種尺寸的SOM模型均進行500次訓練，選擇QE最小的一次結果作為各尺寸的最終結果。通過比較不同尺寸SOM模型的結果，將輸出層神經元網格尺寸選為5×2，輸出層神經元數量為10，該尺寸下500次訓練的部分QE 結果如表7所示（QE 的最小值加粗表示）。

SOM模型的可視化聚類結果如圖6所示。采用K均值聚類法將SOM模型的結果分為2大類，本次聚類效果如圖7所示，大多數輸入向量擁有較大的輪廓值，將輸入數據集分為2個子類是合理的。

表6 輸出層神經元網格尺寸（14組）

表7 部分試驗QE結果（輸出層神經元網格尺寸：5×2）

圖6 SOM模型可視化聚類結果（1987—2018年數據）

圖7 聚類效果圖（1987—2018年數據）

潼關站2個子類的年徑流量及年輸沙量分布如圖8所示。同樣地，依據各子類的年徑流量均值及年輸沙量均值將2個子類分別命名為沙中水中類型及沙少水中類型，各類型的年徑流量均值及年輸沙量均值的出現頻率如表8所示。

采用4.1 節中的處理方法，將1987—2018 年根據時間節點2000 年分為1987—1999 年及2000—2018年2個階段，2種類型在2個階段中的出現頻數（頻率）如表9所示。可以看出，沙中水中類型的出現頻率相對減少，沙少水中類型的出現頻率相對增加。

圖8 潼關站2個子類的年徑流量及年輸沙量分布（1987—2018年數據）

表8 潼關站2種類型的水沙均值及出現頻率（1987—2018年數據）

表9 2種類型在2個階段中的出現頻數（頻率）（1987—2018年數據）

4.3 結果討論比較表4及表9的結果，發現使用1919—2018年百年數據分析時，沙中水中類型的出現頻率先增后減，而使用1987—2018年數據分析時，該類型的出現頻率相對減少。給出百年數據分析結果中4種類型在1919—1960年、1961—1986年、1987—1999年及2000—2018年中的出現頻數（頻率），如表10所示，沙中水中類型在1987—1999年及2000—2018年出現頻率相對減少，與表9結果相符，而在1919—1960年及1961—1986年出現頻率明顯增加，使得從百年尺度上看沙中水中類型的出現頻率先增后減。

表10 4種類型在各時期中的出現頻數（頻率）（1919—2018年數據）

確定水利工程的設計及校核洪水時往往采用洪水頻率分析成果，目前我國選用皮爾遜Ⅲ型曲線作為未知的洪水總體頻率曲線線型［46］，從而根據擬合曲線確定出設計洪水成果，例如將曲線上P=0.1%時的洪峰流量視為尚未發生過的千年一遇洪水洪峰流量。將水沙組合類型中出現頻率最大的一組水沙組合（最高頻水沙組合）下的徑流量及輸沙量視為水沙組合類型特征值，參考洪水頻率分析方法中使用頻率分析成果確定未發生過洪水的思路，認為此特征值可能較接近未來實際情況。4.1節與4.2 節的分析表明，兩種時間尺度（1919—2018 年與1987—2018 年）下出現頻率最大的水沙組合類型（均為沙少水中類型）對應的年徑流量均值及年輸沙量均值分別為249 億m3及5.3 億t 與211.2 億m3及3.1億t，此種差異產生的原因主要是時間尺度不同造成潛在分類閾值不同，即不同的時間尺度上對于黃河水沙豐度的衡量標準也是動態變化的。由王遠見等學者的研究成果［47］，黃河源區1959—2017年的年降雨量有顯著的增加趨勢，由靳少波等學者的研究成果［44］，2007—2016年唐乃亥以上14站平均降雨量較1967—2016年多年均值增加6.2%，鑒于近些年降雨量有所增加的趨勢，綜合4.1節與4.2節結果，為留有一定安全余地，黃河中下游潼關以下河段年徑流量特征值宜選用1919—2018年分析成果，即249億m3。劉曉燕等學者的研究成果［48］表明，1980年代后，龍門、咸陽、張家山、河津及狀頭等水文斷面以上地區單位降雨產沙量（產沙強度）明顯降低，考慮到黃河流域通過水土保持措施使得沙量明顯減少且基本不會大幅度增加的事實，表明未來潼關以下自然狀況下，實際年輸沙量一般不會超過表9 與表10 中2000—2018 年階段沙少水中類型的年輸沙量均值2.4 億t。綜合兩節結果，黃河中下游潼關以下河段年輸沙量特征值選用1987—2018年分析成果即可，即3.1億t。此外，值得注意的是，1987年后沙中水中類型的出現頻率超過30%，為留有一定安全余地，如為在黃河中下游修建水庫時留出充足的淤積庫容，應將沙中水中類型年輸沙量均值7.9 億t（約為8 億t）亦考慮在內，作為黃河干流修建大型工程及下游治理方案論證的沙量特征值。鑒于問題的復雜性，黃河中下游潼關以下河段水量和沙量特征值的確定有待進一步深入研究。

5 結論

本文利用自組織映射-K均值聚類耦合方法對黃河中游潼關水文站1919—2018年的年徑流量及年輸沙量數據進行分析，將潼關站的水沙組合分為沙多水多、沙多水中、沙中水中及沙少水中4 種類型。沙少水中類型的出現頻率最大，為42%，沙多水多類型的出現頻率仍有20%。將研究時段分為水土保持運動達到高潮前（1919—1960年）、水土保持運動達到高潮后至實施退耕還林前（1961—1999年）及實施退耕還林后（2000—2018年）3個階段，發現沙少水中類型主要在實施退耕還林后的時期出現，在3個階段中的出現頻率相對增加；沙多水多及沙多水中類型出現在實施退耕還林前的時期，在3個階段中的出現頻率相對減少；沙中水中類型在3個階段均有出現，出現頻率先增后減。

考慮到近期黃河水沙條件的顯著變化，利用自組織映射-K 均值聚類耦合方法對1986 年后（1987—2018年）潼關站的年徑流量及年輸沙量數據進行分析，將潼關站的水沙組合分為沙中水中及沙少水中2 種類型，沙少水中類型的出現頻率為65.625%。將研究時段分為1987—1999 年及2000—2018年2個階段，發現沙少水中類型出現頻率相對增加，沙中水中類型出現頻率相對減少，與百年分類結果的差異體現出黃河水沙變化的復雜性。

基于自組織映射-K均值聚類耦合方法的潼關站水沙組合分類可較好地反映出實測水沙特性，有助于深入了解黃河中下游水沙變化特點，潼關站水沙組合特征值亦可為黃河中下游水沙變化研究提供參考。