基于時(shí)變?cè)鲆嫠哪Ｐ偷膶?shí)時(shí)預(yù)報(bào)研究

劉慧媛，夏 軍,，鄒 磊，洪 思

(1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072；2. 中國科學(xué)院地理科學(xué)與資源研究所 陸地水循環(huán)及地表過程重點(diǎn)實(shí)驗(yàn)室，北京 100101)

近年來，強(qiáng)降雨洪水頻發(fā)，特別是暴雨洪水誘發(fā)的自然災(zāi)害經(jīng)常發(fā)生，洪澇災(zāi)害已然成為影響社會(huì)經(jīng)濟(jì)健康發(fā)展的一大隱患?，F(xiàn)階段，我國大江大河因?yàn)樾藿舜罅抗こ淘O(shè)施以及采用了先進(jìn)的非工程措施取得了很好的防洪效果，但在中小流域，由于其降雨徑流非線性特性顯著，洪水預(yù)報(bào)問題依然面臨著十分嚴(yán)峻的局面。

目前，洪水預(yù)報(bào)從簡單的概念性水文模型、數(shù)學(xué)物理模型、非線性系統(tǒng)理論模型及與物理過程聯(lián)系的分布式流域水循環(huán)模型等[1]，發(fā)展到后來的與GIS 和 RS 等技術(shù)結(jié)合的復(fù)雜數(shù)字流域水文模型，已取得了快速的發(fā)展和應(yīng)用。國內(nèi)外很多學(xué)者針對(duì)洪水預(yù)報(bào)模型已作了大量的研究工作，如斯坦福流域水文模型[2](SWM)、薩克拉門托模型[3](SAC)、水箱模型[4](TANK)、新安江模型[5]、SHE 模型[6]、SWAT模型[7]、TOPMODEL模型[8]、VIC模型[9]等。我國水文學(xué)者夏軍依據(jù) Volterra 泛函級(jí)數(shù)表達(dá)的降雨徑流非線性系統(tǒng)理論，利用從復(fù)雜水文關(guān)系找出簡單關(guān)系的非線性系統(tǒng)分析，提出了水文非線性系統(tǒng)理論的時(shí)變?cè)鲆婺Ｐ?(Time Variant Gain Model，簡記為TVGM) 。論述了造成降雨徑流系統(tǒng)關(guān)系非線性的主要原因是產(chǎn)流過程中由于土壤濕度(即土壤含水量)的變化會(huì)導(dǎo)致產(chǎn)流量發(fā)生變化[10]。該模型參數(shù)少、結(jié)構(gòu)簡單，能以一種簡單的系統(tǒng)關(guān)系模擬復(fù)雜非線性水文過程，模擬效果良好。

