面平均雨量誤差和產(chǎn)流量誤差修正的對比研究
——以富江流域為例

司 偉，余鴻慧，包為民，瞿思敏，張 乾

(1. 河海大學(xué)水文及水資源學(xué)院，南京 210098；2. 囯網(wǎng)新源水電有限公司富春江水力發(fā)電廠，杭州 311504)

降雨徑流過程涉及蒸發(fā)、下滲、產(chǎn)流和匯流等多個環(huán)節(jié)，影響因素復(fù)雜，實時洪水預(yù)報精度難以滿足要求。預(yù)報模型初始狀態(tài)值誤差、模型輸入誤差、模型結(jié)構(gòu)的不充分性以及參數(shù)非優(yōu)化都會影響預(yù)報結(jié)果的精度[1,2]。所以經(jīng)常使用誤差修正技術(shù)來減少模型在預(yù)報時產(chǎn)生的誤差[3]。

簡單的誤差修正方法是通過流量誤差直接對模型計算結(jié)果(流量過程)進(jìn)行修正[4]。組合卡爾曼濾波(EnKF)或者粒子濾波(PF)方法經(jīng)常通過數(shù)據(jù)同化技術(shù)對于模型初始值或者模型參數(shù)進(jìn)行修正，然而研究結(jié)果表明濾波方法在預(yù)報修正中的效果過于自信，對于模型預(yù)報的不確定性估計不足[5]。

鑒于此，包為民等[6]應(yīng)用系統(tǒng)響應(yīng)理論，首次將單位線作為線性時不變系統(tǒng)，引入實時洪水預(yù)報修正中，建立了一種向信息源頭追溯的反饋修正模型，該方法結(jié)構(gòu)簡單，物理概念清晰。之后司偉等[7]將動態(tài)系統(tǒng)響應(yīng)曲線引入誤差修正中，提出了產(chǎn)流誤差動態(tài)系統(tǒng)響應(yīng)曲線(Dynamic System Response Curve, DSRC) 修正方法，基于系統(tǒng)和微分角度，逐時段對產(chǎn)流量進(jìn)行修正，該方法具有物理基礎(chǔ)且不損失預(yù)見期，通過不同流域的應(yīng)用檢驗[8,9]，效果顯著，明顯優(yōu)于傳統(tǒng)的自回歸修正方法。在隨后的研究中，張小琴、劉可新等對此方法進(jìn)行了改進(jìn)[10,11]，楊姍姍等將此理論方法引入對面平均雨量誤差的修正研究中，通過修正水文模型的重要輸入項(面平均雨量)誤差，來提高水文模型的洪水預(yù)報精度，研究結(jié)果表明此方法對于雨量站密度較低的流域具有顯著的效果[12]。

以上基于系統(tǒng)響應(yīng)理論的洪水預(yù)報誤差修正研究主要分為對面平均雨量誤差進(jìn)行修正和對產(chǎn)流量誤差進(jìn)行修正，兩種方法對于提高洪水預(yù)報精度均有一定效果。然而，以上研究并未對兩種方法在同一流域的應(yīng)用效果進(jìn)行探討，鑒于此，本文作者將分別對兩種修正方法進(jìn)行簡單描述，并對兩種方法應(yīng)用于富春江流域的結(jié)果進(jìn)行對比和分析，希望對于研究人員選取合適的修正方法有一定的指導(dǎo)意義。

1 修正方法介紹

本文中兩種修正方法(面平均雨量誤差修正和產(chǎn)流量誤差修正)均基于水文模型的系統(tǒng)響應(yīng)理論，即將水文模型概化為輸入與輸出之間的響應(yīng)系統(tǒng)，根據(jù)修正內(nèi)容的不同來改變響應(yīng)系統(tǒng)的結(jié)構(gòu)。

1.1 面平均雨量誤差修正方法

如圖1所示，將水文模型作為一般系統(tǒng)考慮，圖中P(t)為面平均降雨系列，E(t)為蒸散發(fā)，Q(t)為計算的流量過程，θ為模型參數(shù)，X(t)是模型狀態(tài)變量，X0為狀態(tài)初始值。圖1中的水文模型系統(tǒng)可以表達(dá)為由輸入、輸出、參數(shù)和狀態(tài)變量組成的系統(tǒng)方程：

Q(t)=f[P(t),E(t),X(t),θ]

(1)

