包為民，顧雨薇，司 偉**，侯 露，盧金利，落全富

(1:河海大學水文水資源學院，南京 210098)(2:臺州市水文站，臺州 318001)(3:青山水庫管理處，杭州 311305)

黃土地區(qū)泥沙問題是目前水文研究的重點問題之一. 研究黃河水沙變化、把握黃河水沙變化規(guī)律，對于進一步完善治黃方略、實施流域水資源配置與管理及重大水利工程布局，意義重大[1]. 目前泥沙模型的研究取得了較為豐碩的成果，主要包括經(jīng)驗相關模型和物理成因模型. 經(jīng)驗相關模型中，Wischmeier WH等[2]提出了著名的通用土壤侵蝕方程USLE. 基于物理機理的流域水沙模型研究較多，湯立群等[3]基于坡面徑流侵蝕量計算公式提出了流域產(chǎn)沙模型；WEPP、LISEM、EUROSEM等模型[4-6]基于土壤侵蝕過程的物理描述，模擬了流域侵蝕過程，在許多國家被成功推廣應用. 已有的大多數(shù)流域產(chǎn)沙模型未能考慮不同流域尺度的泥沙規(guī)律，模型模擬精度以及使用范圍受到了很大限制. 因此，包為民[7]基于流域自然地理特性，在對黃土地區(qū)產(chǎn)沙機制進行概化的基礎上，提出了結(jié)構(gòu)和參數(shù)均具有物理意義的流域水沙模擬概念模型. 該模型經(jīng)黃土地區(qū)不同流域尺度的實際流域和前期大量的試驗流域驗證，可以很好地模擬黃土高原地區(qū)流域的水沙規(guī)律[8].

水沙模擬中重要的動力條件是水流，因此水流模擬的精度直接影響流域產(chǎn)沙和匯沙過程，進而影響水沙模型的模擬精度和效果[9]. 水流模擬的主要技術(shù)手段是依賴于水文模型. 水文模型計算誤差來源有很多，其中包括重要輸入項降雨資料、下墊面條件、流域初始狀態(tài)值、模型本身誤差等[10]. 由于降雨是水文預報模型計算的重要輸入項，因此降雨資料直接影響模型產(chǎn)流，進而影響模型的模擬精度[11]. 在水沙模型計算過程中，基于黃土區(qū)域自然地理特性，流域面平均雨量大多采用“以點帶面”的方式進行計算，即利用雨量站點的觀測資料代表該雨量站所代表的單元面上的降雨，這種方式不可避免的給流域降雨帶來誤差. 同時，部分降雨資料由于觀測等因素存在時段均化現(xiàn)象，伴有一定程度的缺測、漏測和誤測問題，降雨資料由插補展延而來, 造成模型模擬誤差. 因此需要對輸入水沙模型進行計算的流域面平均雨量進行修正，從而提高洪水和泥沙的模擬精度. 國內(nèi)外學者對降雨資料觀測誤差的分析與估計，目前研究較多的是對雷達測雨資料誤差進行的修正，然而以遙測系統(tǒng)為觀測手段，對流域尺度范圍的降雨觀測資料誤差進行分析和研究則開展得較少，在國外幾乎很少涉及[12-15]. 為提高模型模擬精度，司偉等[16-17]首先提出了動態(tài)系統(tǒng)響應曲線修正方法，該方法是向誤差源頭追溯的動態(tài)反饋修正方法，具有物理基礎強和不損失預見期等優(yōu)點[18]. 將該方法應用在新安江模型降雨誤差估計中，對流域面平均雨量進行修正以提高預報精度，取得了非常好的修正效果[10,19-21]. 基于以上研究，為了更精確模擬黃土區(qū)域水沙規(guī)律，本文嘗試將該方法應用在水沙模擬概念模型中，通過矯正降雨誤差，以提高水流和泥沙的模擬精度.

因此在水沙模擬概念模型計算中引入降雨動態(tài)系統(tǒng)響應曲線修正方法，對流域面平均雨量進行修正，建立面平均雨量和流域流量之間的響應關系. 降雨和產(chǎn)匯流過程是流域產(chǎn)沙匯沙過程的動力來源，水流是泥沙運移和輸送的載體，水流模擬精度與流域出口斷面的輸沙率過程息息相關. 鑒于此，本研究嘗試通過對面平均雨量進行修正，將修正后面平均雨量系列輸入模型進行計算重新得到新的流量和輸沙率過程，以達到提高水沙模型模擬精度的目的.

