李月玉，周建奕，蔣汝成，周 密(.云南省水利水電勘測(cè)設(shè)計(jì)研究院，昆明 6500；.云南農(nóng)業(yè)大學(xué)建筑工程學(xué)院，昆明 6500)

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

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

1 研究方法

1.1 K均值聚類分析方法

K均值聚類分析方法通過方差分析來篩選最優(yōu)的分類數(shù)，即定義一個(gè)F統(tǒng)計(jì)量，平均組間平方和與平均組內(nèi)平方和之比，數(shù)值越大說明該特征的組內(nèi)關(guān)系越緊密，而組間關(guān)系越離散，分類相對(duì)也就更為合理。將試驗(yàn)中變化的因素稱為因子，因子在試驗(yàn)中所取的不同狀態(tài)稱為水平，構(gòu)建數(shù)學(xué)模型為：

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

(1)

(2)

(3)

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

1.2 流域洪水實(shí)時(shí)分類修正方法

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

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

2 實(shí)例應(yīng)用

本文以沿渡河流域34場(chǎng)洪水為例進(jìn)行基于K均值聚類分析的流域洪水實(shí)時(shí)分類修正方法的研究，其中1981-1986年之間的30場(chǎng)洪水用于分類，1987年的4場(chǎng)洪水用于預(yù)報(bào)檢驗(yàn)。

2.1 洪水特征指標(biāo)的選取

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

表1 洪水特征指標(biāo)標(biāo)準(zhǔn)化成果表Tab.1 Standardized characteristic index of selected floods

表2 選取的30場(chǎng)洪水的分類結(jié)果Tab.2 Flood clusters results

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

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

2.2 模型參數(shù)的分類率定

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

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

2.3 實(shí)時(shí)洪水分類修正

針對(duì)1987年的4場(chǎng)洪水，計(jì)算其每個(gè)時(shí)段各指標(biāo)到各類中心點(diǎn)(各指標(biāo)平均值構(gòu)成的6維向量)的歐式距離，逐時(shí)段判斷實(shí)際洪水所屬類別,選擇其對(duì)應(yīng)的模型參數(shù)進(jìn)行模型計(jì)算，然后利用計(jì)算值對(duì)初始預(yù)報(bào)值進(jìn)行修正，其修正效果評(píng)定以及與未經(jīng)修正的傳統(tǒng)預(yù)報(bào)結(jié)果的分析比較見表5。

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

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

表6 “19870719”、“19870823”兩場(chǎng)洪水初始特征指標(biāo)值Tab.6 Initial characteristic index of two selected floods

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

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

(1)首先根據(jù)本場(chǎng)洪水初始狀態(tài)及第1個(gè)時(shí)段降水及洪水信息判斷實(shí)際洪水所屬類別為第三類，然后選擇第三類洪水對(duì)應(yīng)的模型參數(shù)進(jìn)行模型計(jì)算；

(2)根據(jù)初始狀態(tài)及第1、第2個(gè)時(shí)段的降水及洪水信息判斷實(shí)際洪水所屬類別為第一類，則選擇第一類洪水對(duì)應(yīng)的模型參數(shù)對(duì)第二時(shí)段的預(yù)報(bào)值進(jìn)行實(shí)時(shí)修正；

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

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

3 結(jié) 語

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

□

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

[2] 瞿思敏，包為民．實(shí)時(shí)洪水預(yù)報(bào)綜合修正方法初探[J]．水科學(xué)進(jìn)展, 2003,14(2):167-172．

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

[4] 何少華．遞推最小二乘與誤差自回歸聯(lián)合實(shí)時(shí)校正方法[J]．水電能源科學(xué)，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] 宋星原．河道洪水實(shí)時(shí)預(yù)報(bào)方法研究[D]．武漢：武漢水利電力大學(xué), 1995．

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

[9] 萬新宇，包為民，荊艷東，等．基于主成分分析的洪水相似性研究[J]．水電能源科學(xué), 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] 包為民．水文預(yù)報(bào)[M]．3版. 北京：中國(guó)水利水電出版社，2006．

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

[14] 李鴻雁,劉寒冰,苑希民,等. 人工神經(jīng)網(wǎng)絡(luò)峰值識(shí)別理論及其在洪水預(yù)報(bào)中的應(yīng)用[J]. 水利學(xué)報(bào), 2002,6(6):15-20.