基于微分響應的流域產(chǎn)流分單元修正方法*

張小琴，吳成城，佘亮亮，包為民

(1:河海大學水文水資源學院，南京 210098)(2:寧波弘泰水利信息科技有限公司，寧波 315000)

洪水預報誤差來源多且影響復雜，洪水預報至今仍是水情預報的難點[1-2]. 實時修正技術(shù)是保證實時洪水預報精度的關(guān)鍵[3]. 現(xiàn)有誤差修正方法大多沒有考慮流域產(chǎn)匯流過程的物理基礎(chǔ)，主要以控制斷面誤差為基礎(chǔ)，將預報對象的誤差與模型分離，本質(zhì)上是“黑箱”分析方法，普遍存在外延性差、預見期損失等問題[4-5]. 洪水預報誤差修正研究中，通常需要考慮模型自變量變化與模型計算結(jié)果誤差之間的關(guān)系. 函數(shù)微分中自變量增量與因變量增量之間的微分響應關(guān)系，能夠較好地描述復雜系統(tǒng)中各影響因素之間的變化關(guān)系[6-7]. 包為民等提出將微分響應關(guān)系用于洪水預報誤差修正研究，具有物理基礎(chǔ)強、不增加參數(shù)、不損失預見期等優(yōu)點[8]. 現(xiàn)有研究主要利用微分響應方法對流域面平均降雨[9-10]、面平均產(chǎn)流[11-13]、面平均狀態(tài)變量[14-16]等進行修正，沒有考慮洪水預報誤差的空間分布差異.

流域降水時空分布和下墊面分布復雜，流域產(chǎn)匯流過程呈高度非線性、時空變化大[17]. 目前洪水預報誤差研究基本上是對各種誤差源等權(quán)看待，對誤差時空變化考慮不足會降低誤差修正效果，影響洪水預報精度. 有研究利用雨量站網(wǎng)密度與洪水預報誤差分配比例之間的定量關(guān)系將洪水預報總誤差劃分為面雨量輸入誤差和模型參數(shù)誤差[18]，利用各單元計算流量與各單元計算流量和的比例對流域面平均雨量估計值進行各單元分配以減小各單元面雨量誤差[19]. 洪水預報誤差修正大多基于流域系統(tǒng)輸出信息，雖然該信息能夠反映計算過程中各階段誤差的最終結(jié)果[20]，但難以區(qū)分反映各階段誤差對洪水預報的影響[21]. 從流域系統(tǒng)輸出的誤差序列中辨析誤差貢獻，考慮誤差時空變化進行實時修正有利于提高洪水預報精度.

本文針對洪水預報空間誤差修正困難的問題，基于水文系統(tǒng)微分響應關(guān)系構(gòu)建了流域產(chǎn)流分單元修正方法，并比較分析了在不同預見期時流域產(chǎn)流分單元修正、流域面平均產(chǎn)流修正和自回歸修正的效果. 該研究對于流域水文要素時空誤差修正、洪水預報精度提高具有重要意義和實用價值.

1 方法介紹

1.1 水文系統(tǒng)微分響應修正方法

水文模型的輸入、狀態(tài)變量、參數(shù)等任一變量的改變，都會引起流域出口斷面流量的變化. 把流域水文模型概化為一個復雜函數(shù)，利用微分建立自變量增量與因變量增量之間的響應關(guān)系：

(1)

式中,Q為模型因變量(模型輸出)；X為模型自變量(輸入變量、狀態(tài)變量、參數(shù)等)；?f/?X為變量X的微分響應.

式(1)中自變量增量即為模型影響因素的誤差，因變量增量即為模型計算結(jié)果的誤差. 式(1)的矩陣形式為：

ΔQ=UΔX+ξ

(2)

采用正則化最小二乘法求解式(2)，得變量X的誤差估計量為[16]：

ΔX=(UTU+βI)-1UTΔQ

(3)

式中,I為單位矩陣,β為正則化系數(shù).

修正后的變量X′為：

X′=X+ΔX

(4)

