周 夢(mèng)，陳 華，郭富強(qiáng)，許崇育

(武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

0 引 言

洪水預(yù)報(bào)是根據(jù)前期和現(xiàn)時(shí)的水文、氣象等信息，揭示和預(yù)測(cè)洪水的發(fā)生及其變化過程的應(yīng)用科學(xué)技術(shù)。進(jìn)行洪水預(yù)報(bào)的水文模型主要是采用觀測(cè)到的歷史水文資料和數(shù)據(jù)，確定模型參數(shù)的合理值，然后對(duì)未來(lái)的水文數(shù)據(jù)進(jìn)行預(yù)測(cè)從而進(jìn)行洪水預(yù)報(bào)。當(dāng)系統(tǒng)的參數(shù)隨時(shí)間變化時(shí)，普通的參數(shù)估計(jì)將無(wú)法得到確定的最優(yōu)值，實(shí)時(shí)校正就是追蹤系統(tǒng)參數(shù)的變動(dòng)過程尋找合適的參數(shù)來(lái)進(jìn)行系統(tǒng)模擬[1]。

在實(shí)時(shí)洪水預(yù)報(bào)系統(tǒng)中，每次預(yù)報(bào)做出之前，實(shí)時(shí)校正模型根據(jù)當(dāng)時(shí)的實(shí)測(cè)信息(新息)對(duì)預(yù)報(bào)模型的參數(shù)、狀態(tài)變量、輸入向量或預(yù)報(bào)值進(jìn)行某種校正，使其更符合客觀實(shí)際?，F(xiàn)在模型洪水預(yù)報(bào)實(shí)時(shí)校正涉及兩方面的研究?jī)?nèi)容：一是建立洪水實(shí)時(shí)預(yù)報(bào)校正模型；二是采用有效的實(shí)時(shí)校正技術(shù)[2]。實(shí)時(shí)洪水預(yù)報(bào)校正技術(shù)有典型的水文模型流量預(yù)報(bào)實(shí)時(shí)校正算法[1]，其在多站點(diǎn)的大流域和流量陡漲陡落的流域局限性比較大；誤差自回歸校正算法[3]關(guān)鍵在于其回歸系數(shù)的確定，確定之后無(wú)需改動(dòng)模型，計(jì)算量較少，但是其使用效果跟預(yù)報(bào)的殘差與前期殘差序列的一致性;衰減記憶最小二乘算法[4]其重點(diǎn)在于實(shí)時(shí)洪水遺忘因子的確定，但是洪水過程隨機(jī)性較大，可以針對(duì)不同洪水選取合適的遺忘因子；卡爾曼濾波算法[5]可以有效濾掉實(shí)測(cè)及預(yù)報(bào)噪聲，關(guān)鍵在于模型噪聲和量測(cè)噪聲方差矩陣的確定；基于貝葉斯理論的自適應(yīng)實(shí)時(shí)校正模型[6]其能夠充分利用先驗(yàn)信息以及實(shí)測(cè)信息得到預(yù)報(bào)的后驗(yàn)分布，關(guān)鍵在于誤差分布的估計(jì)以及模型參數(shù)的確定。這些方法都能根據(jù)實(shí)測(cè)信息及時(shí)更新預(yù)報(bào)數(shù)據(jù)，來(lái)提高預(yù)報(bào)精度。

現(xiàn)有研究中對(duì)于某一流域不同方法的適用性比較較少，并且少有比較各種方法對(duì)于不同洪水評(píng)價(jià)指標(biāo)的表現(xiàn)的文章，所以本文在建溪流域進(jìn)行實(shí)例研究，在新安江水文模型模擬的基礎(chǔ)上比較了誤差自回歸實(shí)時(shí)校正模型、卡爾曼濾波算法以及基于貝葉斯理論的實(shí)時(shí)校正方法的適用性。

1 研究區(qū)域與數(shù)據(jù)

建溪流域位于福建省的北部，流域面積14 787 km2，流域平均降水量1 600～2 600 mm，是福建省鋒面雨高值區(qū)。如圖1所示，建溪七里街水文站以上有崇陽(yáng)溪、南浦溪和松溪三大支流。本文主要研究七里街及其上游的流域，涉及的水文站點(diǎn)有麻沙站，崇安站，管潭站，松溪站，建陽(yáng)站，水吉站，新廠站以及七里街站。

本研究用到的數(shù)據(jù)為2005-2014年的小時(shí)觀測(cè)數(shù)據(jù)，取2005-2013年為率定期，2014年作為檢驗(yàn)期，時(shí)間步長(zhǎng)取1 h，對(duì)場(chǎng)次洪水進(jìn)行計(jì)算。校正主要是對(duì)檢驗(yàn)期新安江模型的模擬數(shù)據(jù)進(jìn)行校正，選取的場(chǎng)次洪水共有10場(chǎng)。