1 水沙模擬概念模型介紹

水沙模擬概念模型是包為民教授在20世紀末研究黃土地區(qū)流域泥沙提出的概念性模型，并被用于模擬蛇家溝和水旺溝流域等20多個小流域以及皇甫和河口鎮(zhèn)-龍門鎮(zhèn)區(qū)間等中大流域的水沙變化，取得了較好的模擬效果. 模型結(jié)構(gòu)簡單，物理概念清晰，適用于黃土區(qū)域干旱半干旱地區(qū)，以滿足水土保持措施效益分析和黃河中游水沙變化原因分析等問題研究的大尺度、長系列模擬的需要，為人類活動和氣候變化對流域水沙變化的影響提供基礎數(shù)據(jù)[22].

該模型根據(jù)水沙在流域上產(chǎn)生、運動各個環(huán)節(jié)機制的差異，分為水流模擬和泥沙模擬兩大部分[23]. 水流模擬分為產(chǎn)流、坡面匯流和溝道匯流3部分；泥沙模擬分為面蝕產(chǎn)沙、溝蝕產(chǎn)沙、坡面匯沙和溝道匯沙4部分. 產(chǎn)流機制采用垂向混合產(chǎn)流模式，坡面匯流采用線性水庫法，溝道匯流采用馬斯京根河道演算法；面蝕產(chǎn)沙采用土壤抗侵蝕能力分布曲線，溝蝕產(chǎn)沙采用拜格諾河道水流懸移質(zhì)泥沙公式，坡面匯沙和溝道匯沙采用泥沙匯集概念模擬公式.

模型輸入為降雨P和蒸發(fā)E資料，主要輸出為流域出口斷面的流量過程和泥沙過程. 模型主要結(jié)構(gòu)和參數(shù)如圖1和表1所示.

圖1 水沙模擬概念模型結(jié)構(gòu)圖Fig.1 Structural chart of conceptual model of water and sediment simulation

表1 水沙模擬概念模型參數(shù)

Tab.1 Parameters of conceptual model of water and sediment simulation

種類參數(shù)符號參數(shù)意義敏感程度取值范圍產(chǎn)流WM/mm流域平均土壤蓄水含量敏感180～350FC/(mm/min)流域平均穩(wěn)定下滲率敏感0.1～10BF流域下滲率分布曲線指數(shù)不敏感0.1～2KF土壤對下滲率影響系數(shù)不敏感0.1～2匯流B土壤蓄水容量分布曲線指數(shù)敏感0.1～1KG自由水箱地下水出流系數(shù)敏感0.1～1CS地面徑流線性水庫消退系數(shù)敏感0.1～0.65CI壤中流線性水庫消退系數(shù)敏感0.65～0.95CG地下徑流線性水庫消退系數(shù)敏感0.95～0.9999KE/min馬斯京根法河段傳播時間敏感0.5～2XE馬斯京根法流量比重系數(shù)敏感0.05～0.55產(chǎn)沙ACM/(kg/m3)坡面水力侵蝕能力系數(shù)敏感0.001～1REMM/kg坡面最大抗侵蝕能力不敏感>0BS抗侵蝕能力分布曲線指數(shù)不敏感0.1～5CM/(kg/m3)斷面最大含沙量不敏感5～2000CSS/(kg/m3)坡面毛溝侵蝕參數(shù)敏感0.1～250匯沙KXD非線性沖刷系數(shù)敏感0.1～0.999KS/min河段輸沙平均傳播時間敏感0.5～3XS河段輸沙權(quán)重系數(shù)敏感0.1～0.55ζ沖淤系數(shù)敏感0.01～1

2 降雨動態(tài)系統(tǒng)響應方法介紹

動態(tài)系統(tǒng)響應曲線修正方法是用于水文預報模型中的誤差修正技術(shù)，以模型輸入或某中間變量與輸出之間的對應關系建立響應曲線. 由于降雨是模型的重要輸入項，降雨誤差會在模型計算中進行傳遞[24]. 因此通過修正流域面平均雨量，提高模型的模擬精度.

本文將該方法用于水沙模擬概念模型中，將模型的水流模擬部分概化為一個系統(tǒng)(圖2)，建立動態(tài)系統(tǒng)響應曲線函數(shù)[17]:

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

