劉 恒

(遼寧省水利水電科學研究院有限責任公司，遼寧 沈陽 110003)

0 引 言

單純依靠改變模型參數(shù)來提高不同類型洪水的預報精度，雖然在一定程度上能提高洪水預報精度；但具有一定的局限性，仍然無法避免單個水文模型結(jié)構(gòu)的不確定性對預報結(jié)果的影響。目前，國內(nèi)外已經(jīng)提出很多種水文模型，分別在不同的地區(qū)和流域發(fā)揮重要作用；但是，每種水文模型都是學者對自然界降雨徑流及洪水運動規(guī)律的一種理解和概化，無法全面描述實際的產(chǎn)匯流過程[1]。這也表明每種模型都存在一定的優(yōu)缺點，不能說明哪種水文模型在任何情況都完全優(yōu)于其他模型[2-3]。既然單個模型存在一定的缺陷，那么能不能將幾種模型的預報結(jié)果進行組合，盡可能地減少單一模型結(jié)構(gòu)不確定性對預報結(jié)果的影響，降低模型結(jié)構(gòu)帶來的誤差?；谶@一想法，國內(nèi)外學者提出了多模型組合預報的概念。即，綜合多個模型的預報結(jié)果，提高整體預報精度。目前，對組合洪水預報方法有很多，早期的水文模型綜合研究采用簡單平均法、加權(quán)平均法、模糊數(shù)學法、神經(jīng)網(wǎng)絡模型、貝葉斯平均法等方法。這些方法彌補了單個模型的缺陷，充分發(fā)揮每個模型的優(yōu)勢[4]。為進一步提高預報精度，將同時考慮模型參數(shù)不確定性與模型結(jié)構(gòu)不確定性對預報結(jié)果的影響，本文將在多組概念性水文模型的基礎上，采用神經(jīng)網(wǎng)絡、遺傳算法、貝葉斯、蒙特卡羅等算法，將單模型的分類洪水預報與多模型的組合預報結(jié)合起來，建立分類洪水組合預報模型，綜合考慮各模型的優(yōu)勢，從而提高預報精度。關(guān)于概念性水文模型、基于BP神經(jīng)網(wǎng)絡的洪水分類，以及基于遺傳算法的參數(shù)率定相關(guān)論文比較多，本文在此不對相關(guān)理論進行贅述。

1 基于貝葉斯模型平均法的組合預報

1.1 貝葉斯平均法基本原理

貝葉斯模型平均法(BMA)是由A. E. Rafery等人提出的，利用多模型進行概率預報的統(tǒng)計處理方法[5]，是一種基于貝葉斯理論，將模型本身的不確定性考慮在內(nèi)的統(tǒng)計分析方法。作為一種能夠提高模擬結(jié)果的有效綜合計算方法，貝葉斯模型平均方法是通過綜合幾個模型預報值的后驗分布來推斷預報量的更可靠概率分布分析工具。它不僅能提供一個綜合的預報值，同時可以獲取描述模擬結(jié)果的不確定區(qū)間，為用戶通過更多的參考依據(jù)，被國內(nèi)外學者廣泛應用于眾多研究領域[6-9]。貝葉斯平均法克服了其他預測模型的一系列確定，包括未考慮主觀先驗信息、沒有提取各模型正確的預測信息、沒有考慮模型不確定性等。其基本表達式如下

(1)

式中，yBMA為BMA法的組合預測值；yMk為單個模型Mk的預測值；t為變量；P(Mk|D)是給定數(shù)據(jù)D情況下，模型的后驗概率。BMA法的實質(zhì)：以單個模型的后驗概率為權(quán)重，對單一模型的預測值進行加權(quán)平均得到貝葉斯模型均值。

1.2 基于貝葉斯平均法的組合預報

近年來，BMA方法被應用到洪水預報領域，將多個水文模型的預報結(jié)果進行組合，這對于提高我國水文預報的水平具有重要的現(xiàn)實意義和科學價值[10-11]。實現(xiàn)對各水文模型的預報流量序列進行加權(quán)平均，獲得組合預報值和預報區(qū)間[12]，具體模型構(gòu)建過程如下：