圖1 建溪流域圖Fig.1 Jianxi river basin

2 模型及方法介紹

2.1 新安江模型

新安江模型是由趙人俊教授領(lǐng)導(dǎo)的研究組在20世紀(jì)70年代提出的模擬徑流模型，在80年代由二水源模型發(fā)展成為了三水源模型，適用于濕潤(rùn)地區(qū)與半濕潤(rùn)地區(qū)。三水源新安江模型蒸散發(fā)計(jì)算采用三層蒸散發(fā)模型；產(chǎn)流計(jì)算采用蓄滿產(chǎn)流模型；用自由水蓄水庫(kù)結(jié)構(gòu)將總徑流劃分為地表徑流、壤中流和地下徑流三種；流域匯流計(jì)算采用線性水庫(kù)；河道匯流采用馬斯京根分段連續(xù)演算或滯后演算法[7]。本次率定模型所用的優(yōu)化算法是SCE-UA算法(Shuffled Complex Evolution),SCE-UA算法是一種全局優(yōu)化算法，適合新安江模型這類高維參數(shù)的全局優(yōu)化[8-10]。

新安江模型主要有15個(gè)參數(shù)需要率定[11],分段馬斯京根法對(duì)于不同的河段，參數(shù)k和x需要各自率定。現(xiàn)根據(jù)歷史蒸發(fā)、降雨以及徑流資料，可以用全局優(yōu)化方法率定參數(shù)，得到最優(yōu)的模擬模型。

2.2 AR算法

洪水預(yù)報(bào)誤差序列具有拖尾的特性，且為白噪聲，故洪水預(yù)報(bào)誤差序列為平穩(wěn)時(shí)間序列，可采用時(shí)間序列AR模型進(jìn)行實(shí)時(shí)校正[12]。其主要是對(duì)模型殘差序列采用殘差自回歸估計(jì)式：

et+L=c1et+c2et-1+…+cpet-p+1+ξt+L

(1)

那么預(yù)報(bào)結(jié)果的校正式為：

(2)

其中：

et=Q(t)-QC(t)

2.3 Bayesian算法

(3)

假定天然徑流過程滿足自回歸滑動(dòng)模型，先驗(yàn)分布用自回歸滑動(dòng)平均模型進(jìn)行模擬；假定實(shí)測(cè)徑流和預(yù)報(bào)徑流服從線性關(guān)系，似然函數(shù)使用AR模型進(jìn)行模擬。

2.4 卡爾曼濾波算法

卡爾曼濾波理論的基本思想即采用信號(hào)與噪聲的狀態(tài)空間模型，利用前一時(shí)刻的估計(jì)和觀測(cè)值來(lái)更新對(duì)狀態(tài)變量的估計(jì)，求現(xiàn)時(shí)刻的估計(jì)值[13]。其采用一個(gè)狀態(tài)方程和一個(gè)量測(cè)方程來(lái)完整的描述線性動(dòng)態(tài)過程，前者描述系統(tǒng)狀態(tài)向量隨時(shí)間的變化規(guī)律，后者采用函數(shù)關(guān)系反映量測(cè)向量與系統(tǒng)狀態(tài)向量之間的依賴關(guān)系[14]。其遞推公式既可以得到濾波估計(jì)值，又可以得到誤差的方差陣，即可以完成自身的誤差分析[15]。其公式可以表示為：

狀態(tài)方程：

Xt+1=FtXt+BtUt+wt

(4)

量測(cè)方程：

Yt=HtXt+vt

(5)

式中：Xt為狀態(tài)向量；Yt為測(cè)量向量；Ut為系統(tǒng)的輸入；wt、vt分別為狀態(tài)噪聲和測(cè)量噪聲，且wt、vt為獨(dú)立一致的正態(tài)隨機(jī)變量；Bt、Ht為傳輸矩陣；Ft為系統(tǒng)從t到t+1時(shí)刻的n×n維狀態(tài)轉(zhuǎn)移矩陣[6]。

2.5 建溪流域洪水預(yù)報(bào)模型

本研究主要用到新安江模型及3種實(shí)時(shí)校正模型，如圖2所示，先根據(jù)歷史資料建立建溪流域新安江模型，確定匯流參數(shù)，以及區(qū)間的新安江模型參數(shù)。再用新安江模型進(jìn)行檢驗(yàn)期的流量模擬以及預(yù)報(bào)，預(yù)報(bào)時(shí)假設(shè)預(yù)報(bào)時(shí)刻以后流域降雨為零。預(yù)報(bào)值會(huì)分別經(jīng)過3種校正方法進(jìn)行校正，校正前后的值會(huì)通過分段馬斯京根法匯流到七里街，區(qū)間用新安江模型進(jìn)行預(yù)報(bào)，最后可以匯總得到七里街的流量預(yù)報(bào)。再次用相應(yīng)的方法進(jìn)行校正。

