劉可新，李 匡，王紫微，梁犁麗，徐海卿

(中國水利水電科學研究院北京中水科水電科技開發(fā)有限公司，北京 100038)

1 研究背景

水文模型主要研究的是時不變的離線系統(tǒng)，大都采用觀測到的歷史水文資料，先確定好模型參數(shù)，然后用于未來的洪水預報。這種方式在實時洪水預報系統(tǒng)中時常得不到滿意的結果[1]。這時往往要采用誤差修正以取得較好的預報結果，因此誤差修正技術一直是洪水預報研究的重點，其中包括傳統(tǒng)的自回歸模型及其改進方法[2,3]、卡爾曼濾波技術[4]，后來的集合卡爾曼濾波[5]、神經網絡技術等[6]。近些年來，諸多學者在這方面做了大量工作，相繼提出了抗差分析[7-8]、綜合修正方法[9]等多種修正方法，均獲得了較好的修正效果。然而，傳統(tǒng)的誤差修正方法基本上都存在物理基礎不強、損失預見期和修正效果不理想等缺陷[10]。鑒于此，文獻[10]應用系統(tǒng)響應理論，將單位線作為線性時不變系統(tǒng)，引入實時洪水預報修正中，建立了一種向信息源頭追溯的反饋修正模型，繼文獻[10]之后，文獻[11]將動態(tài)系統(tǒng)響應曲線引入誤差修正中，提出了產流誤差動態(tài)系統(tǒng)響應曲線(Dynamic System Response Curve, DSRC)修正方法(以下簡稱DSRC方法)，對產流量進行修正，該方法具有物理基礎且不損失預見期，應用于實際流域，效果十分顯著，隨后系統(tǒng)響應修正方法得以發(fā)展與應用[12-17]。

但DSRC方法有時不穩(wěn)定，沒有考慮到流量測量誤差對產流修正穩(wěn)定性的影響，忽視了產流修正系列對流量測量誤差的敏感性，因而修正后的產流系列時常出現(xiàn)“震蕩”現(xiàn)象，本質上是反演計算的不適定問題。鑒于此，本文深入分析了流量測量誤差對產流修正系列的影響程度，通過正則化方法大大降低了產流修正系列對流量測量誤差的敏感性，提出了產流動態(tài)系統(tǒng)響應正則化(DSRR)修正方法(以下簡稱DSRR方法) 。

2 修正方法介紹

2.1 DSRC方法討論

流域出口的實測流量包含信息最全面，計算過程中各階段的誤差最終都將反映在流量上，因此，提取這些誤差并追溯其產生源頭是誤差修正的關鍵所在[18]。DSRC方法假設誤差來源于產流以上部分，以產流系統(tǒng)響應曲線作為誤差追溯其源頭的橋梁，從而將誤差合理分配到各時段產流并進行修正。

如前所述，DSRC方法是應用中的反問題，這類問題往往是不適定的[19]。所謂適定性是指解的存在性，解的唯一性和解的穩(wěn)定性。應用檢驗表明DSRC方法往往存在不穩(wěn)定問題，這一定程度上影響了修正效果。所謂穩(wěn)定性是指解對觀測值的連續(xù)依賴性，下面通過簡單例子解釋該性質。

線性方程組：

現(xiàn)在考慮觀測向量的微小誤差對解向量的影響，考察線性方程組：

解得x=(1 1)T。

從以上例子看出，觀測向量出現(xiàn)微小變化時，解向量卻變化很大，這就是解不穩(wěn)定的現(xiàn)象，其中方程組的系數(shù)矩陣A稱為病態(tài)矩陣，病態(tài)矩陣的條件數(shù)往往很大(A的條件數(shù)為40 002)。反演問題中往往存在病態(tài)矩陣，并且觀測向量的誤差難以避免，因此即使模型合理有效，也很難解得滿意的解向量。

DSRC方法同樣存在相應穩(wěn)定性問題，解向量對應于產流修正系列，觀測向量對應于流量誤差系列，實測流量存在誤差是難以避免的，因此流量誤差系列并非真值，存在一定的波動，于是產流誤差系列往往得不到滿意的結果，其不穩(wěn)定表現(xiàn)為產流修正系列時程上的“振蕩”現(xiàn)象。DSRC方法對產流修正系列的穩(wěn)定性要求較高，因為誤差修正是外延過程，產流修正系列不穩(wěn)定會嚴重影響修正效果。鑒此，提出了產流動態(tài)系統(tǒng)響應正則化修正方法，將不穩(wěn)定的反演問題轉換成鄰近的穩(wěn)定的反演問題，降低了產流修正系列對流量誤差系列的敏感程度，提高了產流修正系列的穩(wěn)定性，從而優(yōu)化了誤差修正效果。

2.2 DSRR方法介紹

