基于K均值聚類分析的流域洪水實時分類修正

李月玉，周建奕，蔣汝成，周 密(.云南省水利水電勘測設計研究院，昆明 6500；.云南農業大學建筑工程學院，昆明 6500)

隨著社會經濟的不斷發展，及時準確的實時洪水預報對確保水庫安全，有效減輕下游洪水威脅和災害，實現整個流域的防洪安全有重要的現實意義。但是由于觀測資料誤差、流域水文規律變化、流域水文模型簡化等問題，使得實時洪水預報系統存在著較多的誤差[1-3]。傳統的洪水預報誤差修正(如：誤差自回歸法、最小二乘法和卡爾曼濾波算法等[4-7])主要運用當前的水情信息，對歷史上的水情信息利用不夠充分，這不僅減少了信息來源量，也浪費了大量歷史信息，同時對預報的精度亦會有一定影響[8]。近年來，對歷史水文信息的數據挖掘尤其是對洪水分類的研究取得了很大進展[9]。在對預報誤差的分析中，常會發現許多洪水的誤差是十分相似的，例如，對于前期土壤含水量較高，降雨范圍高度集中，降雨強度大大超過平均情況的這類型降雨，如果模型仍按平均情況處理，則會使地面徑流估計偏小，匯集速度過慢，使洪峰估計偏小。而對于前期比較干旱，降雨范圍大，但降雨強度比較溫和的這類降雨則正好相反。

本文結合沿渡河流域，采用最為經典，應用也極為廣泛的K均值聚類分析方法對歷史洪水進行聚類[10,11]，在新安江三水源模型的基礎上對不同聚類的洪水分別進行參數率定，在實時洪水預報誤差修正中根據實時的降雨和洪水信息逐時段判斷實時洪水所屬類別，然后根據判別結果采用對應類的模型參數進行模型計算，利用計算所得流量值對初始預報值進行修正，從而將歷史信息運用到誤差修正中，很大程度上減小了模型參數帶來的誤差，提高了預報精度。

1 研究方法

1.1 K均值聚類分析方法

K均值聚類分析方法通過方差分析來篩選最優的分類數，即定義一個F統計量，平均組間平方和與平均組內平方和之比，數值越大說明該特征的組內關系越緊密，而組間關系越離散，分類相對也就更為合理。將試驗中變化的因素稱為因子，因子在試驗中所取的不同狀態稱為水平，構建數學模型為：

F=[SSA/(K-1)][SSE/(n-K)]

(1)

(2)

(3)

式中：K為水平數；ni為第i個水平下的樣本容量；SSA、SSE分別為組間、組內離差平方和。

1.2 流域洪水實時分類修正方法

由于洪水現象十分復雜，變化頻繁，不同類型的洪水具有不同的產匯流規律，如果采用一組水文預報模型參數對全流域洪水進行預報，必然存在很大的預報誤差。根據這一產生誤差的原因，本文建立了基于K均值聚類分析的流域洪水實時分類修正方法，依據輸入的實時降雨及洪水信息，通過計算洪水特征指標到各類中心點的歐式距離對洪水進行實時快速在線分類，然后根據面臨時刻洪水的所屬類別在模型參數庫中選擇相應類別的模型參數，最后利用選取的水文預報模型參數進行模型計算，并利用計算出的流量值對初始預報值進行修正。實時洪水分類修正流程如圖1所示。

圖1 實時洪水分類修正流程圖Fig.1 Flow chart for classified correction of real-time flood forecasting

2 實例應用

本文以沿渡河流域34場洪水為例進行基于K均值聚類分析的流域洪水實時分類修正方法的研究，其中1981-1986年之間的30場洪水用于分類，1987年的4場洪水用于預報檢驗。

2.1 洪水特征指標的選取

聚類前需選擇洪水指標，選擇原則為：①對洪水發生過程有顯著影響的特征因子。②在洪水發生前能獲取這些指標的數值，這樣方能在實際預報中使用[12]。故本文選擇前3天累積降雨量P3 d、前10天累積降雨量P10 d、雨期最大雨強Imax、雨型系數β、雨期累積降雨量P總、起漲流量Q0共6個洪水特征指標。為消除不同指標間的量綱差異，需對數據進行標準化處理，標準化后的數據見表1，分類結果見表2。