然而絕大多數(shù)中小流域地形復(fù)雜，下墊面條件空間分布不均，一般模型難以精準(zhǔn)概述研究區(qū)域的整體情況。除此之外，由于人類活動(dòng)的影響，很多地區(qū)其自然環(huán)境發(fā)生了一定程度的變化，使得長時(shí)期內(nèi)流域產(chǎn)匯流規(guī)律變化較大[11]。而現(xiàn)如今，廣泛使用的水文預(yù)報(bào)模型是由歷史實(shí)測降雨和徑流資料率定得到的，以大多數(shù)洪水能取得好的預(yù)報(bào)結(jié)果為優(yōu)選準(zhǔn)則，尋求最佳模型參數(shù)，其反映的是研究區(qū)域長期下墊面條件的平均情況，并不能準(zhǔn)確描述每一場洪水的徑流過程[12]。除此以外，歷史資料的缺失或觀測誤差及模型結(jié)構(gòu)的不確定性等問題，也會(huì)使得傳統(tǒng)模型在預(yù)報(bào)時(shí)，常常無法避免地出現(xiàn)誤差[13]。因而要想直接從模型下手，構(gòu)建結(jié)構(gòu)簡單、參數(shù)較少、精度夠好的預(yù)報(bào)方案具有較大的難度。所以為了提高預(yù)報(bào)精度，有必要在構(gòu)建水文模型之外采取一定的實(shí)時(shí)校正技術(shù)來修正預(yù)報(bào)誤差。實(shí)時(shí)校正就是實(shí)時(shí)地計(jì)算水文預(yù)報(bào)模型最新出現(xiàn)的預(yù)報(bào)誤差，并以此為依據(jù)，對(duì)模型參數(shù)、狀態(tài)、預(yù)報(bào)輸出值進(jìn)行修正，及時(shí)地對(duì)系統(tǒng)本身的動(dòng)態(tài)變化進(jìn)行跟蹤并利用新的觀測信息，從而使預(yù)報(bào)系統(tǒng)迅速適應(yīng)現(xiàn)實(shí)的狀況，以此來提高水文預(yù)報(bào)的精度。目前常用的實(shí)時(shí)校正方法可以分為兩類：第一類是對(duì)模型內(nèi)部參數(shù)進(jìn)行跟蹤校正，此方法一般精度較高，但同時(shí)需要明確流域產(chǎn)匯流規(guī)律，對(duì)實(shí)測資料要求高，故而這類實(shí)時(shí)校正方法難以在實(shí)測資料匱乏的地區(qū)廣泛適用。第二類是利用預(yù)報(bào)誤差序列，直接對(duì)預(yù)報(bào)輸出值進(jìn)行校正，此方法在有效消減誤差的同時(shí)盡可能簡化誤差處理部分，對(duì)實(shí)測資料要求低，水文模型與實(shí)時(shí)校正模型之間具有一定的獨(dú)立性，可與多種預(yù)報(bào)模型聯(lián)合使用，可應(yīng)用范圍較廣。

本文基于時(shí)變?cè)鲆嫠哪Ｐ?，采用具有可變遺忘因子的遞推最小二乘法對(duì)模型計(jì)算輸出值進(jìn)行實(shí)時(shí)校正，目的在于提升模型預(yù)報(bào)效果，提高中小流域洪水預(yù)報(bào)的精度。時(shí)變?cè)鲆婺Ｐ捅旧砭哂袇?shù)少、結(jié)構(gòu)簡單的優(yōu)點(diǎn)，不僅能較好地模擬產(chǎn)匯流過程，且對(duì)無資料地區(qū)也具有很強(qiáng)的適用性。在實(shí)時(shí)校正方面，目前有很多研究采用的是第一類校正方法，即對(duì)模型參數(shù)進(jìn)行校正，需明確水文模型的物理機(jī)理，可供研究的模型有限，如新安江模型[13]、水箱模型[14]等，實(shí)時(shí)校正過程計(jì)算起來也相對(duì)復(fù)雜。對(duì)時(shí)變?cè)鲆婺Ｐ陀?jì)算結(jié)果進(jìn)行實(shí)時(shí)校正，是在結(jié)構(gòu)簡單的模型基礎(chǔ)上，采用相對(duì)簡化的實(shí)時(shí)修正技術(shù)，計(jì)算量大大減少，會(huì)一定程度地提高模型的時(shí)效性。除此以外，研究表明該方法能夠取得突出的模擬效果，因而不失為一種運(yùn)行高效且計(jì)算精度良好的預(yù)報(bào)方案，對(duì)于防洪減災(zāi)具有重要意義。

1 水文非線性時(shí)變?cè)鲆婺Ｐ?/h2>

時(shí)變?cè)鲆婺Ｐ褪窍能娊淌谟?989-1999年期間，在參加愛爾蘭國立大學(xué)的國際河川水文預(yù)報(bào)研討班時(shí)提出的一種非線性系統(tǒng)模型。傳統(tǒng)的水文線性理論假定系統(tǒng)的增益(Gain factor)為一固定常數(shù)。然而夏軍教授，通過分析全球60多個(gè)不同氣候區(qū)域的流域?qū)崪y長序列水文資料后，發(fā)現(xiàn)其增益并非常數(shù)，它與土壤濕度、流域下墊面特性以及氣候特性等影響因素息息相關(guān)，并表現(xiàn)出地球水文系統(tǒng)固有的指數(shù)關(guān)系[10]。本文采用改進(jìn)后的多水源水文時(shí)變?cè)鲆婺Ｐ蚚15]，將產(chǎn)匯流分別分為地表和地下兩部分，以下為具體計(jì)算公式。