1.2.1概率密度函數(shù)

設Q為預報變量，D=[X，Y]為實測數(shù)據(jù)(X為輸入資料，Y為實測的流量資料)，f=[f1，f2，…，fk]是K個模型預報的集合。根據(jù)貝葉斯模型平均法的理論，BMA的概率要表示為

(2)

首先將模型預報值和實測流量都用Box-Cox轉(zhuǎn)化方法進行正態(tài)轉(zhuǎn)化，可以方便BMA計算。因此，綜合預報量Q的后驗分布實際上是以后驗模型概率p(fk|D)為權(quán)重，對所有模型的后驗分布pk(Q|fk，D)進行加權(quán)的一個平均值。

1.2.2期望與方差計算

BMA平均預報值是單個模型預報值的加權(quán)平均結(jié)果。如果單個模型預報值和實測流量均服從正態(tài)分布，則BMA平均預報值

(3)

(4)

1.2.3正態(tài)轉(zhuǎn)換

由于實測和預報流量序列都存在一定的不確定性，不符合正態(tài)分布的規(guī)律，如果在原始數(shù)據(jù)基礎上直接推求概率密度函數(shù)，將該計算帶來很大的麻煩。因此，首先要將這些數(shù)據(jù)采用一定的轉(zhuǎn)換規(guī)則將其轉(zhuǎn)為服從正態(tài)分布的數(shù)據(jù)。本文采用Box-Cox(博克斯-卡克斯)轉(zhuǎn)換法，即

(5)

式中，odt是為t時刻未轉(zhuǎn)化的原始數(shù)據(jù)；tdt是為轉(zhuǎn)化后的數(shù)據(jù)；λ是為Box-Cox轉(zhuǎn)化系數(shù)。因為EM算法中需要把模型預報值結(jié)合計算，所以對所有的模型預報值和實測流量，都采用同一個轉(zhuǎn)化系統(tǒng)λ。

1.2.4期望最大化算法

(6)

(1)初始化。設Iter=0

(7)

(8)

(2)計算初始似然值。即

(9)

(3)計算隱藏變量。設Iter=Iter+1

(10)

(4)計算權(quán)重

(11)

(5)計算模型預報誤差

(12)

(6)計算似然值

(13)

(7)檢驗收斂性。如果l(θ)(Iter)-l(θ)(Iter-1)小于或等于預先設定的允許誤差，就停止；l(θ)(Iter)-l(θ)(Iter-1)小于或等于預先設定的允許誤差，就停止；否則回到步驟(3)。

2 基于神經(jīng)網(wǎng)絡和貝葉斯平均法的分類組合預報

基于神經(jīng)網(wǎng)絡和貝葉斯平均法的分類多模型組合預報的基本原理是將洪水分類預報與多模型組合預報相結(jié)合。首先，基于BP神經(jīng)網(wǎng)絡模型實現(xiàn)洪水的分類，將洪水劃分為高水、中水、低水三類；然后，針對不同類型洪水采用遺傳算法進行參數(shù)率定，獲取每種類型洪水對應的最優(yōu)或近似參數(shù)組，實現(xiàn)分類洪水預報。在分類洪水預報的基礎上，本文選取新安江、大伙房、遼寧指數(shù)3個概念性水文模型(工程水文學)，使用BMA方法將3模型針對不同類型洪水的預報結(jié)果進行綜合，在綜合之前需要對實測和預報序列進行正態(tài)轉(zhuǎn)換。然后以實測序列隸屬于某個水文模型后驗概率為權(quán)重，對所選單項模型方法的預報變量進行加權(quán)平均，實現(xiàn)了不同水文模型預報的合成及概率預報。將模型本身參數(shù)與結(jié)構(gòu)的不確定性考慮在內(nèi)，最終獲得分類洪水最優(yōu)參數(shù)的多模型組合預報值及概率區(qū)間。分類多模型組合預報包括模型構(gòu)建和模型計算兩個環(huán)節(jié)。