式中：Q=[q1,q2,q3,…,qn]T為流域出口斷面流量過程；P=[p1,p2,p3,…,pm]T為面平均雨量向量；E=[e1,e2,e3,…,en]T為蒸發(fā)；X=[x1,x2,x3,…,xn]T為狀態(tài)變量向量；θ為參數(shù)向量。

圖1 一般水文模型系統(tǒng)示意圖Fig.1 A systems diagram of general hydrological model

概念性模型的絕大多數(shù)參數(shù)都具有明確的物理意義，是流域水文特征的反應(yīng)，模型參數(shù)一旦率定之后一般在短期之內(nèi)不會改變。模型狀態(tài)變量取決于模型的輸入輸出項[13]，因此，此方法研究中認(rèn)為模型參數(shù)不變化，面平均雨量的變化是引起流域出口斷面流量變化的主要因素。上述系統(tǒng)方程可以簡化為：

Q(P)=f(P)

(2)

對于流域水文模型，把降雨輸入作為自變量求全微分，可得到式(2)的微分關(guān)系：

(3)

(4)

式中：Q(PO)為由降雨系列PO計算得到的流量系列。假設(shè)使用觀測流量系列長度為L，Q(t)=[Q1,Q2,Q3,…,QL]T。代入式(4)可得：

(5)

則式(5)的矩陣形式為：

Q(P)=Q(PO)+SΔP+W

(6)

式中：ΔP=[Δp1,Δp2,Δp3,…,Δpm]T即為面平均雨量誤差的估計值；W=[w1,w2,w3,…,wL]T為流量觀測隨機(jī)誤差項，一般為服從零均值分布的為白噪聲向量；S為系統(tǒng)響應(yīng)矩陣，其表達(dá)式如(7)所示：

(7)

將式(7)代入(3)得到由降雨變化量與流量變化量表達(dá)的系統(tǒng)響應(yīng)關(guān)系：

dQ=SdP

(8)

上式S矩陣中的每一項都可以用下式差分近似求解。

(9)

式中：Q(p1,…,pi,…)為由修正前的面平均雨量系列計算得到的流量系列；Q(p1,…,pi+Δpi,…)為由面平均雨量系列增加ΔPi之后計算得到的流量系列。

上式中當(dāng)i不變，而j從1到L變化時，這L項差分值就是S矩陣中的一列。而此列正是雨量pi的單位變化量所對應(yīng)的系統(tǒng)響應(yīng)系列，本研究中稱之為雨量pi所對應(yīng)的系統(tǒng)響應(yīng)曲線，S矩陣中的每一列均用同樣方法求得。

對于一場洪水，假如有m個時段雨量，其相應(yīng)的雨量誤差依時間順序為eP1,eP2,…,ePm，假如有L個時段的流量系列，流量的誤差序列為eQ1,eQ2,…,eQL，則降雨誤差與流量誤差之間的關(guān)系可以表達(dá)為式(10)：

(10)

上式表達(dá)為向量矩陣形式為：

EQ=S·EP

(11)

式(11)表達(dá)了降雨誤差與流量計算誤差的響應(yīng)關(guān)系。則面平均雨量誤差與流量誤差之間的系統(tǒng)響應(yīng)關(guān)系為：

EQ=S·EP+ζ

(12)

面平均雨量誤差的最小二乘估計值為：

E′P=(STS)-1STEQ

(13)

式中：ζ為殘差向量；E′P為降雨誤差的動態(tài)系統(tǒng)響應(yīng)曲線估計向量。由式(13)計算得到的面平均雨量誤差系列加上原來的面平均雨量系列重新使用預(yù)報模型計算，便可得到通過降雨系統(tǒng)響應(yīng)曲線反饋修正以后的流域出口斷面流量過程。

1.2 產(chǎn)流量誤差修正方法

此修正方法將預(yù)報模型的匯流部分作為響應(yīng)系統(tǒng)，如圖2中虛線框內(nèi)系統(tǒng)，通過計算時段產(chǎn)流量所對應(yīng)的系統(tǒng)響應(yīng)曲線來修正產(chǎn)流量，進(jìn)而修正出口斷面流量過程。圖2中方框內(nèi)SERM表示蓄滿產(chǎn)流，TLSE表示三層蒸散發(fā)，SORC表示分水源，F(xiàn)C表示匯流，MMSSR表示馬斯京根河道分?jǐn)噙B續(xù)演算。E為蒸散發(fā)，P為降雨，R為產(chǎn)流量，RS、RI和RG分別為劃分以后的3種徑流成分，Q為坡地匯流后的流量，QT為流域出口的總徑流過程。