1.1 產(chǎn) 流

(1)在地表產(chǎn)流中，可利用系統(tǒng)增益G(t)與毛雨X(t)的乘積來表示有效凈雨Rs的值。其中，通過單一線性水庫系統(tǒng)來計(jì)算流域土壤前期影響雨量API的值，進(jìn)而利用產(chǎn)流參數(shù)與API來計(jì)算系統(tǒng)增益G(t)。

Rs(t)=G(t)X(t)=[g1+g2API(t)]X(t)

(1)

式中：g1和g2均為地表產(chǎn)流參數(shù)。

(2)在計(jì)算地下產(chǎn)流量時(shí)，有效凈雨Rg可表示為產(chǎn)流參數(shù)與API的乘積。

Rg(t)=g3API(t)

(2)

式中：g3為地下產(chǎn)流參數(shù)。

1.2 匯 流

(1)地表匯流采用Nash單位線進(jìn)行計(jì)算，但獲取到的凈雨資料一般為離散形式，因而通常是將瞬時(shí)單位線u(0,t)轉(zhuǎn)換為時(shí)段單位線q(Δt,t)，再由式(3)即可計(jì)算得到地面匯流量Qs。

(3)

(2)地下水匯流量Qg本文選用線性水庫演算法。

Qg2=Rg(1-KKG)U+KKGQg1

(4)

(5)

式中：KKG為產(chǎn)流參數(shù)；U為折算系數(shù)；F為流域面積，km2；Δt為時(shí)間步長，本文為1 h。

2 模型參數(shù)識(shí)別原理和方法

系統(tǒng)識(shí)別指的是通過觀測一個(gè)系統(tǒng)或過程的輸入-輸出關(guān)系來確定該系統(tǒng)過程的數(shù)學(xué)模型[10]，也叫參數(shù)率定。本文采用SCE-UA算法[16]進(jìn)行模型參數(shù)率定。該算法綜合了隨機(jī)搜索、確定性搜索和生物進(jìn)化等方法的優(yōu)點(diǎn)，在算法中引入種群概念[17]。使復(fù)合形點(diǎn)在可行域內(nèi)被隨機(jī)生成和演化，通過對(duì)各個(gè)復(fù)合體的定期洗牌和重組來確保每個(gè)復(fù)合體獲得的信息在整個(gè)空間范圍內(nèi)得到共享，以產(chǎn)生最優(yōu)點(diǎn)從而達(dá)到率定參數(shù)的目的。大量實(shí)驗(yàn)研究成果表明，該算法優(yōu)化效果良好，能夠快速、一致、穩(wěn)定地收斂到全局最優(yōu)解。

3 實(shí)時(shí)校正及方法

科學(xué)技術(shù)的進(jìn)步大大提高了水文預(yù)報(bào)的預(yù)測精度，延長了預(yù)見期。但對(duì)于許多洪水事件，特別是暴雨洪水，由于其復(fù)雜多變的特征，現(xiàn)有的洪水預(yù)報(bào)技術(shù)還不完善，還不能滿足社會(huì)對(duì)洪水預(yù)報(bào)的及時(shí)性和可靠性的要求。絕大多數(shù)水文模型都是利用歷史資料確定好模型參數(shù)，然后用于未來洪水預(yù)報(bào)，即研究的是時(shí)不變的離線系統(tǒng)。這種方式在洪水預(yù)報(bào)預(yù)警中常常得不到令人滿意的結(jié)果。洪水預(yù)報(bào)誤差主要來源于以下幾個(gè)方面：①初始狀態(tài)誤差：模型的初始值由于前期資料的缺乏具有不確定性，因而會(huì)產(chǎn)生一定誤差；②模型結(jié)構(gòu)誤差：任何一種模型都不能完整描述復(fù)雜多變的降雨徑流關(guān)系，目前還沒有天衣無縫的模型結(jié)構(gòu)；③參數(shù)識(shí)別的誤差：率定得到的參數(shù)與其真值之間總會(huì)有一定的誤差，任何一種率定方法都不能得到完美的率定值。且參數(shù)會(huì)因外部條件的改變而發(fā)生變化，多場洪水資料率定的結(jié)果并不能保證會(huì)適用于未來每一場洪水；④觀測誤差：受客觀條件影響，觀測值會(huì)有一定的誤差。