2.1 模型構(gòu)建

(1)洪水分類模型構(gòu)建?；贐P神經(jīng)網(wǎng)絡構(gòu)建洪水分類模型，采用歷史典型洪水對模型結(jié)構(gòu)進行率定與檢驗，最終構(gòu)建出基于BP神經(jīng)網(wǎng)絡的在線洪水分類模型[13]。

(2)不同類型洪水模型參數(shù)率定。基于遺傳算法對不同類型洪水(高水、中水、低水)的不同模型(新安江模型、大伙房模型、遼寧指數(shù)模型)進行參數(shù)率定，從而獲得每種類型洪水對應不同模型參數(shù)。

(3)權(quán)重計算。對實測流量與預報流量正態(tài)轉(zhuǎn)換，采用BMA方法(見圖1)分別對每種類型洪水3個水文模型的預報流量序列進行分析計算，采用期望最大化算法(EM)估計BMA中的未知參數(shù)，從而獲得獲取3種類型洪水對應的3組權(quán)重[14]。

圖1 分類組合預報流程

2.2 模型計算

(1)在線分類。輸入實測降雨及洪水數(shù)據(jù)，基于BP神經(jīng)網(wǎng)絡在線洪水分類模型對實時洪水進行分類，確定該場次洪水所屬類別。

(2)分類預報。根據(jù)洪水類別，分別采用大伙房模型、新安江模型、遼寧指數(shù)模型3個水文模型進行預報，可獲得3組預報流量序列。

(3)組合預報。根據(jù)每種類型洪水對應的權(quán)重組合，對新安江、大伙房、遼寧指數(shù)3種模型的預報流量序列進行組合，從而獲得每種類型洪水的預報均值和預報區(qū)間。

3 基于蒙特卡羅法的不確定性區(qū)間估計

(14)

(3)重復步驟(1)和(2)M次。M是為生成不確定性區(qū)間取樣的總個數(shù)，令M=1 000。BMA在任意時刻t的1 000個樣本值由上述方法取樣得到以后，將他們從小到大排序，BMA的90%預報不確定性區(qū)間就是5%和95%分位數(shù)之間的部分。

4 分類洪水組合預報模型實例驗證

4.1 基礎資料

選取大伙房水庫1960年～2016年的25場歷史典型洪水作為本研究的基礎資料，并將洪水分成3類，包括5場高水、9場中水、11場低水。選擇20場典型洪水資料作為模型的率定期，包括4場高水、7場中水、9場低水。另外5場典型洪水資料進行模型的檢驗和比較，其中包括1場高水、2場中水、2場低水。歷史洪水信息見表1。

4.2 驗證與分析

4.2.1正態(tài)轉(zhuǎn)換

采用P-P圖來判斷實測流量和不同模型(新安江模型、大伙房模型、遼寧指數(shù)模型)的預報流量是否符合正態(tài)分布。P-P圖是一種判斷數(shù)據(jù)序列是否符合正態(tài)分布的重要手段，它是以期望的累積概率作為縱坐標，以觀測的累積概率作為橫坐標，把樣本值以散點的形式落在直角坐標系中。如果樣本序列服從正態(tài)分布，則樣本值應集中落在對角線的附近。如果實測流量和預報流量不服從正態(tài)分布，采用Box-Cox轉(zhuǎn)換法對實測流量和預報流量進行正態(tài)轉(zhuǎn)換。本文采用Minitab工具實現(xiàn)數(shù)據(jù)的Box-Cox轉(zhuǎn)換，采用SPSS工具實現(xiàn)P-P圖的繪制。以“20130816”實測洪水為例，轉(zhuǎn)換前和轉(zhuǎn)換后的實測數(shù)據(jù)正態(tài)頻率曲線見圖2和圖3。

圖2 轉(zhuǎn)化前“20130816”實測流量

表1 歷史洪水信息

圖3 轉(zhuǎn)化后“20130816”實測流量