近些年來，正則化方法已有較多的研究成果[20,21]。所謂正則化方法，是指用一些與原問題相鄰近的適定問題的解去逼近原問題的解，因此如何建立“鄰近”的問題是正則化方法的核心。Tikhcmov正則化方法是較成熟的一種正則化方法，其主要思想是通過增加一個罰函數(shù)項構建適定的鄰近問題。DSRC方法在有些情況下存在穩(wěn)定性問題，此時流量觀測誤差對產流修正系列影響較嚴重，因此本次研究將Tikhcmov正則化方法應用于DSRC方法中，提出了產流動態(tài)系統(tǒng)響應正則化修正方法。

為研究方便，以新安江模型產流以下的部分作為系統(tǒng)，產流系列作為系統(tǒng)輸入，出口流量作為系統(tǒng)輸出，系統(tǒng)結構如圖1。

圖1 系統(tǒng)示意

圖1是新安江模型示意圖，其中E0表示實際蒸散發(fā)，P表示降雨，R表示產流量，Q表示出口流量，TLOC表示三層蒸發(fā)，NSM表示蓄滿產流，SRC表示分水源，F(xiàn)T表示匯流，BC表示反演計算，虛線框中的部分表示系統(tǒng)。

根據(jù)文獻[11]，只考慮產流對出口流量的影響，那么圖1所示系統(tǒng)可以表示為：

Q(t)=Q(R,θ,t)

(1)

式中：R=[r1,r2,…,rn]T表示產流系列；Q(R,θ,t)即實測流量系列Q0。

根據(jù)動態(tài)系統(tǒng)響應理論，建立產流與流量之間的關系，可用式(2)表示：

Q(R,θ,t)=Q(RC,θ,t)+UΔR+E

(2)

式中：RC=[rC1,rC2,…,RCn]T為計算得到的產流量初始系列；ΔR=[Δr1,Δr2,Δr3,…]T為待估計的產流修正系列；Q(RC,θ,t)為計算流量初始系列；E=[e1,e2,…,eL]T為流量觀測隨機誤差項；U為動態(tài)系統(tǒng)響應矩陣，其表示各時段產流量與出口流量之間的關系。

DSRC方法根據(jù)無約束最小二乘原理，所求ΔR滿足式(3)：

Q(RC,θ,t)-UΔR]

(3)

于是有：

ΔR=(UTU)-1UT[Q(R,θ,t)-Q(RC,θ,t)]

(4)

實踐發(fā)現(xiàn)，UTU矩陣有時十分病態(tài)甚至趨于奇異，這時流量觀測誤差對ΔR的影響很大，甚至導致估計嚴重失真。

鑒于此，引進了Tikhcmov正則化方法，提出了DSRR方法，在式(3)中加入罰函數(shù)項，如式(5)：

(5)

式中:α為權重系數(shù)。

經推導得式(6)：

ΔR=(UTU+αI)-1UT[Q(R,θ,t)-Q(RC,θ,t)]

(6)

式中：I為單位矩陣。

DSRR方法較DSRC方法更穩(wěn)定表現(xiàn)在兩個方面，一方面式(5)是約束最小二乘法的目標函數(shù)，由兩部分組成，誤差項和罰函數(shù)項，其中誤差項表示修正后計算流量與實測流量的接近程度，其值越小表示越接近，罰函數(shù)項表示穩(wěn)定程度，若要滿足其值較小，那么ΔR不會出現(xiàn)“振蕩”現(xiàn)象；另一方面從式(6)中可以看出，當UTU病態(tài)甚至奇異時，加上αI項可解決該問題，降低ΔR對流量觀測誤差的敏感程度，消除估計嚴重失真的現(xiàn)象。α作為權重系數(shù)，在式(5)中作用是平衡誤差項與罰函數(shù)項，在式(6)中作用是保證降低UTU病態(tài)程度，而又滿足與原問題鄰近。

修正后計算產流系列R′C為：

R′C=RC+ΔR

(7)

將R′C代入模型重新計算出流過程：

Q′=Q(R′C,θ,t)

(8)

Q′為修正后的計算流量系列。

3 實際流域應用與檢驗

3.1 評定指標的選取

為了較好地反映修正效果及修正的穩(wěn)定性，本文采用以下4個指標作為評價指標，其中洪量相對誤差、洪峰相對誤差和Ns系數(shù)反應修正效果，條件數(shù)反應修正的穩(wěn)定性。各指標的計算公式如下：

(1)洪量相對誤差：

(9)

(2)洪峰相對誤差：

(10)

(3)Nash-Sutcliffe系數(shù)：

(11)

(4)條件數(shù)：

cond(A)2=‖A‖2‖A-1‖2

(12)

3.2 實際案例檢驗