為了對(duì)系統(tǒng)本身的動(dòng)態(tài)變化進(jìn)行跟蹤和利用新的觀測信息，就需要對(duì)水文預(yù)報(bào)模型采用“在線”識(shí)別的方式進(jìn)行實(shí)時(shí)校正。實(shí)時(shí)校正的核心技術(shù)是利用“新信息”對(duì)系統(tǒng)模型進(jìn)行實(shí)時(shí)地修正。目前在洪水預(yù)報(bào)中，常用的實(shí)時(shí)校正算法有：誤差自回歸校正算法；卡爾曼濾波法；遞推最小二乘法；等。誤差自回歸算法[18]在實(shí)際模擬過程中首先需要確定模型的階數(shù)，并選用合適的方法進(jìn)行系數(shù)估計(jì)。而卡爾曼濾波法[18]需要將模型線性化，校正模型建立不易，且易受狀態(tài)向量、觀測向量選擇的影響，受噪聲協(xié)方差矩陣影響較大，可能出現(xiàn)校正效果不穩(wěn)定的問題。故本文采用具有可變遺忘因子的遞推最小二乘法對(duì)模型計(jì)算結(jié)果進(jìn)行實(shí)時(shí)校正。

最小二乘法的實(shí)質(zhì)就是讓實(shí)測數(shù)據(jù)和預(yù)報(bào)數(shù)據(jù)之間的距離平方和最小。而遞推最小二乘法就是在它的基礎(chǔ)上改進(jìn)的，省去了求逆的步驟，其核心是根據(jù)現(xiàn)時(shí)輸入的信息來實(shí)時(shí)更新預(yù)報(bào)誤差的權(quán)重，以達(dá)到實(shí)時(shí)校正的目的。將可變的遺忘因子λ引入遞推最小二乘法中，能夠避免出現(xiàn)“數(shù)據(jù)飽和”現(xiàn)象，同時(shí)能夠快速跟蹤系統(tǒng)變化[19]。

算法的計(jì)算步驟如下：

(1)給定初始值θ(0)、P(0)、R。

(2)已知{y(t),t=0,1,…,T}為實(shí)時(shí)觀測的流量序列與模型模擬流量序列的殘差序列。

(3)設(shè)t=0。

(4)校正后的殘差序列：

(6)

矩陣X為：

X(t)=[y(t-1),y(t-2),…,y(t-m)]

(7)

式中：n為模型階數(shù)。

(5)進(jìn)行第t+1次采樣得y(t+1)。

(6)計(jì)算增益矩陣：

(8)

(7)計(jì)算遺忘因子：

(9)

當(dāng)λ(t+1)<λmin時(shí)，λ(t+1)=λmin；當(dāng)λ(t+1)>λmax時(shí)，λ(t+1)=λmax；

使得λ控制在一定范圍之內(nèi)，并且有0<λ<1， 愈小，對(duì)系統(tǒng)變化的跟蹤能力就越強(qiáng)，但估計(jì)誤差會(huì)變小。

(8)修正參數(shù)估計(jì)：

(10)

(9)將校正后的殘差序列疊加到模型模擬流量序列上得到校正后的流量序列，若序列全部校正完畢則結(jié)束，否則令t=t+1，轉(zhuǎn)到(4)。

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

本文選取湖北省3個(gè)中小流域(高家堰、西河驛、漁洋關(guān))作為研究區(qū)域，搜集整理各流域多年次洪資料(時(shí)間步長為1 h)，應(yīng)用TVGM模型進(jìn)行次洪模擬，并進(jìn)行實(shí)時(shí)校正。將多場洪水按要求劃分為率定期和檢驗(yàn)期，首先在率定期進(jìn)行參數(shù)率定，得到優(yōu)選的模型參數(shù)；繼而在檢驗(yàn)期對(duì)模型進(jìn)行檢驗(yàn)，并在檢驗(yàn)期實(shí)行實(shí)時(shí)校正，將校正后的結(jié)果與初始洪水預(yù)報(bào)結(jié)果作對(duì)比分析。表1所示為研究區(qū)域的基本信息。