圖2 新安江模型概化系統(tǒng)示意圖Fig.2 The generalized system diagram of XAJ model

則產(chǎn)流量和模型預(yù)報流量之間的系統(tǒng)響應(yīng)函數(shù)可以表示為：

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

(14)

與面平均雨量誤差修正方法類似，由式(13)可得到產(chǎn)流量的修正量的計算式為：

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

(15)

式中：U為產(chǎn)流量系統(tǒng)響應(yīng)矩陣；RC為修正前計算產(chǎn)流量系統(tǒng)，則修正后的產(chǎn)流量系列為：

R′C=RC+ΔR

由式(16)計算出修正后的產(chǎn)流量R′C。 將修正后的產(chǎn)流量系列重新用新安江模型進(jìn)行計算，即可得到修正的流域出口斷面的計算流量過程QU。

2 應(yīng)用對比

2.1 流域介紹

富春江流域面積為18 866 km2，流域內(nèi)共有28個雨量站，多年平均降雨量為1 659.2 mm。其中蘭溪水文站是富春江干流上的重要入庫控制站，控制總流域面積的85%，控制范圍內(nèi)有24個雨量站，對于整個流域的洪水控制和管理起著重要作用。富春江水庫為調(diào)節(jié)能力較小的日調(diào)節(jié)水庫，承擔(dān)著防洪灌溉、發(fā)電、城市供水、航運(yùn)以及平衡下游河道水環(huán)境容量和應(yīng)對突發(fā)水質(zhì)事件等重要功能，因此入庫洪水預(yù)報的精度非常重要。富春江水庫上游的兩個重要支流分別是金華江和衢江，兩條主要支流上均有徑流式電站，其放水流量不規(guī)律，尤其是衢江小溪灘徑流式電站以上還有大型水庫，流域內(nèi)的水利工程設(shè)施給富春江水庫的入庫洪水預(yù)報造成很大的困難。由于富春江流域面積較大，而流域內(nèi)雨量站個數(shù)相對較少，雨量站密度較低(673.8 km2/個)，造成雨量站資料代表性不足。對于空間集總式或半分布式水文模型，流域的大小以及流域內(nèi)雨量站密度會對面平均雨量計算有較大影響。

2.2 應(yīng)用結(jié)果對比

鑒于富春江流域洪水預(yù)報誤差的兩個主要方面的影響，本文將兩種不同的預(yù)報誤差修正方法應(yīng)用于富春江流域，并對應(yīng)用結(jié)果進(jìn)行對比分析。本文使用富春江流域2010-2014年間的14場歷史洪水進(jìn)行了驗證。兩種方法修正后的計算結(jié)果如表1所示。

由表1的數(shù)據(jù)可知，經(jīng)過兩種方法的應(yīng)用14場歷史洪水的預(yù)報精度均有明顯提高。在面平均雨量誤差修正方法中，平均NS系數(shù)由修正前的0.769提高到0.848，提高了0.079；徑流深相對誤平均值由修正前的-1.5%減小到-0.6%，減小了0.9%；洪峰相對誤平均值由修正前的-1.1%減小到-0.3%，減小了0.8%；對于產(chǎn)流量誤差修正方法，平均NS系數(shù)由修正前的0.769提高到0.853，提高了0.084；徑流深相對誤平均值由修正前的-1.5%減小到-0.6%，減小了0.9%；洪峰相對誤平均值由修正前的-1.1%減小到-0.1%，減小了1.0%，結(jié)果顯示兩種修正方法均有明顯的效果。

為了分析不同量級洪水修正效果的代表性，本文中將歷史洪水分級。由表1可知，14場歷史洪水的平均次洪降雨量為206.9 mm，故本文中將降雨量大于平均次洪雨量的歷史洪水歸為大洪水(1～5號洪水)，其余歸為中小洪水(6～14號洪水)。驗證兩種方法對于不同量級洪水修正效果提高的幅度，本文中使用如下指標(biāo)：

δNS=(NS-M-NS-X)/(NS-C-NS-X)

(17)

式中：NS-M為面平均雨量修正方法的效果；NS-C為產(chǎn)流量修正方法的效果；NS-X為未經(jīng)修正的新安江模型計算結(jié)果。其值大于1則表明面平均雨量誤差的修正效果比產(chǎn)流量修正方法修正效果好，小于1則后者修正效果優(yōu)于前者。修正效果對比如圖3所示。