本次研究將灘坑流域作為實驗流域，灘坑流域位于浙江省甌江流域的小溪支流上，流域面積為3 330 km2，流域內有15個雨量站。灘坑流域受亞熱帶季風影響，氣候溫和，雨量豐沛，四季變化明顯，日照充足。流域徑流以降水補給為主，降水量年際變化大，豐枯水年相互交替出現(xiàn)，流域年平均降水量一般在1 500～2 100 mm，多年平均水面蒸發(fā)量為969.9 mm。流域內降水量空間分布不均，年分配有顯著差異，全年大致可分為3-4月春雨，5-6月梅雨，7-9月臺風雷雨和10-次年2月干季4個時期。流域多年平均流量120 m3/s，多年平均年徑流量37.8 億m3，徑流年內分配不均勻，全年以6月份來水量最大，占全年水量的19.7%。流域水系及站網布設見圖2。

圖2 灘坑流域示意圖

本次研究以灘坑流域1985—2004年19場洪水作為檢驗資料，使用新安江模型對流域的歷史洪水進行計算和修正。為全面深入分析DSRR方法的修正效果，選取其中一場洪水(洪號：31980618)作為典型，分別采用DSRC方法和DSRR方法，以不損失預見期為條件，對該場洪水進行修正，其效果見圖3和圖4。由圖3看出，DSRR方法修正后的時段產流平穩(wěn)變化，沒有出現(xiàn)陡漲陡落的情況，并且修正前后產流過程線形狀基本保持一致，說明修正效果符合實際情況；而應用DSRC方法，修正后時段產流雖然前期比較平穩(wěn)，但后期出現(xiàn)了較明顯的不穩(wěn)定現(xiàn)象，而且隨著時段的增加，不穩(wěn)定越來越嚴重，最后甚至出現(xiàn)了“振蕩”現(xiàn)象，表現(xiàn)為相鄰時段產流量交替變化，修正效果明顯不符合產流規(guī)律，脫離實際情形。出現(xiàn)這一狀況的原因考慮有以下幾點：一方面矩陣UTU條件數(shù)(1.2×109)過大，是病態(tài)矩陣，使得產流修正系列ΔR對流量觀測誤差過于敏感，這就意味著微小的觀測誤差同樣會造成ΔR劇烈變化，而流量觀測誤差是難以避免的，于是采用DSRC方法會導致震蕩現(xiàn)象，而采用DSRR方法，矩陣UTU+αI的條件數(shù)(780)大大減小了，這就降低了ΔR對流量觀測誤差的敏感性，增強了方法的穩(wěn)定性；另一方面DSRR方法中加入了罰函數(shù)項，約束產流修正系列，使其在一定范圍內變化；此外，以不損失預見期為條件，使得可用的實測流量信息有限，有效信息不足以彌補由于測量誤差產生的干擾信息的影響。由4中不難看出，DSRR方法得到的洪水過程線比較光滑，而DSRC方法則出現(xiàn)了局部的陡漲陡落，這與兩種方法產流修正效果的穩(wěn)定程度一致；對洪水的修正效果，DSRR方法整體上優(yōu)于DSRC方法，這是由于新方法使產流修正更加穩(wěn)定，修正結果更加符合實際情況，而DSRC方法雖注重洪水過程線擬合程度，但忽視了流量測量誤差對產流修正結果的嚴重影響。

圖3 DSRR方法與DSRC方法產流修正效果比較

圖4 DSRR方法與DSRC方法典型洪水修正效果對比

灘坑流域19場洪水模擬及修正結果見表1。由表1看出，DSRR方法的修正效果十分顯著，各場洪水的預報精度明顯提高，平均洪量相對誤差、洪峰相對誤差和峰現(xiàn)時差分別降低了3.8%、8.2%和1.0h，平均確定性系數(shù)提高到0.970；DSRR方法的應用效果明顯優(yōu)于DSRC方法，各項指標均有所改善；同時發(fā)現(xiàn)采用DSRC方法，每場洪水條件數(shù)均很大，平均達到1.1×109，說明該方法修正結果的穩(wěn)定性較差，而DSRR方法將條件數(shù)降低了7個數(shù)量級，平均為773，極大地降低了產流修正系列對流量測量誤差的敏感程度，提高了修正結果的穩(wěn)定性；此外還發(fā)現(xiàn)，由于各場洪水的降雨蒸發(fā)等水文條件不同，采用DSRR方法使修正結果趨于穩(wěn)定時，不同洪水其條件數(shù)各異。

表1 灘坑流域19場洪水實時修正結果