流域水文系統(tǒng)是一個非線性系統(tǒng)，考慮非線性問題線性化的影響，采用逐步迫近進行反演求解[22]：

X(k+1)=X(k)+ΔX(k)

(5)

式中, 上標k為迫近次數(shù),X(k)為第k次迭代的變量估計值, ΔX(k)為第k次迭代的變量誤差估計值.

將修正后的變量重新輸入模型進行計算可得到修正后的模型結(jié)果. 基于水文系統(tǒng)微分響應關(guān)系建立誤差修正方法，利用了產(chǎn)匯流模型結(jié)構(gòu)的物理基礎(chǔ)，概念清晰，計算簡便.

1.2 流域產(chǎn)流分單元微分響應修正方法

流域水文模型中每個子過程可以概化為一個系統(tǒng)，多個密切聯(lián)系的子過程也可以概化為一個系統(tǒng). 新安江模型包括蒸散發(fā)計算、產(chǎn)流計算、分水源計算和匯流計算[23]，可以將蒸散發(fā)計算和產(chǎn)流計算作為一個系統(tǒng)(記為E_R系統(tǒng))，水源劃分和匯流計算作為一個系統(tǒng)(記為R_Q系統(tǒng)). 產(chǎn)流計算是對降雨在數(shù)量上的再分配，產(chǎn)流量是影響水量平衡的重要因素，在洪水預報中起著重要作用. 對于產(chǎn)流微分響應修正，可通過R_Q系統(tǒng)利用流域出口流量反演修正產(chǎn)流(圖1).

圖1 新安江模型產(chǎn)流修正示意圖(P為降雨、EM為水面蒸發(fā)、R為產(chǎn)流、Q為流域出口流量)Fig.1 Runoff correction strategy for the Xin’anjiang model

采用新安江模型進行流域洪水預報，通常需進行流域劃分. 假設(shè)流域劃分為s個單元流域，流域出口流量對各單元流域產(chǎn)流的響應可以表達為：

Q=f(Ri,t)(i=1～s)

(6)

式中,i為單元流域編碼,s為單元流域數(shù)目,Ri=[ri,1,…,ri,j,…,ri,n]T為第“i”單元流域產(chǎn)流.

考慮誤差空間分布，構(gòu)建流域產(chǎn)流分單元微分響應修正(又稱流域產(chǎn)流分塊修正或時空產(chǎn)流修正)，記為SDR，如圖2a所示，計算步驟如下：

圖2 基于微分響應的流域產(chǎn)流修正流程圖Fig.2 Scheme of runoff correction based on differential response

(1)計算j時刻各單元流域產(chǎn)流量Ri(i=1～s)，將i單元流域j時刻計算產(chǎn)流量ri,j加上1個單位(其他時刻其他單元流域產(chǎn)流量保持不變)代入R_Q系統(tǒng)計算得到新的流量系列，減去原來的流量系列，得到i單元流域U矩陣j列uj.

(2)在不同時刻(j=1～n)重復步驟(1)即可得到i單元流域U矩陣所有列.

(3)計算產(chǎn)流誤差估計值ΔR，根據(jù)i單元流域面積占全流域面積的比例將估計的產(chǎn)流誤差分配給i單元流域，得到i單元產(chǎn)流誤差修正量為ΔR′i=ΔR·ηi，其中ηi為i單元流域面積占全流域面積的比例.

(4)重復步驟(1)～(3)，可得到各單元流域U矩陣及相應的產(chǎn)流修正量.

已有研究采用流域面平均產(chǎn)流修正[11-13]，記為AMR，如圖2b所示，計算流程如下：

(1)計算j時刻各單元流域產(chǎn)流量Ri(i=1～s)，將各單元流域j時刻計算產(chǎn)流量加上1個單位(其他時刻產(chǎn)流量保持不變)代入R_Q系統(tǒng)計算得到新的流量系列，減去原來的流量系列，得到全流域U矩陣j列uj.

(2)在不同時刻(j=1～n)重復步驟(1)就得到U所有列.