由圖2、3可以看出，“20130816”實測流量在轉(zhuǎn)換前是高度不服從正態(tài)分布的，經(jīng)過Box-Cox轉(zhuǎn)換后，數(shù)據(jù)明顯接近正態(tài)分布。

4.2.2參數(shù)估計

表2 不同流量級別的模型參數(shù)估計值

由表2可見，在高水情況下，大伙房模型的權(quán)重系數(shù)明顯高于新安江模型和遼寧指數(shù)模型的權(quán)重系數(shù)，分別為0.72 、0.17、0.11，說明基于超滲產(chǎn)流機理的大伙房模型針對高水具有較好的預報效果；在中水情況下，新安江模型、大伙房模型和遼寧指數(shù)模型的權(quán)重系數(shù)基本相同，分別為0.36、0.33 、0.31，說明3個模型針對中水的預報效果基本相同；在低水情況下，新安江模型、遼寧指數(shù)模型的權(quán)重系數(shù)要高于大伙房模型的權(quán)重系數(shù)，分別為0.37、0.21、0.42 ，說明基于蓄滿產(chǎn)流機理的新安江模型、遼寧指數(shù)模型針對低水具有較好的預報效果。其原因為：在中國北方地區(qū)很多地方是超滲產(chǎn)流與蓄滿產(chǎn)流交替產(chǎn)生，當降雨量或雨強較大時，以超滲產(chǎn)流為主；當降雨量或雨強較小時，以蓄滿產(chǎn)流為主。

表3 分類組合預報率定期預報精度分析

注：NXAJ為分類洪水新安江模型、NDHF為分類洪水大伙房模型、NLNZS為分類洪水遼寧省指數(shù)模型、NBMA為分類洪水貝葉斯平均法。

4.2.3評定指標

本文采用相對誤差及確定性系數(shù)兩個指標來分別評定單個模型分類預報結(jié)果和BMA組合預報的預報精度，如上所述；同時，采用覆蓋率作為評定指標，對BMA組合預報的預報區(qū)間的優(yōu)良性進行評定。覆蓋率(CR)是指預報區(qū)間覆蓋實測流量數(shù)據(jù)的比率，為最常用的預報區(qū)間評價指標。計算公式如下

(15)

式中，At為t時刻預報區(qū)間是否覆蓋實測流量，如果覆蓋實測流量則為1，否則為0，NT為預報時間序列長度。CR值越大，表示預報區(qū)間覆蓋率越高。

4.2.4預報值對比分析

采用上述估計的參數(shù)對率定和檢驗的歷史洪水，分別采用單個模型(新安江模型、大伙房模型和遼寧指數(shù)模型)和分類貝葉斯平均法進行洪水模擬，獲取每場洪水的預報結(jié)果(洪峰、峰現(xiàn)時間和預報區(qū)間)，并采用洪峰相對誤差、峰現(xiàn)時間誤差和確定性系數(shù)對模擬結(jié)果進行評定，率定洪水及檢驗洪水特征值統(tǒng)計結(jié)果絕對值見表3和表4。

按照洪水類型(高水、中水和低水)對率定期和檢驗期預報結(jié)果絕對值進行統(tǒng)計，統(tǒng)計結(jié)果如表5所示。

由表3、4可以看出，在率定期和檢驗期的25場洪水中，分類組合預報的洪峰相對誤差比單個模型小，洪峰相對誤差小于等于20%的洪水次數(shù)為22，合格率達到88%；峰現(xiàn)時間誤差小于等于1個時段的洪水次數(shù)為21，合格率為84%。因此，分類組合預報的洪峰預報標準可達到甲級標準，顯著提高了不同類型洪水的預報精度。

從表5統(tǒng)計結(jié)果可看出，除極少數(shù)的洪水以外，分類洪水貝葉斯平均法預報效果優(yōu)于單個模型預報效果。在率定期和檢驗期的25場歷史典型洪水中，總體平均洪峰相對誤差絕對值從12.79、10.9、12.77降到7.05，總體平均峰現(xiàn)時間誤差從1.21、1.03、1.17 h降到0.64 h，總體平均確定性系數(shù)從0.85、0.86、0.83提高到0.90。這表明分類洪水貝葉斯平均法可以應用于大伙房水庫上游流域，效果整體上優(yōu)于單個模型情況，方法可行、有效。