由圖3結(jié)果可以看出，對于5場大洪水，面平均雨量修正方法的效果比產(chǎn)流量誤差修正的效果好，由于富春江流域雨量站密度較低，降雨量較大時，采用雨量站單點雨量代替雨量站單元的面平均雨量進(jìn)行預(yù)報時，由雨量代表性誤差造成的預(yù)報誤差更加明顯和突出。對于中小洪水，產(chǎn)流量修正方法的效果比面平均雨量誤差修正方法效果好，是因為次洪雨量小，由雨量誤差引起的流量誤差相對減小，而產(chǎn)流量誤差中包含了其他方面誤差，因此產(chǎn)流量修正方法對于中小洪水的修正效果較好。但是從修正效果的平均值來看，產(chǎn)流量誤差修正方法的修正效果優(yōu)于面平均雨量誤差修正方法。

表1 富春江流域14場歷史洪水的實測值和修正前后的結(jié)果對比Tab.1 The observed values, simulated results and updated results of the 14 flood events in Fuchunjiang basin

注：表中P表示降雨，RO表示實測徑流深，QOP表示實測洪峰流量，ΔRT表示徑流深相對誤差，ΔQP表示洪峰流量相對誤差，NS表示Nash-Sutcliffe系數(shù)。

圖3 兩種方法修正效果對比Fig.3 Improvement comparison of the two updating method

3 結(jié) 語

本文將基于系統(tǒng)響應(yīng)曲線的面平均雨量誤差修正方法與產(chǎn)流量誤差修正方法應(yīng)用于富春江流域的相同歷史洪水中，通過對應(yīng)用結(jié)果的對比和分析，兩種方法修正后的平均NS系數(shù)分別由修正前的0.769提高到0.848和0.853，分別提高了0.079 和0.084，均具有較為明顯的修正效果。

然而通過不同量級洪水修正效果的對比得出，在對于大洪水進(jìn)行修正應(yīng)用時，面平均雨量修正方法的修正效果較好，對于中小洪水，產(chǎn)流量修正方法的效果優(yōu)于面平均雨量方法的效果。另外，面平均雨量修正方法適用于所有概念性水文模型，因此在次洪降雨量較大且雨量站密度較低的大流域推薦使用面平均雨量修正方法。在應(yīng)用于次洪降雨量較小的中小型流域時，尤其當(dāng)使用的水文模型為集總式產(chǎn)流模型(如新安江模型、HYMOD等)時，建議使用產(chǎn)流量修正方法。希望本文的研究結(jié)果對于研究者選取適當(dāng)?shù)膶崟r洪水預(yù)報誤差修正方法具有一定的指導(dǎo)意義。

□

[1] Komma J, Bloschl G, Reszler C. Soil moisture updating by ensemble kalman filtering in real-time flood forecasting[J]. Journal of Hydrology, 2008,357(3-4): 228-242.

[2] Lee W L, Eugene L P. Accuracy of precipitation measurements for hydrologic modeling[J]. Water Resources Research, 1974,10(4):857-863.

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

[4] Marco F, Anna B, Silvia B, et al. Forecasting discharges at the downstream end of a river reach through two simple Muskingum based procedures[J]. Journal of Hydrology, 2011,399:335-352.

[5] Moradkhani H, Dechant C M, Sorooshin S. Evolution of ensemble data assimilation for uncertainty quantification using the particle filter-markov Chain monte carlo method[J]. Water Resources Research, 2012,48(12), W12520, DOI:10.1029/2012WR012144.

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

[7] Bao W, Si W, Qu S. Flow updating in real-time flood forecasting based on runoff correction by dynamic system response curve[J]. Journal of Hydrologic Engineering, 2014,(19):747-756.

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

[9] 闕家駿, 包為民. 產(chǎn)流誤差的動態(tài)系統(tǒng)響應(yīng)曲線修正方法應(yīng)用[J]. 三峽大學(xué)學(xué)報(自然科學(xué)版), 2014,36(6):6-9.

[10] 張小琴, 劉可新, 包為民, 等. 產(chǎn)流誤差比例系數(shù)的系統(tǒng)響應(yīng)修正方法[J]. 水科學(xué)進(jìn)展, 2014,25(6):789-96.

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

[12] 楊姍姍, 包為民, 楊小強(qiáng), 等. 微分響應(yīng)在降雨誤差修正中的應(yīng)用[J]. 中國農(nóng)村水利水電, 2015,(1):75-9.

[13] Hamid M, Soroosh S, Hoshin G, et al. Dual state-parameter estimation of hydrological models using ensemble Kalman filter[J]. Advances in Water Resources, 2005,28:135-147.