(3)計算產(chǎn)流量誤差估計值ΔR.

(4)根據(jù)單元流域面積占全流域面積的比例將產(chǎn)流誤差估計值分配給各單元流域，得到各單元產(chǎn)流誤差修正量ΔR′i.

流域產(chǎn)流分單元微分響應修正(SDR)和流域面平均產(chǎn)流微分響應修正(AMR)的不同點在于微分響應矩陣的構(gòu)建和產(chǎn)流估計誤差的分配. AMR在求解微分響應矩陣時所有單元流域產(chǎn)流量均增加1個單位，所求得的微分響應反映了流域出口流量對流域平均產(chǎn)流變化的響應. SDR分別計算不同單元流域?qū)奈⒎猪憫ㄟ^各單元流域微分響應提取流域出口斷面流量中包含的空間(單元流域)信息，增加了信息利用. AMR微分響應矩陣為2維數(shù)組Um×n，SDR微分響應矩陣為3維數(shù)組Um×n×s. AMR是將產(chǎn)流估計誤差按面積比例分配給各個單元流域，而SDR是將利用單元流域微分響應求得的誤差估計值分配給相應的單元流域.

2 應用檢驗

2.1 實例概況

將構(gòu)建的微分響應誤差修正方法用于修正三水源新安江模型計算產(chǎn)流. (1)比較SDR和AMR的修正效果. (2)比較SDR、AMR與二階自回歸(AR)在采用不同預見期時的修正效果.

選擇具有不同雨量站點密度和不同洪水特點的淮河大坡嶺流域和閩江七里街流域為研究區(qū)域. 流域圖如圖3，流域信息如表1. 采用泰森多邊形法，將大坡嶺流域劃分為13個單元流域，將七里街流域劃分為41個單元流域，每個單元流域有一個雨量站.

圖3 研究區(qū)域流域圖Fig.3 Location of the study basins表1 流域信息Tab.1 Informations of two basins

采用三水源新安江模型進行流域洪水模擬. 大坡嶺流域采用1987-2009年14場洪水資料，從起漲到洪峰出現(xiàn)的平均歷時為36.3 h，其中10場用于新安江模型參數(shù)率定，4場用于檢驗. 七里街流域采用1988-1999年20場洪水資料，從起漲到洪峰出現(xiàn)的平均歷時為77.9 h，其中15場用于新安江模型參數(shù)率定，5場用于檢驗.

流域名稱流域面積/km2雨量站個數(shù)多年平均降雨量/mm雨量站控制面積/(km2/站)資料系列 歷史洪水場次大坡嶺 164013918126.1519872009年14七里街14787411851360.6619881999年20

微分響應修正采用逐步迫近求解，以誤差均方差為目標函數(shù)判斷是否進行下一次迭代. 通過微分響應修正時所使用實測流量資料長度設(shè)置了一系列不同預見期，預見期越長所用實測流量時段數(shù)越少. 采用徑流深相對誤差(RRE)、洪峰相對誤差(RPE)、誤差均方差(RMSE)和Nash-Sutcliffe系數(shù)(NSE)[24]評估計算結(jié)果.

2.2 結(jié)果分析

表2 大坡嶺流域洪水模擬結(jié)果Tab.2 Flood simulation results in the Dapoling Basin

表3 七里街流域洪水模擬結(jié)果Tab.3 Flood simulation results in the Qilijie Basin

圖4 大坡嶺流域20050829洪水過程線Fig.4 Hydrographs for the flood event 20050829 in the Dapoling Basin

圖5 大坡嶺流域不同預見期時洪水修正結(jié)果比較Fig.5 Mean evaluation measures with different lead times in the Dapoling Basin

圖6 七里街流域19930630洪水過程線Fig.6 Hydrographs for the flood event 19930630 in the Qilijie Basin

圖7 七里街流域不同預見期時洪水修正結(jié)果比較Fig.7 Mean evaluation measures with different lead times in the Qilijie Basin

2.3 討論

