產(chǎn)流誤差的動態(tài)系統(tǒng)響應曲線修正方法應用

闕家駿 包為民，2

（1.河海大學 水文水資源學院，南京 210098；2.河海大學 水文水資源與水利工程科學國家重點實驗室，南京 210098）

洪水預報是防洪非工程措施中的一項重要內(nèi)容，準確及時的洪水預報對防洪減災、水庫調(diào)度興利、工農(nóng)業(yè)安全生產(chǎn)等意義重大［1］.在實際流域中，由于存在許多誤差因素，如設備故障，導致資料缺測或不合理的觀測數(shù)據(jù)；水利工程、農(nóng)田蓄放水誤差；流域水文規(guī)律的變化；水文規(guī)律簡化誤差，即模型結(jié)構(gòu)誤差，故需利用實時系統(tǒng)能獲得的觀測信息和一切能利用的其他信息對預報誤差進行實時修正，以彌補流域水文模型的不足［2］.洪水預報實時修正方法有很多種［3－6］，如AR模型、動態(tài)系統(tǒng)響應曲線、卡爾曼濾波、遞推最小二乘估計等.AR模型結(jié)構(gòu)簡單、計算簡便，但前提是流量誤差的前后時段具有較強相關(guān)性，否則修正效果不甚理想.卡爾曼濾波較為復雜，可容納的信息量大，但實時洪水預報常難以提供修正所需的足夠信息，這就大大限制了方法的使用，導致其修正效果與AR模型效果類似.本文根據(jù)東圳流域的實際情況建立概念性水文模型，分別采用DSRC和二階AR模型對預報結(jié)果進行實時修正，并將修正結(jié)果與預報結(jié)果進行比較分析.

1 流域概況

東圳水庫位于木蘭溪支流延壽溪中游，是一座以灌溉為主，結(jié)合防洪、發(fā)電、供水、養(yǎng)殖等綜合利用的多年調(diào)節(jié)的大型水庫，流域面積321km2，流域內(nèi)植被較好，平均年降水量1 740mm，多年平均入庫水量3.02億m3，水庫供應莆田市的生活用水和工農(nóng)業(yè)用水，特別是湄洲灣北岸工業(yè)發(fā)展用水，是莆田市的生命之源［7］.

流域?qū)賮啛釒ШＱ笮约撅L氣候，10月～次年3月為旱季，4～9月為雨季，雨量充沛，累積雨量約占全年雨量的80%.其中4月下旬～6月下旬為梅雨季，常出現(xiàn)持續(xù)連綿的陰雨，洪水特點是歷時長、洪量大、但洪峰一般較小，累積總雨量約占全年降雨量的1／3，稱前汛；7～9月為臺風季節(jié)，臺風夾帶降雨一般歷時短，雨量集中、強度大、洪水特點是洪峰較大，稱后汛.洪水陡漲陡落，匯流較快，洪水過程一般為3d左右.

2 東圳流域概念性水文模型

根據(jù)東圳流域的地理位置、降水、氣候、下墊面條件和洪水過程線分析，屬于典型的濕潤地區(qū)，擬選用三水源新安江模型.模型分為蒸散發(fā)、產(chǎn)流、分水源和匯流計算4個層次結(jié)構(gòu)，其中蒸散發(fā)計算為3層蒸散發(fā)模式；產(chǎn)流計算為蓄滿產(chǎn)流模式；分水源計算為自由水蓄水庫結(jié)構(gòu)；坡面匯流計算為單位線或線性水庫；河道匯流計算為馬斯京根法.模型將流域分成若干個單元流域，對劃分好的每塊單元流域分別進行蒸散發(fā)、產(chǎn)流、水源劃分和匯流計算，得到單元流域出口的流量過程；對單元流域出口的流量過程進行出口以下的河道匯流計算，得到該單元流域在全流域出口的流量過程；將每塊單元流域在全流域出口的流量過程線性疊加，即為全流域出口總的流量過程，模型結(jié)構(gòu)如圖1所示［2］.圖中下標i為單元站點；P為降水量；IMP為各單元流域不透水面積比；Ew為蒸散發(fā)能力；EU、EL、ED分別為上層、下層、深層蒸散發(fā)量；E為實際蒸散發(fā)量；WU、WL、WD分別為上層、下層、深層張力水蓄量；R為蓄滿產(chǎn)流計算的總徑流量；RS、RI、RG分別為地表徑流、壤中流、地下徑流；I和Q分別為馬斯京根河道演算的入流和出流.

圖1 三水源新安江模型結(jié)構(gòu)圖

3 實時修正方法

3.1 產(chǎn)流誤差的動態(tài)系統(tǒng)響應曲線修正方法