(1)

圖2 水沙模型概化系統(tǒng)Fig.2 Generalization system of conceptual model of water and sediment simulation

在本文中，只考慮降雨對水沙模型流量的響應過程，因此降雨動態(tài)系統(tǒng)響應曲線函數(shù)簡化為：

Q(P)=f(P)

(2)

所謂降雨動態(tài)系統(tǒng)響應曲線，即對某時段的降雨Pi增加一個單位值，將新的降雨系列輸入模型進行計算得到的流量與原始計算流量之間的差值.

把降雨輸入作為自變量對流量求全微分，并將觀測流量系列代入式中，則其矩陣表達式為：

Q(P)≈QC(PO)+UΔP+W

(3)

式中，ΔP是需求解的降雨系列誤差(mm)；W是經(jīng)過實時修正后的預報系統(tǒng)殘差(m3/s)；U為系統(tǒng)響應矩陣.

U矩陣如下：

(4)

計算降雨動態(tài)系統(tǒng)響應曲線的步驟如下：

(1)假設流域有降雨系列PO，實測流量為QO，PO輸入模型計算得到流域原始計算流量系列為QC.

(2)流域各站點降雨系列PO中第i個時段Pi，在其余時段各站點降雨均不變的情況下，增加一個單位值，得到新的降雨系列PCi.

(3)將新的降雨系列PCi輸入模型進行計算，得到新的計算流量過程QCi.

(4)將新的流量過程QCi與原始計算流量QC相減，其差值為Pi的系統(tǒng)響應曲線，屬于U矩陣中的第i列.

(5)依次計算U矩陣中的每一列，最終建立降雨系列PO的動態(tài)系統(tǒng)響應曲線. 運用最小二乘法，計算得到降雨的修正誤差值：

ΔP=(UTU)-1UT(QO-QC)

(5)

(6)修正后的降雨系列:

PU=PO+ΔP

(6)

輸入模型進行計算，新的流域出口斷面流量和輸沙率過程即該流域流量和輸沙率的修正值.

3 應用檢驗

檢驗方法修正效果，具體采用以下4個指標來衡量修正效果：

(1)徑流深相對誤差ΔR，本文取誤差絕對值：

(7)

(2)產(chǎn)沙量相對誤差ΔS，本文取誤差絕對值：

(8)

(3)納什效率系數(shù)NS：

(9)

(4)納什效率系數(shù)提高幅度INS：

(10)

3.1 理想案例檢驗

在理想案例中，假設模型的輸入、輸出、參數(shù)和初始狀態(tài)值已知. 假設一個降雨系列PO，將其輸入水沙模型計算得到的流量和輸沙率QO和SO，作為理想案例的實測系列. 隨機給定一組服從零均值分布的降雨誤差值ΔP，誤差值不超過降雨系列的±80%，將加上誤差值的新的降雨系列PC輸入模型計算得到計算流量QC和計算輸沙率SC. 其中，

PC=PO+ΔP

(11)

對PC系列降雨進行動態(tài)系統(tǒng)響應修正，修正得到降雨誤差值為-ΔP′，修正的降雨系列為:

PU=PC-ΔP′=PO+ΔP-ΔP′

(12)

PU輸入模型計算得到修正后流量QU和輸沙率SU. 分別將修正降雨誤差值-ΔP′與給定降雨誤差值ΔP進行對比，比較兩者相關性，檢驗該方法是否能夠精確反演降雨誤差；分別將修正后流量、輸沙率和計算流量、輸沙率與實測流量、輸沙率進行對比分析，檢驗修正后的模型模擬精度是否得到提高.

在本研究中，理想案例降雨系列采用曹坪流域1962年7月23日場次洪水降雨資料，流域面積187 km2，有30個雨量站，洪水匯流演算馬斯京根演算法河段數(shù)為10. 給定隨機降雨誤差值建立理想案例，表2為理想案例的徑流深和產(chǎn)沙量修正情況，圖3為理想案例給定降雨誤差和修正降雨誤差值對比圖，圖4為理想案例降雨數(shù)據(jù)的修正效果，圖5為理想案例水流和泥沙修正效果.

表2 理想案例流量和輸沙率修正效果