表4 分類組合預報檢驗期預報精度分析

表5 分類組合預報不同洪水類型統(tǒng)計結(jié)果

為了更形象地說明BMA在水文模型綜合預報中的效果，選擇“20130816”洪水作為典型洪水，分別繪制降雨、實測流量，以及NXAJ模型、NDHF模型、NLNZS模型和NBMA模型的預報流量過程線圖，洪水模擬過程如圖4所示。

圖4 20130816洪水單模型預報及NBMA組合預報過程

從圖4可見：分類組合預報模型比單個模型能夠更好地模擬每場洪水的洪峰流量、峰現(xiàn)時間及洪水過程。分類組合預報的洪水過程與實測水位擬合度明顯優(yōu)于單個水文模型，分類組合預報模型集合了單個模型的優(yōu)勢，其預報結(jié)果更貼合實測值?？傮w上，分類組合預報可以達到令人滿意的預報精度，無論是洪峰流量、峰現(xiàn)時間還是洪水過程的模擬，其預報效果都要優(yōu)于單個模型。這表明分類組合預報可有效提高預報精度，該方法是可行的。

4.2.5預報區(qū)間對比分析

除均值預報外，分類貝葉斯平均法還可以給出預報區(qū)間。采用上述估計的參數(shù)對率定期和檢驗期的25場歷史典型洪水，采用分類貝葉斯平均法進行洪水模擬，獲取每場洪水的預報區(qū)間，并采用覆蓋率作為指標來分析比較NBMA組合預報的90%預報區(qū)間與實測值之間的對比，從而分析比較其預報不確定性區(qū)間的優(yōu)良性。分類貝葉斯平均法在整個流量序列90%置信度的置信區(qū)間率定期和檢驗期預報結(jié)果分別見表6、7。

按照洪水類型對率定期和檢驗期分析結(jié)果進行統(tǒng)計，統(tǒng)計結(jié)果如表8所示。

由表6～8可知，分類貝葉斯平均法計算的預報區(qū)間對實測值的覆蓋率在高水部分時效果較好，實測洪峰流量均在置信區(qū)間范圍內(nèi)，覆蓋率均在80%以上，平均覆蓋率為88.65%;但是中水和低水部分的平均覆蓋率較低，分別為52.38%和33.46%。

為了更形象地說明NBMA在水文模型綜合預報中的效果，以“20130816”洪水為例，分別采用NBMA方法和單個模型進行計算，繪制出降雨過程、實測流量、預報均值、90%置信區(qū)間的過程圖，從而分析預報區(qū)間的優(yōu)良性(見圖5)。

從圖5預報區(qū)間分布圖可以看出，“20130816”洪水的絕大多數(shù)的實測數(shù)據(jù)均落于90%置信區(qū)間內(nèi)，區(qū)間對實測值的覆蓋度為88.32%。基于NBMA的預報方法定量的提供了洪水流量范圍，使決策人員定量地考慮預報的不確定性，為防洪規(guī)劃提供了依據(jù)。

表6 分類組合預報模型率定期預報區(qū)間分析 m3/s

表7 分類組合預報模型檢驗期預報區(qū)間分析 m3/s

表8 不同洪水類型預報區(qū)間統(tǒng)計結(jié)果

圖5 NBMA模擬序列90%不確定性置信區(qū)間

4.3 結(jié)果分析

綜上所述，貝葉斯模型加權(quán)平均方法是通過綜合幾個模型預報值的后驗分布來推斷預報量的更可靠概率分布分析工具。它不僅能提供一個綜合的預報值，還能提供一個綜合的預報區(qū)間。本文采用3個水文模型，統(tǒng)一用遺傳算法率定參數(shù)，得到9組不同的預報值用于BMA方法的綜合；然后對BMA和組成其3個模型從預報精度和預報不確定區(qū)間這兩方面進行了詳細的比較和分析，得到如下結(jié)論：