表1 所研究中小流域基本信息

4.1 高家堰流域次洪模擬

該流域設(shè)有高家堰水文站，位于東經(jīng)111°03′，北緯30°36′，地處長陽土家族自治縣高家堰鎮(zhèn)，河名為丹水，流域面積334 km2。測驗(yàn)河段順直長400 m，河床全部由卵石組成，斷面沖淤頻繁。兩岸為高山，植被覆蓋度良好。右岸灘栽種有柳林，左岸為公路，路基由塊石砌坎。本站控制流域多年平均降水量為1 567.2 mm，歷時(shí)實(shí)測年最大降水量達(dá)2 307.4 mm。該流域內(nèi)有4個(gè)配套雨量站，它們分別是堡子、白沙驛、染房坪和高家堰雨量站。圖1展示了流域水系及站點(diǎn)分布。

(1)模型參數(shù)率定值如表2所示。

表2 參數(shù)率定值

(2)檢驗(yàn)期各場次洪水模型模擬結(jié)果如圖2所示。

圖1 高家堰流域水系及站點(diǎn)分布圖

4.2 西河驛流域次洪模擬

西河驛流域位于東經(jīng)115°26′～115°51′，北緯30°14′～30°38′，河名蘄水，流域面積1 971 km2，干流總長120 km。西河驛水文站位于距河口23 km處，該站斷面以上流域集水面積占流域總面積的91%，為1 800 km2。流域地勢為東北高，西南低，上游為山區(qū)，中、下游為丘陵、平原區(qū)。當(dāng)?shù)貙賮啛釒Ъ撅L(fēng)性濕潤氣候，流域多年平均年降雨量為1 384.6 mm，年最小降雨量出現(xiàn)于1968年，為800.6 mm，最大降雨量出現(xiàn)在1954年，為2 292.9 mm，多年平均徑流系數(shù)為0.51。本流域內(nèi)有9個(gè)配套雨量站，它們分別是大同、花園、劉河、獅子口、獅子堰、田橋、望天畈，張家榜、株林河雨量站。圖3展示了西河驛流域水系及站點(diǎn)分布狀況。

(1)模型參數(shù)率定值如表3所示。

表3 參數(shù)率定值

圖2 高家堰流域各場次洪水模擬效果圖

(2)檢驗(yàn)期各場次洪水模型模擬結(jié)果如圖4所示。

4.3 漁洋關(guān)流域次洪模擬

漁洋關(guān)水文站位于湖北省宜昌市五峰土家族自治縣漁洋關(guān)鎮(zhèn)，東經(jīng)111°04′，北緯30°10′，屬清江一級(jí)支流漁洋河的代表站，集水面積465 km2。測驗(yàn)河段較順直，長約300 m，河床由卵石組成，沖淤變化較小常年保持穩(wěn)定，水位流量關(guān)系控制條件良好。本流域內(nèi)有4個(gè)配套雨量站，它們分別是王家沖、長樂坪、高橋、漁洋關(guān)。圖5展示了漁洋關(guān)流域水系及站點(diǎn)分布狀況。

(1)模型參數(shù)率定值如表4所示。

圖3 西河驛站流域水系及站網(wǎng)分布圖

參數(shù)名稱g1g2g3NKKKG參數(shù)值0.35 0 0.58 1.67 20.00 0.80

(2)檢驗(yàn)期各場次洪水模型模擬結(jié)果如圖6所示。

4.4 總結(jié)與分析

表5給出了各流域檢驗(yàn)期實(shí)時(shí)校正前后結(jié)果對(duì)比。

表5 各流域校正前后對(duì)比表