大坡嶺流域和七里街流域所采用的歷史洪水從起漲到洪峰出現(xiàn)的平均歷時差別較大，能用于誤差修正的實測流量個數(shù)不同，微分響應修正均獲得了較好的修正效果. 在七里街流域部分雨量站雨量資料缺乏，采用鄰近站資料替代，雨量資料時空誤差較大，通過模型計算引起的產(chǎn)流時空誤差也較大. 對于時空誤差較大的洪水，AMR修正效果有限，SDR修正效果較好. 實例中微分響應修正以誤差均方差作為判斷條件進行逐步迫近反演求解，SDR修正和AMR修正的平均迭代次數(shù)均小于10次. 經(jīng)誤差修正后，誤差均方差減小了，NSE系數(shù)提高了，但并不能保證所有洪水的洪峰和徑流深模擬精度都提高，可以考慮采用兼顧洪峰和洪水過程的多目標函數(shù)作為判斷條件.

微分響應修正是基于流域出口流量與被修正量之間的微分響應關(guān)系建立的反饋修正方法，理論上可用于修正模型任一自變量，實際應用中通常用于修正主要誤差. 在洪水預報中可用于修正降雨、產(chǎn)流、參數(shù)、狀態(tài)變量等，修正效果取決于被修正變量是否是主要誤差來源. 產(chǎn)流計算受諸多誤差影響，模型輸入、狀態(tài)變量等誤差經(jīng)蒸散發(fā)計算和產(chǎn)流計算會傳播累積影響計算產(chǎn)流量，產(chǎn)流修正很有必要. 在大坡嶺流域預見期小于2 h，在七里街流域預見期小于4 h時，AR修正效果優(yōu)于SDR，主要是因為AR修正是直接對所有誤差影響下的流域出口流量進行修正. 但AR修正效果取決于輸出結(jié)果誤差之間的相關(guān)性，隨著預見期增大修正效果明顯變差.

結(jié)合新安江模型的產(chǎn)流微分響應修正是將分水源和匯流部分作為響應系統(tǒng)，利用流域出口流量誤差信息反饋修正產(chǎn)流量，將修正后的產(chǎn)流重新輸入模型進行后續(xù)的分水源、坡地匯流和河道匯流計算，進而修正流域出口斷面流量. 流域出口斷面流量是所有誤差影響的最終結(jié)果，對其進行時空誤差辨識很難. SDR修正建立了各單元流域產(chǎn)流與流域出口斷面流量之間的微分響應關(guān)系，利用各單元流域微分響應關(guān)系可以從流域出口斷面流量中提取有效的空間信息用于產(chǎn)流量修正. 基于微分響應的流域分單元誤差修正思路可為分布式流域水文模型誤差修正研究提供一種有效途徑，進一步可通過水文系統(tǒng)各要素引起的微分響應關(guān)系研究水文要素時空變化，成果應用可以有效提高實時洪水預報精度.

3 結(jié)論

1)基于微分響應構(gòu)建了流域產(chǎn)流分單元修正方法(SDR)，利用匯流系統(tǒng)物理基礎(chǔ)建立了各單元流域產(chǎn)流誤差與流域出口流量誤差之間的響應關(guān)系，采用正則化最小二乘法結(jié)合逐步迫近求解產(chǎn)流修正估計量，對產(chǎn)流誤差修正量進行時空分配實現(xiàn)流域產(chǎn)流分單元修正.

2)在淮河大坡嶺流域和閩江七里街流域結(jié)合新安江模型進行方法檢驗，結(jié)果表明：采用產(chǎn)流微分響應修正后的模型效果優(yōu)于原水文模型效果；SDR效果總體優(yōu)于AMR；隨著預見期增大，SDR和AMR的修正效果優(yōu)于AR.

3)微分響應分單元修正能通過水文系統(tǒng)從流域出口信息中提取有效的空間信息用于誤差修正，考慮了誤差的空間分布，使用簡單，不改變原模型結(jié)構(gòu)，可與不同水文模型結(jié)合用于修正不同時空誤差，對于時空誤差研究具有重要意義.