(1)相對誤差和確定性系數(shù)兩個指標的值表明，BMA方法在一定程度上提高了預報精度，對于整個流量序列，BMA的預報值比單個模型的預報值都高；對于不同洪水類型，BMA的預報精度在高水和低水部分都高于單個模型。

(2)利用BMA90%不確定性區(qū)間對模型結(jié)構(gòu)不確定性進行評價。

綜合而言，對BMA和組成其單個模型比較分析表明，BMA方法降低了單一水文預報結(jié)果的不確定性，能夠在一定程度上提高預報精度。更重要的是，BMA能夠提供可靠的預報不確定性區(qū)間，為防洪決策提供參考依據(jù)。

5 結(jié) 論

本文針對水文模型結(jié)構(gòu)的不確定性，基于貝葉斯平均法開展多模型組合預報研究。綜合考慮模型參數(shù)及模型結(jié)構(gòu)不確定性共同對預報結(jié)果造成的影響，在洪水分類預報研究成果基礎上，創(chuàng)新性構(gòu)建了基于BP神經(jīng)網(wǎng)絡和貝葉斯平均法的洪水分類組合預報模型。以大伙房水庫為實例驗證區(qū)域，選取1960年～2016年的25場歷史典型洪水作為本研究的基礎資料，實現(xiàn)對模型的率定與檢驗。結(jié)論如下：

(1)通過貝葉斯平均法推求新安江模型、大伙房模型、遼寧指數(shù)模型3個模型對應不同類型洪水的預報流量序列的均值、方差，利用后驗概率作為每個模型所占權(quán)重，反映了每個模型的預報效果。分析表明，在高水情況下，大伙房模型的權(quán)重系數(shù)明顯高于新安江模型和遼寧指數(shù)模型的權(quán)重系數(shù)，分別為0.72、0.17、0.11，說明基于超滲產(chǎn)流機理的大伙房模型針對高水具有較好的預報效果；在中水情況下，新安江模型、大伙房模型和遼寧指數(shù)模型的權(quán)重系數(shù)基本相同，分別為0.36、0.33、0.31，說明三個模型針對中水的預報效果基本相同；在低水情況下，新安江模型、遼寧指數(shù)模型的權(quán)重系數(shù)要高于大伙房模型的權(quán)重系數(shù)，分別為0.37、0.42、0.21，說明基于蓄滿產(chǎn)流機理的新安江模型、遼寧指數(shù)模型針對低水具有較好的預報效果。分析原因：在中國北方地區(qū)很多地方是超滲產(chǎn)流與蓄滿產(chǎn)流交替產(chǎn)生，當降雨量或雨強較大時，以超滲產(chǎn)流為主，當降雨量或雨強較小時，以蓄滿產(chǎn)流為主。

(2)按照洪水類型對率定期和檢驗期預報結(jié)果進行統(tǒng)計的結(jié)果表明，分類組合預報效果明顯優(yōu)于單個模型預報效果，平均洪峰相對誤差絕對值從10.43、9.43、10.15降到6.36，平均峰現(xiàn)誤差從1.19 h、1.02 h、1.6 h降到0.54 h，平均確定性系數(shù)從0.85、0.86、0.83提高到0.90。選擇“20130816”洪水作為典型洪水進行過程分析，檢驗表明組合預報結(jié)果擬合度明顯優(yōu)于單個水文模型，分類組合預報模型集合了單個模型的優(yōu)勢，預報結(jié)果更貼合實測值。貝葉斯平均法不僅能提高預報精度，還能推求出性質(zhì)更為優(yōu)良的預報區(qū)間，對預報不確定性的定量評估，提高預報的可靠性。分析表明，實測值的覆蓋率在高水部分時效果較好，實測洪峰流量均在置信區(qū)間范圍內(nèi)，覆蓋率均在80%以上，平均覆蓋率為88.65%，中水和低水部分的平均覆蓋率較低，分別為52.38%和33.46%。