圖2 建溪流域洪水預(yù)報(bào)模型結(jié)構(gòu)圖Fig.2 Technology roadmap

3 結(jié)果與討論

3.1 模型模擬結(jié)果

根據(jù)子流域區(qū)間降雨蒸發(fā)數(shù)據(jù)以及上游流量匯流，新安江模型可以對(duì)子流域流量數(shù)據(jù)進(jìn)行模擬，然后采用分段馬斯京根法進(jìn)行七里街匯流計(jì)算，區(qū)間繼續(xù)采用新安江模型進(jìn)行模擬，得到的預(yù)報(bào)結(jié)果如表1所示。

表1 率定期及檢驗(yàn)期新安江模型模擬結(jié)果Tab.1 Simulated results for calibration and validation period

表1列出了新安江模型模擬的結(jié)果，率定期為2005-2013年，檢驗(yàn)期為14年?？梢钥吹?，率定期的洪水平均納西效率系數(shù)都達(dá)到了0.9以上，子流域的模擬情況普遍比整個(gè)大流域的模擬情況要好。檢驗(yàn)期最小的納西效率系數(shù)都不到0.9，子流域的平均納西效率系數(shù)還是保持在0.9以上。對(duì)于各場(chǎng)洪水模擬的結(jié)果來(lái)看，80%的子流域洪水模擬納西效率系數(shù)都高于0.9，對(duì)于個(gè)別洪水，有的子流域降雨很少，洪峰不明顯，所以模擬效果不夠好，但是對(duì)整個(gè)流域的誤差影響較小。

3.2 各校正方法校正結(jié)果比較分析

子流域的預(yù)測(cè)及校正完成后，用校正后的預(yù)測(cè)數(shù)據(jù)，采用分段馬斯京根法進(jìn)行七里街匯流計(jì)算，區(qū)間繼續(xù)采用新安江模型進(jìn)行模擬預(yù)報(bào)，得到的預(yù)報(bào)結(jié)果如表2所示。

由表2可以看到，3種校正方法都相較于新安江模型模擬值明顯提高了精度，納西效率系數(shù)能增加0.15左右，水量平衡誤差也大大降低。隨著預(yù)見期的增加，預(yù)報(bào)精度降低，就研究洪水的平均水平來(lái)說，對(duì)于前幾個(gè)小時(shí)的預(yù)見期，貝葉斯方法的納西效率系數(shù)和水量平衡誤差這兩個(gè)指數(shù)都優(yōu)于其余兩種方法。在預(yù)見期增加后，AR模型校正的水量相對(duì)誤差比較小，在預(yù)見期大于9之后，卡爾曼濾波的納西效率系數(shù)較高，而貝葉斯方法的相對(duì)誤差較小。綜合來(lái)看，貝葉斯方法在納西效率系數(shù)和水量誤差方面在各場(chǎng)洪水的總體水平上上表現(xiàn)較好。

為了具體分析各種方法不同預(yù)見期的表現(xiàn)，現(xiàn)在對(duì)于每場(chǎng)洪水的納西效率系數(shù)進(jìn)行統(tǒng)計(jì)，找出每場(chǎng)洪水最高和最低的納西效率系數(shù)如表3所示。

各種方法在不同預(yù)見期時(shí)納西效率系數(shù)最高和最低的頻率如圖3所示，由柱狀圖的分布情況可以明顯發(fā)現(xiàn)，貝葉斯方法對(duì)短預(yù)見期表現(xiàn)最好，模擬出納西效率系數(shù)比另外兩種方法好的頻率很高，但是隨著預(yù)見期的增加，貝葉斯方法的優(yōu)勢(shì)呈降低趨勢(shì)，在預(yù)見期超過6 h時(shí)已呈劣勢(shì)。相反地，卡爾曼濾波校正方法對(duì)短預(yù)見期的校正效果不如貝葉斯方法和AR方法，但是在預(yù)見期擴(kuò)大時(shí)，開始顯現(xiàn)出明顯的優(yōu)勢(shì)。這說明卡爾曼濾波校正技術(shù)比貝葉斯方法更適合長(zhǎng)預(yù)見期的洪水。而AR方法很穩(wěn)定，頻數(shù)分布比較均勻，對(duì)各預(yù)見期表現(xiàn)沒有明顯變化趨勢(shì)。綜合比較，在預(yù)見期較短時(shí)，可以選擇貝葉斯方法，而預(yù)見期變長(zhǎng)之后，可以選擇卡爾曼濾波方法。