1)RO為實測徑流深；2)RC為模型原始預報徑流深；3)ΔR1為模型原始預報徑流深與實測徑流深之間相對誤差；4)RU為模型修正后徑流深；5)ΔR2為模型修正徑流深與實測徑流深之間相對誤差；6)SO為實測產(chǎn)沙量；7)SC為模型原始計算產(chǎn)沙量；8)ΔS1為模型原始計算產(chǎn)沙量與實測產(chǎn)沙量之間相對誤差；9)SU為模型修正后產(chǎn)沙量；10)ΔS2為模型修正產(chǎn)沙量與實測產(chǎn)沙量之間相對誤差.

從表2可以看出，應用降雨動態(tài)系統(tǒng)響應曲線后，理想案例水流部分徑流深誤差由8.29%降低到2.67%，納什效率系數(shù)由0.989提高到0.997；泥沙部分產(chǎn)沙量相對誤差由6.96%降低到4.77%，納什效率系數(shù)由0.990提高到0.996. 從圖3可以看出，給定降雨誤差與修正得到的降雨誤差之間NS系數(shù)為0.821，相關系數(shù)為0.908，二者非常接近. 從圖4可以看出，相比較有誤差的降雨系列PC，修正后的降雨系列PU更接近于降雨基準數(shù)據(jù)PO系列，這說明本文中所采用的降雨誤差估計方法能夠準確地估計出給定的降雨誤差，精確反演降雨誤差. 從圖5可以看出，修正后的流量和輸沙率更接近流量真值和輸沙率真值，流量和輸沙率擬合過程都有所改善. 這說明降雨動態(tài)系統(tǒng)響應曲線修正方法能夠同時提高水沙模型的水流和泥沙模擬精度，可以進一步用于實際流域的應用檢驗.

圖3 理想案例給定雨量誤差ΔP及修正雨量誤差ΔP′值(圖中1個時刻為模型計算時間尺度0.5 h)Fig.3 Given rainfall error ΔP and corrected rainfall error ΔP′ of ideal case

圖4 理想案例修正前后降雨值(圖中1個時刻為模型計算時間尺度0.5 h)Fig.4 Rainfall values before and after revision of ideal case

圖5 理想案例修正效果(圖中1個時刻為模型計算時間尺度0.5 h)Fig.5 Correction effect of ideal case

3.2 實際流域應用

本文選取黃土地區(qū)曹坪流域作為實驗流域. 曹坪流域位于黃河大理河流域中下游地區(qū)，流域大部分為黃土丘陵溝壑區(qū)域，植被稀疏，水土流失嚴重. 是典型的大陸性季風氣候區(qū)，夏季炎熱潮濕多暴雨，冬季寒冷干燥. 降雨年內(nèi)變化大，降雨多集中在7-9月. 流域多年平均降雨量447 mm，多年平均氣溫7.8～9.6℃，流域泥沙侵蝕模數(shù)為2.2萬t/(km2·a).

研究流域曹坪(岔巴溝)為干旱半干旱氣候，產(chǎn)流模式符合混合產(chǎn)流模式，可以運用水沙概念模擬模型進行計算. 用于計算修正的資料系列為1961-1989年的13場洪水資料，模型計算時間尺度為0.5 h，流域水系及站網(wǎng)布設見圖6.

圖6 曹坪流域水系及站網(wǎng)布設Fig.6 Layout of water system and station network in Caoping Basin

水沙模型經(jīng)函數(shù)曲面參數(shù)率定方法[25]進行參數(shù)率定后，對降雨進行動態(tài)系統(tǒng)響應修正，檢驗模型的修正效果. 模型參數(shù)率定情況見表3. 13場洪水水流和泥沙修正結(jié)果分別見表4和表5.

表3 水沙模型曹坪流域參數(shù)率定

表4 曹坪流域水流修正效果

表5 曹坪流域泥沙修正效果

圖7 曹坪流域13場洪水納什效率系數(shù)修正前后對比Fig.7 Contrast of NS coefficient before and after correction of 13 floods in Caoping Basin

圖8 第8場次洪水修正效果(圖中1個時刻為模型計算時間尺度0.5 h)Fig.8 Correction effect of flood 8