表1 洪水特征指標標準化成果表Tab.1 Standardized characteristic index of selected floods

表2 選取的30場洪水的分類結果Tab.2 Flood clusters results

為進一步驗證分類的合理性，對四類洪水的洪水特征值進行統計，結果見表3。由表可看出，各類洪水在洪峰流量、平均降雨量以及平均徑流深上都有明顯的差異，這說明將歷史洪水分為四類能很好地對洪水特征進行區別。

表3 分類洪水的特征Tab.3 Classification characteristics of four types of flood

2.2 模型參數的分類率定

根據表2的聚類分析結果，結合各類洪水的洪水特征，對每一類洪水分別率定一組新安江模型參數[13]。分類后的模型參數率定結果見表4。表中K、WM、WUM、WLM、B、C均為新安江模型的產流參數，分別為蒸發折算系數、流域平均蓄水容量、流域上層蓄水容量、流域下層蓄水容量、流域蓄水容量分布曲線指數、流域蒸發擴散系數；SM、EX、KI、KG為新安江模型的分水源參數，分別為流域自由水平均蓄水容量、流域自由水分布曲線指數、壤中流出流系數、地下水出流系數；CS、CI、CG、KE、XE均為新安江模型的匯流參數，分別為地面線性水庫匯流系數、壤中流線性水庫匯流系數、地下水線性水庫匯流系數、馬斯京根法河段傳播時間、馬斯京根法流量比重系數。

表4 新安江模型參數分類率定結果Tab.4 Parameter estimation for Xin'anjiang model based on flood classification

2.3 實時洪水分類修正

針對1987年的4場洪水，計算其每個時段各指標到各類中心點(各指標平均值構成的6維向量)的歐式距離，逐時段判斷實際洪水所屬類別,選擇其對應的模型參數進行模型計算，然后利用計算值對初始預報值進行修正，其修正效果評定以及與未經修正的傳統預報結果的分析比較見表5。

表5 沿渡河流域洪水分類修正方案與傳統預報分析比較表Tab.5 Comparison of classified correction method with traditional prediction method for real-time flood forecasting in Yandu River basin

以沿渡河“19870719”和“19870823”兩場洪水為例，說明本文提出的方法進行洪水分類實時修正的詳細過程，兩場洪水初始特征指標值見表6。為了避免因線條過多造成圖形繁瑣不易辨析，在圖2、圖3中僅繪出了依據1個時段、2個時段降雨及洪水信息，采用傳統方法預報的洪水過程，以及分類實時修正后的流量過程。對“19870719”這場洪水進行實時修正預報的具體過程如下。

表6 “19870719”、“19870823”兩場洪水初始特征指標值Tab.6 Initial characteristic index of two selected floods

圖2 “19870719”洪水的實測、修正前后的放大的洪峰流量過程線Fig.2 Comparisons among measured data, peak discharge prediction before and after the correction of Flood 19870719

圖3 “19870823”洪水的實測、修正前后的放大的洪峰流量過程線Fig.3 Comparisons among measured data, peak discharge prediction before and after the correction of Flood 19870823

(1)首先根據本場洪水初始狀態及第1個時段降水及洪水信息判斷實際洪水所屬類別為第三類，然后選擇第三類洪水對應的模型參數進行模型計算；

(2)根據初始狀態及第1、第2個時段的降水及洪水信息判斷實際洪水所屬類別為第一類，則選擇第一類洪水對應的模型參數對第二時段的預報值進行實時修正；

(3)按照前述步驟繼續下去直到本場洪水結束，如圖2(“19870719”)。計算結果顯示“19870719”這場洪水實測洪峰流量值為819(m3/s)，修正前的預報值為1 106(m3/s)，修正后的值為935(m3/s)。按照同樣的方法對“19870823”這場洪水進行實時分類修正，其計算結果顯示其實測洪峰值為556(m3/s)，修正前的預報值為465(m3/s)，修正后的值557(m3/s)，如圖3(“19870823”)所示。通過比較可知本文介紹的實時洪水分類修正方法能有效的提高洪水在洪峰部分的預報精度。