表2 不同預(yù)見期新安江模型與各校正方法的平均納西效率系數(shù)以及平均水量誤差Tab.2 Average NSE and RE of data forecasted by Xinanjiang model and three correction methods in different prediction period

表3 各校正方法校正之后的納西效率系數(shù)比較Tab.3 Comparison of NSE after correction of three correction methods

圖3 各種方法校正之后得到最高和最低納西效率系數(shù)的頻率比較Fig.3 Frequency comparison of the highest and lowest NSE after correction of three methods

圖4 各方法校正之后對(duì)應(yīng)不同洪峰流量的最高最低納西效率系數(shù)頻率比較Fig.4 Frequency comparison of maximum and minimum NSE corresponding to different peak flow after different methods correction[16] Maximum NSE Minimum NSE

除了考慮不同校正方法跟預(yù)見期的關(guān)系，還考慮到不同方法的校正效果可能跟洪峰流量有一定的關(guān)系。子流域以及整個(gè)流域各場(chǎng)洪水的洪峰流量以及6個(gè)預(yù)見期時(shí)各種校正方法校正效果如圖4所示，可見AR方法在6個(gè)預(yù)見期時(shí)表現(xiàn)較好，得到最優(yōu)納西效率系數(shù)的次數(shù)比較穩(wěn)定，貝葉斯方法對(duì)于低流量洪水表現(xiàn)較好，卡爾曼濾波比較適合高流量的洪水。

洪水預(yù)報(bào)中，對(duì)洪水洪峰的預(yù)測(cè)也十分重要，現(xiàn)在對(duì)不同預(yù)見期預(yù)測(cè)的洪水的洪峰誤差進(jìn)行統(tǒng)計(jì)，得出比較圖如圖5所示。

圖5 各校正方法不同各預(yù)見期的洪峰誤差比較Fig.5 Comparison of flood peak errors at different forecast periods with different correction methods

由圖5可以看到，各個(gè)校正方法在前兩個(gè)預(yù)見期預(yù)報(bào)精度都比較高，洪峰誤差都在±3%以內(nèi)，在預(yù)見期延長(zhǎng)時(shí)，誤差都在擴(kuò)大，三種方法校正效果從圖中可以看到，卡爾曼濾波和貝葉斯方法，較于AR模型洪峰誤差較小，并且從分布區(qū)間和中位數(shù)來(lái)看，卡爾曼濾波表現(xiàn)更好，洪峰誤差較小。

圖6以20140725這場(chǎng)洪水為例，展示了不同預(yù)見期3種方法預(yù)報(bào)的結(jié)果，可以看到，隨著預(yù)見期的增加，預(yù)報(bào)值偏離實(shí)測(cè)值越來(lái)越遠(yuǎn)。在3個(gè)預(yù)見期時(shí)，3種方法的校正曲線比較接近實(shí)測(cè)值，但是已經(jīng)不再光滑，在到12個(gè)預(yù)見期時(shí)，已經(jīng)偏離較遠(yuǎn)了，這是因?yàn)?2個(gè)預(yù)見期，在用新安江模型預(yù)報(bào)時(shí)，假設(shè)沒有降雨，12個(gè)預(yù)見期沒有降雨這在雨量方面的誤差過大，這個(gè)時(shí)候的預(yù)報(bào)精度已經(jīng)不能很好的保證了，只能做一個(gè)大致的參考。

4 結(jié) 語(yǔ)

建溪流域雨量豐沛，氣候濕潤(rùn)，新安江模型本身就可以較好地模擬洪水過程。AR模型，卡爾曼濾波以及貝葉斯方法3種校正方法都能明顯提高流量模擬的精度，對(duì)于不同的評(píng)判指標(biāo)，3種校正方法校正效果有所不同。對(duì)于水量相對(duì)誤差，貝葉斯方法校正后效果最好，AR模型與卡爾曼濾波方法校正結(jié)果不相上下。AR模型和卡爾曼濾波方法校正之后的洪水過程納西效率系數(shù)保持很高，卡爾曼濾波方法對(duì)于長(zhǎng)預(yù)見期或者流量較大的洪水表現(xiàn)得比較好，可以得到更好的納西效率系數(shù)，而貝葉斯方法對(duì)于流量較小的洪水表現(xiàn)得更好，AR模型對(duì)于各個(gè)量級(jí)的洪水都能有比較好的表現(xiàn)。綜上，3種方法都能在實(shí)際應(yīng)用中達(dá)到較為理想的效果，根據(jù)不同的預(yù)報(bào)需求，可以選擇相應(yīng)的方法進(jìn)行校正。

圖6 不同預(yù)見期3種方法預(yù)測(cè)值與實(shí)測(cè)比較Fig.6 Comparison between predicted and measured values of three methods in different forecast periods