從表4洪水修正效果來看，應用降雨動態(tài)系統(tǒng)響應曲線后，修正后的徑流深相對誤差與修正前相比都有所減小，納什效率系數(shù)均有所提高. 平均徑流深相對誤差由2.48%降低到0.62%，納什效率系數(shù)由0.689提高到0.810，納什效率系數(shù)提高幅度為17.56%. 從表5泥沙修正效果來看，產(chǎn)沙量經(jīng)修正后相對誤差全部降低，納什效率系數(shù)全部提高. 平均產(chǎn)沙量相對誤差從1.71%降低到1.28%，納什效率系數(shù)從0.643提高到0.745，納什效率系數(shù)提高幅度為15.86%. 從圖7洪水的納什效率系數(shù)修正效果來看，修正后的納什效率系數(shù)均比修正前有顯著提升. 選取第8場次洪水的流量和輸沙率修正結(jié)果查看修正效果，從圖8可以看出，洪水修正后的流量和輸沙率更接近流量真值和輸沙率真值，流量和輸沙率擬合過程也有所改善. 這說明降雨動態(tài)系統(tǒng)響應曲線修正方法應用于水沙模型后，能夠減小模型徑流深和產(chǎn)沙量的相對誤差，提高納什效率系數(shù)，最終提高模型的模擬精度.

在現(xiàn)今應用的黃河流域土壤侵蝕產(chǎn)沙模型中，流域的復雜性和資料的制約，造成模型模擬結(jié)果精度欠佳的不足[26]. 流域面平均雨量的誤差，影響了模型的模擬精度，使模型初始模擬效果沒有想象中理想. 用于降雨動態(tài)系統(tǒng)響應修正的13場洪水降雨大多分為兩種情況：一是由于原始降雨資料存在降雨時段均化、缺測、漏測等問題，部分降雨數(shù)據(jù)由插補展延而來導致降雨誤差. 這種類型降雨輸入模型進行計算后，相比實際流量和輸沙率過程，計算洪水流量和輸沙率過程相對呈現(xiàn)“矮胖型”，洪水緩漲緩落，洪峰、沙峰偏低，模型模擬精度不高. 二是由于原始資料存在誤測等問題，某些時刻降雨數(shù)據(jù)明顯偏大或者偏小，造成降雨數(shù)據(jù)誤差. 這部分降雨數(shù)據(jù)輸入模型后也會造成計算誤差，使計算流量和輸沙率過程與實際過程不相符合，模型模擬精度較低. 降雨是模型計算的重要輸入項，也是模型誤差的重要來源之一. 因此，在使用水沙模型模擬黃土區(qū)域水沙時，對于有誤差的降雨資料，尤其針對以上兩種降雨分布情況的洪水場次，可以應用降雨動態(tài)系統(tǒng)響應曲線修正方法，通過修正面平均雨量，以提高模型模擬精度.

4 結(jié)論

水沙模擬概念模型物理概念清楚，實用性強，適用于黃土區(qū)域干旱半干旱地區(qū). 應用于黃土地區(qū)流域，合理地模擬了流域水沙變化規(guī)律. 降雨是水沙模型的重要輸入項，降雨誤差能夠直接造成模型誤差，從而影響模型的模擬精度. 因此，本文將降雨動態(tài)系統(tǒng)響應曲線修正方法與水沙模型相結(jié)合，通過修正面平均雨量，將修正后降雨輸入模型進行計算重新得到新的流量和輸沙率過程，以達到減小模型誤差的目的. 結(jié)果表明，在對降雨進行修正后，水沙模型水流徑流深更接近實際值，流量過程納什效率系數(shù)提高；產(chǎn)沙量誤差降低，輸沙率過程納什效率系數(shù)提高. 說明水沙模型應用降雨動態(tài)系統(tǒng)響應曲線修正方法后，能夠提高模型的模擬精度，最終水流和泥沙模擬平均提高幅度分別是17.56%和15.86%，其修正效果是顯著的.

降雨動態(tài)系統(tǒng)響應曲線修正方法能夠很好地提高水沙模型的模擬精度，但是由于水文系統(tǒng)畢竟是非線性系統(tǒng)，影響因素比較復雜，降雨誤差因素并非模型誤差的唯一來源. 考慮到本文僅修正面平均雨量一個變量，因此，對于如何細化降雨誤差，以及綜合考慮降雨和其他影響因素之間量化關系有待進一步研究. 根據(jù)流域降雨分布情況，選擇降雨動態(tài)系統(tǒng)響應曲線修正方法應用于水沙模型模擬中，可以為研究黃河水沙變化規(guī)律提供一定的技術(shù)指導和理論支撐.