產(chǎn)流誤差的動態(tài)系統(tǒng)響應曲線修正方法是將計算值與實測值之間的誤差（流量誤差）歸咎于產(chǎn)流誤差，將圖1中蓄滿產(chǎn)流之后的部分概化為一個系統(tǒng)，假定流域出口斷面的流量過程變化是由于產(chǎn)流量R的變化所引起的系統(tǒng)響應，即

式中，Q（t）為流域出口的總徑流過程；R＝［r1，r2，…，rm］T為產(chǎn)流量系列；θ為模型參數(shù)；t為時間.式（1）的微分表達式為

式中，Rc＝［rc1，rc2，…，rcm］T為計算得到的產(chǎn)流量初始系列值.

假設樣本系列長度為m，則實測流量為Q（t）＝［Q1，Q2，…，Qm］T；計算流量為Q（Rc，θ，t）＝［Q（Rc，θ，1），Q（Rc，θ，2），…，Q（Rc，θ，m）］T，由式（2）可得：

式中，ΔR＝［Δr1，Δr2，…，Δrn］T為產(chǎn)流量誤差；E＝［ξ1，ξ2，…，ξm］T為經(jīng)實時修正后的預報系統(tǒng)殘差，一般為白噪聲向量；U為動態(tài)系統(tǒng)響應曲線矩陣，其表達式為

式（4）U中的每一項都可以用下式進行差分求解：

對式（4）進行最小二乘估計可得：

預報結(jié)果的修正式為

式中，為修正后產(chǎn)流量.

3.2 AR模型

AR模型假定各時段的流量誤差存在相關(guān)性，可用P階的線性關(guān)系來進行模擬.對于誤差系列｛ε1，ε2，…，εt，…，εt＋L，…｝，可建立式（8）表達式.

預報結(jié)果的修正式為

假設樣本系列長度為m，實測流量為Q（t）＝［Q1，Q2，…，Qm］T；計算流量為Qc（t）＝［Qc1，Qc2，…，Qcm］T，可得模型誤差系列｛ε1，ε2，…，εm｝，由式（8）可得矩陣形式為

式中：

對式（10）進行最小二乘估計可得：

4 應用檢驗

在將實時修正方法應用于實際流域后，修正效果用如下3項指標來比較衡量.

1）洪峰相對誤差

2）次洪徑流深相對誤差：

3）確定性系數(shù)：

本文選用東圳流域1965～1999年共16場歷史洪水的降雨、蒸發(fā)和流量資料，利用三水源新安江模型進行計算，模型參數(shù)見表1.所有參數(shù)在計算和修正過程中保持不變.

表1 東圳流域使用的三水源新安江模型參數(shù)

其中KC為蒸散發(fā)折算系數(shù)；WM為平均張力水容量；UM、LM分別為上層、下層張力水容量；B為張力水蓄水容量曲線方次；C為深層蒸散發(fā)折算系數(shù)；SM為表層自由水蓄水容量；EX為表層自由水蓄水容量曲線方次；KI、KG分別為表層自由水蓄水庫對壤中流、地下水的日出流系數(shù)；CS、CI、CG分別為地表水、壤中流、地下水消退系數(shù)；KE、XE為馬斯京根法演算參數(shù)［2］.

表2 新安江模型計算結(jié)果、AR模型和DSRC實時修正結(jié)果

表3 新安江模型計算結(jié)果、AR模型和DSRC實時修正結(jié)果比較

5 結(jié) 論

準確實時的洪水預報是東圳流域防洪安全的重要保證，實時修正模型能夠有效提高概念性水文模型的預報精度.本文分別應用DSRC和AR模型對東圳流域的預報結(jié)果進行實時修正，分析結(jié)果表明，相比于AR模型，DSRC對于洪峰的修正效果更佳，對徑流深的修正效果也較好，這對于水庫流域防洪抗災的意義重大，值得進一步地研究和推廣應用.

［1］ 劉國緯.論防洪減災非工程措施的定義與分類［J］.水科學進展，2003，14（1）：98－103.

［2］ 包為民.水文預報［M］.北京：中國水利水電出版社，2006.

［3］ 包為民，王 浩，趙 超，等.AR模型參數(shù)的抗差估計研究［J］.河海大學學報：自然科學版，2006，34（3）：258－261.

［4］ 司 偉，包為民，瞿思敏.洪水預報產(chǎn)流誤差的動態(tài)系統(tǒng)響應曲線修正方法［J］.水科學進展，2013，24（4）：497－503.

［5］ 郭 磊，趙英林.基于誤差自回歸的洪水實時預報校正算法的研究［J］.水電能源科學，2002，20（3）：25－27.

［6］ 瞿思敏，包為民.實時洪水預報綜合修正方法初探［J］.水科學進展，2003，14（2）：167－171.

［7］ 劉金英.東圳水庫生態(tài)環(huán)境與水資源保護的探討［J］.水利科技，2003，2（2）：14－15.