從圖2可知修正前的洪峰流量值要明顯的大于實測的流量值，而從圖3又可以看出洪峰流量的實測值明顯大于修正前的預報值。這是因為傳統的洪水預報方法通過率定一組水文模型參數來尋求流域徑流形成的一般性或平均化規律，未能根據實際降雨和洪水信息實時地考慮洪水情況。因此在遇到特殊大洪水或者小洪水的時候，傳統的預報值往往出現偏小或者偏大的情況。本文針對這一問題，根據實時的雨情和洪水信息，選擇合理的模型參數進行實時修正，提高了流域洪水預報的精度，尤其是對洪峰流量的預報精度的改善更為顯著。因此，利用基于K均值聚類分析進行流域洪水實時分類修正的方法是提高整個流域洪水預報精度的有效方法。

3 結 語

洪水預報的主要目的是為水庫及庫區下游提供防洪調度服務，因此，如何預報得到準確的洪峰流量和峰現時間是實際生產中最為關注的問題[14]。在洪水預報誤差分析中，常會發現許多相同類型的洪水會有相似的誤差特性。例如，臺風或雷暴雨型洪水，都是由于降雨范圍高度集中，降雨強度大大超過平均情況，而模型仍按平均情況處理，自然就會使產流估計偏小，匯集速度過慢，使洪峰估計偏小。那么不同場次的這種類型洪水，引起誤差的因素都是高強度和高集中，具有相似性。本文提出的基于K均值聚類分析進行流域洪水實時分類修正的方法彌補了傳統預報方法歷時信息量利用不足，以及遇到特殊情形時會出現較大誤差的缺點，同時也避免了自回歸法在洪峰附近修正效果不佳的問題，從而提高了流域洪水預報精度，尤其是洪峰流量的預報精度。

□

[1] 葛守西．一般線性匯流模型實時預報方法的初步探討[J]．水利學報, 1985,(4):1-9．

[2] 瞿思敏，包為民．實時洪水預報綜合修正方法初探[J]．水科學進展, 2003,14(2):167-172．

[3] 瞿思敏，包為民，石 朋，等．AR模式誤差修正方程參數抗差估計[J]．河海大學學報, 2003,31(5):497-500．

[4] 何少華．遞推最小二乘與誤差自回歸聯合實時校正方法[J]．水電能源科學，1996,(2):78-83．

[5] Fortescue T R, Kershenbaum L S,Yolsue B E. Implementation of Self-turning Regulators with Variable Forgetting Factors [J]. Automatica, 1981,17(6):831-835.

[6] Fread D L, Ming J, Kalman A. Filter Enhanced Real Time Dynamic Flood Routine Model[C]∥Proceedings of XXV Congress of Int．Assoc. For Hydrau. Res. Tokyo,1993.

[7] 宋星原．河道洪水實時預報方法研究[D]．武漢：武漢水利電力大學, 1995．

[8] 包為民．洪水預報信息利用問題研究與討論[J]. 水文, 2006，26(2):102-104．

[9] 萬新宇，包為民，荊艷東，等．基于主成分分析的洪水相似性研究[J]．水電能源科學, 2007,25(5):36-39．

[10] Kanungo T, Mount D M, Netanyahu N S, et al. A Local Search Approximation algorithm for k-means clustering [J]. Computational Geometry, 2004,28(2/3):89-112.

[11] Elkan C. Using the Triangle Inequality to Accelerate k-means[C]∥ /Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003). Menlo Park: AAAI Press, 2003:147-153.

[12] 包為民．水文預報[M]．3版. 北京：中國水利水電出版社，2006．

[13] 王佩蘭，趙人俊. 新安江模型(三水源)參數的客觀優選方法[J]. 河海大學學報, 1989,17(4):65-69.

[14] 李鴻雁,劉寒冰,苑希民,等. 人工神經網絡峰值識別理論及其在洪水預報中的應用[J]. 水利學報, 2002,6(6):15-20.