綜合上述結(jié)果來看，時(shí)變?cè)鲆婺Ｐ驮诟呒已摺⑽骱芋A、漁洋關(guān)3個(gè)流域有著較好的模擬效果。檢驗(yàn)期總體納什效率系數(shù)均值為0.78；峰現(xiàn)時(shí)間預(yù)報(bào)合格率在3個(gè)流域分別為乙級(jí)、未達(dá)丙級(jí)、甲級(jí)；洪峰流量預(yù)報(bào)效果相對(duì)較差，在漁洋關(guān)流域達(dá)乙級(jí)標(biāo)準(zhǔn)，在高家堰和西河驛流域均未達(dá)丙級(jí)，尤其在西河驛流域，由圖4(c)、圖4(e)、圖4(h)可以看出，當(dāng)流域面積較大，洪水洪峰流量偏大時(shí)，模型預(yù)報(bào)值與實(shí)測值間的偏差會(huì)較大。采用具有可變遺忘因子的遞推最小二乘法對(duì)模型進(jìn)行實(shí)時(shí)校正后，預(yù)報(bào)效果有了顯著的提高。由圖2(c)、圖4(d)、圖4(e)、圖4(g)、圖4(h)可以看出，校正后的峰現(xiàn)時(shí)間都更加接近實(shí)測值；由圖2(a)、圖2(f)、圖4(c)、圖4(e)、圖4(g)、圖4(h)可以看出，洪峰流量預(yù)報(bào)效果也有了明顯的改善，尤其在圖4(c)這種陡漲陡落的情況下，經(jīng)實(shí)時(shí)校正后的洪峰流量預(yù)報(bào)誤差減小了50%以上，洪水漲落部分也更加貼合實(shí)測值。檢驗(yàn)期3個(gè)流域總體納什效率系數(shù)均值達(dá)0.95，各流域峰現(xiàn)時(shí)間預(yù)報(bào)合格率以及洪峰流量預(yù)報(bào)合格率也都達(dá)到乙級(jí)及以上標(biāo)準(zhǔn)。耦合實(shí)時(shí)預(yù)報(bào)技術(shù)后，模型預(yù)報(bào)精度得到顯著提升，可應(yīng)用于生產(chǎn)部門的作業(yè)預(yù)報(bào)。因此，該預(yù)報(bào)方案對(duì)于流域洪水預(yù)報(bào)具有一定的指導(dǎo)意義。

圖4 西河驛流域各場次洪水模擬效果圖

圖5 漁洋關(guān)流域DEM圖及站點(diǎn)分布

5 結(jié) 語

洪水預(yù)報(bào)是流域防洪調(diào)度決策系統(tǒng)的關(guān)鍵組成部分，是實(shí)現(xiàn)科學(xué)、高效、可靠防洪調(diào)度的基礎(chǔ)。本文基于多水源時(shí)變?cè)鲆婺Ｐ?，采用具有可變遺忘因子的遞推最小二乘法對(duì)模型預(yù)報(bào)結(jié)果實(shí)行實(shí)時(shí)校正。

上述預(yù)報(bào)方案經(jīng)過湖北省3個(gè)流域(高家堰、西河驛、漁洋關(guān))次洪資料的檢驗(yàn)，獲得了滿意的結(jié)果。經(jīng)實(shí)時(shí)校正后，多水源時(shí)變?cè)鲆婺Ｐ皖A(yù)報(bào)結(jié)果更加接近實(shí)測洪水過程，尤其是顯著地改善了當(dāng)洪峰流量較大時(shí)，模型預(yù)報(bào)值偏小的狀況。各流域納什效率系數(shù)達(dá)甲級(jí)標(biāo)準(zhǔn)，峰現(xiàn)時(shí)間預(yù)報(bào)合格率以及洪峰流量預(yù)報(bào)合格率也都達(dá)到乙級(jí)及以上標(biāo)準(zhǔn)。該預(yù)報(bào)方案簡單高效，適用性強(qiáng)，能夠一定程度地提高水文預(yù)報(bào)的精度和時(shí)效性，在實(shí)際應(yīng)用中有良好的推廣前景。

圖6 漁洋關(guān)流域各場次洪水模擬效果圖