本文同時對權重系數(shù)α做了初步探討，在不損失預見期條件下，繪制了不同洪水α與產流修正個數(shù)的關系圖(圖5)，同時在產流修正個數(shù)逐漸增加條件下，繪制了同一場洪水(洪號：31980618)α的變化趨勢圖(圖6)。由圖5看出，不同洪水的α與產流修正個數(shù)具有一定的變化一致性，但趨勢并不顯著，也就是說對于不同的洪水， 的取值與產流修正個數(shù)有一定正相關關系，但還要視情況而定，這主要是由于不同洪水水文條件各異所致；由圖6看出，同一場水，α隨產流修正個數(shù)增加而變大，趨勢十分顯著。

4 結 語

本文深入分析了DSRC方法，提出了DSRR方法，對DSRR方法中α的取值做了初步探討，并對兩種方法的修正效果進行了比較，得出以下幾點結論：

(1)DSRC方法在有些情況下，產流修正系列對流量測量誤差的敏感性較強，難以保證修正結果的穩(wěn)定性；

(2)產流修正系列對流量測量誤差的敏感性可用條件數(shù)量化表示，條件數(shù)越大敏感性往往越大；

(3)DSRR方法降低了產流修正系列對流量測量誤差的敏感性，其修正結果的穩(wěn)定性及修正效果明顯優(yōu)于DSRC方法；DSRR方法中α的取值視洪水而定，不同洪水產流修正個數(shù)對α的取值有一定的參考作用，同一場洪水，隨著產流修正時段的增加，α逐漸變大。

□

[1] 田 雨, 雷曉輝, 蔣云鐘, 等. 洪水預報實時校正技術研究綜述[J]. 人民黃河, 2011,33(3):25-28.

[2] 包為民. 水文預報[M]. 北京: 中國水利水電出版社, 2006.

[3] 周 軼, 菅浩然, 李致家, 等. 基于遞推最小二乘改進算法的洪水預報模型研究[J]. 河海大學學報,2007,35(1):77-80.

[4] 葛守西, 程海云, 李玉榮. 水動力學模型卡爾曼濾波實時校正技術[J]. 水利學報, 2005,36(6):1-9.

[5] CLARK M P, RUPP D E,WOODS R A, et al. Hydrological data assimilation with the ensemble Kalman filter: use of streamflow observation to update states in a distributed hydrological model[J]. Advances in Water Rsources, 2008,31:1 309-1 324.

[6] 朱星明, 盧長娜, 王如云, 等.基于人工神經網絡的洪水水位預報模型[J].水利學報,2005,36(7):806-811.

[7] 包為民, 嵇海祥, 胡其美, 等. 抗差理論及在水文學中的應用[J]. 水科學進展, 2003,14(4):528- 532.

[8] 趙 超, 洪華生, 包為民, 等. 實時洪水抗差預報系統(tǒng)研究[J]. 水文, 2008,28(2):27-29.

[9] 瞿思敏, 包為民. 實時洪水預報綜合修正方法初探[J]. 水科學進展, 2003,14(2):167-171.

[10] 包為民, 司 偉, 沈國華, 等. 基于單位線反演的產流誤差修正[J]. 水科學進展, 2012,23(3):315-322.

[11] 司 偉, 包為民, 瞿思敏. 洪水預報產流誤差的動態(tài)系統(tǒng)響應曲線修正方法[J]. 水科學進展, 2013,24(4):497-503.

[12] Wei Si, Weimin Bao, and Hoshin V. Gupta. Updating real-time flood forecasts via the dynamic system response curve method[J]. Water Resources Research, 2015,(7).

[13] 張小琴, 劉可新, 包為民, 等. 產流誤差比例系數(shù)的系統(tǒng)響應修正方法[J]. 水科學進展, 2014,25(6):789-796.

[14] 包為民, 闕家駿, 賴善證, 等. 洪水預報自由水蓄量動態(tài)系統(tǒng)響應修正方法[J]. 水科學進展, 2015,26(3):365-371.

[15] 楊姍姍, 包為民, 楊小強, 等.微分響應在降雨誤差修正中的應用[J]. 中國農村水利水電, 2015,(1):75-79.

[16] 劉可新, 包為民, 李佳佳, 等. 基于系統(tǒng)響應理論的分水源誤差修正[J]. 水電能源科學, 2014,32(11):44-47.

[17] 劉可新, 包為民, 賴善證, 等. 動態(tài)系統(tǒng)響應曲線修正方法在烏溪溝流域的應用[J]. 中國農村水利水電, 2014,(12):24-26.

[18] 包為民. 洪水預報信息利用問題與討論[J]. 水文, 2006, 26(2):18-21.

[19] 顧勇為, 歸慶明, 邊少鋒, 等. 地球物理反問題中兩種正則化方法的比較[J]. 武漢大學學報·信息科學版, 2005,30(3):238-241.

[20] 董國志. 反問題的正則化方法及其計算[D]. 長沙: 湖南師范大學，2012.

[21] 楊 帆.三類不適定問題的正則化方法研究[D].蘭州:蘭州大學，2014.