包為民,顧雨薇,司 偉**,侯 露,盧金利,落全富
(1:河海大學(xué)水文水資源學(xué)院,南京 210098)(2:臺(tái)州市水文站,臺(tái)州 318001)(3:青山水庫(kù)管理處,杭州 311305)
黃土地區(qū)泥沙問題是目前水文研究的重點(diǎn)問題之一. 研究黃河水沙變化、把握黃河水沙變化規(guī)律,對(duì)于進(jìn)一步完善治黃方略、實(shí)施流域水資源配置與管理及重大水利工程布局,意義重大[1]. 目前泥沙模型的研究取得了較為豐碩的成果,主要包括經(jīng)驗(yàn)相關(guān)模型和物理成因模型. 經(jīng)驗(yàn)相關(guān)模型中,Wischmeier WH等[2]提出了著名的通用土壤侵蝕方程USLE. 基于物理機(jī)理的流域水沙模型研究較多,湯立群等[3]基于坡面徑流侵蝕量計(jì)算公式提出了流域產(chǎn)沙模型;WEPP、LISEM、EUROSEM等模型[4-6]基于土壤侵蝕過程的物理描述,模擬了流域侵蝕過程,在許多國(guó)家被成功推廣應(yīng)用. 已有的大多數(shù)流域產(chǎn)沙模型未能考慮不同流域尺度的泥沙規(guī)律,模型模擬精度以及使用范圍受到了很大限制. 因此,包為民[7]基于流域自然地理特性,在對(duì)黃土地區(qū)產(chǎn)沙機(jī)制進(jìn)行概化的基礎(chǔ)上,提出了結(jié)構(gòu)和參數(shù)均具有物理意義的流域水沙模擬概念模型. 該模型經(jīng)黃土地區(qū)不同流域尺度的實(shí)際流域和前期大量的試驗(yàn)流域驗(yàn)證,可以很好地模擬黃土高原地區(qū)流域的水沙規(guī)律[8].
水沙模擬中重要的動(dòng)力條件是水流,因此水流模擬的精度直接影響流域產(chǎn)沙和匯沙過程,進(jìn)而影響水沙模型的模擬精度和效果[9]. 水流模擬的主要技術(shù)手段是依賴于水文模型. 水文模型計(jì)算誤差來源有很多,其中包括重要輸入項(xiàng)降雨資料、下墊面條件、流域初始狀態(tài)值、模型本身誤差等[10]. 由于降雨是水文預(yù)報(bào)模型計(jì)算的重要輸入項(xiàng),因此降雨資料直接影響模型產(chǎn)流,進(jìn)而影響模型的模擬精度[11]. 在水沙模型計(jì)算過程中,基于黃土區(qū)域自然地理特性,流域面平均雨量大多采用“以點(diǎn)帶面”的方式進(jìn)行計(jì)算,即利用雨量站點(diǎn)的觀測(cè)資料代表該雨量站所代表的單元面上的降雨,這種方式不可避免的給流域降雨帶來誤差. 同時(shí),部分降雨資料由于觀測(cè)等因素存在時(shí)段均化現(xiàn)象,伴有一定程度的缺測(cè)、漏測(cè)和誤測(cè)問題,降雨資料由插補(bǔ)展延而來, 造成模型模擬誤差. 因此需要對(duì)輸入水沙模型進(jìn)行計(jì)算的流域面平均雨量進(jìn)行修正,從而提高洪水和泥沙的模擬精度. 國(guó)內(nèi)外學(xué)者對(duì)降雨資料觀測(cè)誤差的分析與估計(jì),目前研究較多的是對(duì)雷達(dá)測(cè)雨資料誤差進(jìn)行的修正,然而以遙測(cè)系統(tǒng)為觀測(cè)手段,對(duì)流域尺度范圍的降雨觀測(cè)資料誤差進(jìn)行分析和研究則開展得較少,在國(guó)外幾乎很少涉及[12-15]. 為提高模型模擬精度,司偉等[16-17]首先提出了動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法,該方法是向誤差源頭追溯的動(dòng)態(tài)反饋修正方法,具有物理基礎(chǔ)強(qiáng)和不損失預(yù)見期等優(yōu)點(diǎn)[18]. 將該方法應(yīng)用在新安江模型降雨誤差估計(jì)中,對(duì)流域面平均雨量進(jìn)行修正以提高預(yù)報(bào)精度,取得了非常好的修正效果[10,19-21]. 基于以上研究,為了更精確模擬黃土區(qū)域水沙規(guī)律,本文嘗試將該方法應(yīng)用在水沙模擬概念模型中,通過矯正降雨誤差,以提高水流和泥沙的模擬精度.
因此在水沙模擬概念模型計(jì)算中引入降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法,對(duì)流域面平均雨量進(jìn)行修正,建立面平均雨量和流域流量之間的響應(yīng)關(guān)系. 降雨和產(chǎn)匯流過程是流域產(chǎn)沙匯沙過程的動(dòng)力來源,水流是泥沙運(yùn)移和輸送的載體,水流模擬精度與流域出口斷面的輸沙率過程息息相關(guān). 鑒于此,本研究嘗試通過對(duì)面平均雨量進(jìn)行修正,將修正后面平均雨量系列輸入模型進(jìn)行計(jì)算重新得到新的流量和輸沙率過程,以達(dá)到提高水沙模型模擬精度的目的.
水沙模擬概念模型是包為民教授在20世紀(jì)末研究黃土地區(qū)流域泥沙提出的概念性模型,并被用于模擬蛇家溝和水旺溝流域等20多個(gè)小流域以及皇甫和河口鎮(zhèn)-龍門鎮(zhèn)區(qū)間等中大流域的水沙變化,取得了較好的模擬效果. 模型結(jié)構(gòu)簡(jiǎn)單,物理概念清晰,適用于黃土區(qū)域干旱半干旱地區(qū),以滿足水土保持措施效益分析和黃河中游水沙變化原因分析等問題研究的大尺度、長(zhǎng)系列模擬的需要,為人類活動(dòng)和氣候變化對(duì)流域水沙變化的影響提供基礎(chǔ)數(shù)據(jù)[22].
該模型根據(jù)水沙在流域上產(chǎn)生、運(yùn)動(dòng)各個(gè)環(huán)節(jié)機(jī)制的差異,分為水流模擬和泥沙模擬兩大部分[23]. 水流模擬分為產(chǎn)流、坡面匯流和溝道匯流3部分;泥沙模擬分為面蝕產(chǎn)沙、溝蝕產(chǎn)沙、坡面匯沙和溝道匯沙4部分. 產(chǎn)流機(jī)制采用垂向混合產(chǎn)流模式,坡面匯流采用線性水庫(kù)法,溝道匯流采用馬斯京根河道演算法;面蝕產(chǎn)沙采用土壤抗侵蝕能力分布曲線,溝蝕產(chǎn)沙采用拜格諾河道水流懸移質(zhì)泥沙公式,坡面匯沙和溝道匯沙采用泥沙匯集概念模擬公式.
模型輸入為降雨P(guān)和蒸發(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ù)符號(hào)參數(shù)意義敏感程度取值范圍產(chǎn)流WM/mm流域平均土壤蓄水含量敏感180~350FC/(mm/min)流域平均穩(wěn)定下滲率敏感0.1~10BF流域下滲率分布曲線指數(shù)不敏感0.1~2KF土壤對(duì)下滲率影響系數(shù)不敏感0.1~2匯流B土壤蓄水容量分布曲線指數(shù)敏感0.1~1KG自由水箱地下水出流系數(shù)敏感0.1~1CS地面徑流線性水庫(kù)消退系數(shù)敏感0.1~0.65CI壤中流線性水庫(kù)消退系數(shù)敏感0.65~0.95CG地下徑流線性水庫(kù)消退系數(shù)敏感0.95~0.9999KE/min馬斯京根法河段傳播時(shí)間敏感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河段輸沙平均傳播時(shí)間敏感0.5~3XS河段輸沙權(quán)重系數(shù)敏感0.1~0.55ζ沖淤系數(shù)敏感0.01~1
動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法是用于水文預(yù)報(bào)模型中的誤差修正技術(shù),以模型輸入或某中間變量與輸出之間的對(duì)應(yīng)關(guān)系建立響應(yīng)曲線. 由于降雨是模型的重要輸入項(xiàng),降雨誤差會(huì)在模型計(jì)算中進(jìn)行傳遞[24]. 因此通過修正流域面平均雨量,提高模型的模擬精度.
本文將該方法用于水沙模擬概念模型中,將模型的水流模擬部分概化為一個(gè)系統(tǒng)(圖2),建立動(dòng)態(tài)系統(tǒng)響應(yī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
在本文中,只考慮降雨對(duì)水沙模型流量的響應(yīng)過程,因此降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線函數(shù)簡(jiǎn)化為:
Q(P)=f(P)
(2)
所謂降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線,即對(duì)某時(shí)段的降雨P(guān)i增加一個(gè)單位值,將新的降雨系列輸入模型進(jìn)行計(jì)算得到的流量與原始計(jì)算流量之間的差值.
把降雨輸入作為自變量對(duì)流量求全微分,并將觀測(cè)流量系列代入式中,則其矩陣表達(dá)式為:
Q(P)≈QC(PO)+UΔP+W
(3)
式中,ΔP是需求解的降雨系列誤差(mm);W是經(jīng)過實(shí)時(shí)修正后的預(yù)報(bào)系統(tǒng)殘差(m3/s);U為系統(tǒng)響應(yīng)矩陣.
U矩陣如下:
(4)
計(jì)算降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線的步驟如下:
(1)假設(shè)流域有降雨系列PO,實(shí)測(cè)流量為QO,PO輸入模型計(jì)算得到流域原始計(jì)算流量系列為QC.
(2)流域各站點(diǎn)降雨系列PO中第i個(gè)時(shí)段Pi,在其余時(shí)段各站點(diǎn)降雨均不變的情況下,增加一個(gè)單位值,得到新的降雨系列PCi.
(3)將新的降雨系列PCi輸入模型進(jìn)行計(jì)算,得到新的計(jì)算流量過程QCi.
(4)將新的流量過程QCi與原始計(jì)算流量QC相減,其差值為Pi的系統(tǒng)響應(yīng)曲線,屬于U矩陣中的第i列.
(5)依次計(jì)算U矩陣中的每一列,最終建立降雨系列PO的動(dòng)態(tài)系統(tǒng)響應(yīng)曲線. 運(yùn)用最小二乘法,計(jì)算得到降雨的修正誤差值:
ΔP=(UTU)-1UT(QO-QC)
(5)
(6)修正后的降雨系列:
PU=PO+ΔP
(6)
輸入模型進(jìn)行計(jì)算,新的流域出口斷面流量和輸沙率過程即該流域流量和輸沙率的修正值.
檢驗(yàn)方法修正效果,具體采用以下4個(gè)指標(biāo)來衡量修正效果:
(1)徑流深相對(duì)誤差ΔR,本文取誤差絕對(duì)值:
(7)
(2)產(chǎn)沙量相對(duì)誤差ΔS,本文取誤差絕對(duì)值:
(8)
(3)納什效率系數(shù)NS:
(9)
(4)納什效率系數(shù)提高幅度INS:
(10)
在理想案例中,假設(shè)模型的輸入、輸出、參數(shù)和初始狀態(tài)值已知. 假設(shè)一個(gè)降雨系列PO,將其輸入水沙模型計(jì)算得到的流量和輸沙率QO和SO,作為理想案例的實(shí)測(cè)系列. 隨機(jī)給定一組服從零均值分布的降雨誤差值ΔP,誤差值不超過降雨系列的±80%,將加上誤差值的新的降雨系列PC輸入模型計(jì)算得到計(jì)算流量QC和計(jì)算輸沙率SC. 其中,
PC=PO+ΔP
(11)
對(duì)PC系列降雨進(jìn)行動(dòng)態(tài)系統(tǒng)響應(yīng)修正,修正得到降雨誤差值為-ΔP′,修正的降雨系列為:
PU=PC-ΔP′=PO+ΔP-ΔP′
(12)
PU輸入模型計(jì)算得到修正后流量QU和輸沙率SU. 分別將修正降雨誤差值-ΔP′與給定降雨誤差值ΔP進(jìn)行對(duì)比,比較兩者相關(guān)性,檢驗(yàn)該方法是否能夠精確反演降雨誤差;分別將修正后流量、輸沙率和計(jì)算流量、輸沙率與實(shí)測(cè)流量、輸沙率進(jìn)行對(duì)比分析,檢驗(yàn)修正后的模型模擬精度是否得到提高.
在本研究中,理想案例降雨系列采用曹坪流域1962年7月23日?qǐng)龃魏樗涤曩Y料,流域面積187 km2,有30個(gè)雨量站,洪水匯流演算馬斯京根演算法河段數(shù)為10. 給定隨機(jī)降雨誤差值建立理想案例,表2為理想案例的徑流深和產(chǎn)沙量修正情況,圖3為理想案例給定降雨誤差和修正降雨誤差值對(duì)比圖,圖4為理想案例降雨數(shù)據(jù)的修正效果,圖5為理想案例水流和泥沙修正效果.
表2 理想案例流量和輸沙率修正效果
1)RO為實(shí)測(cè)徑流深;2)RC為模型原始預(yù)報(bào)徑流深;3)ΔR1為模型原始預(yù)報(bào)徑流深與實(shí)測(cè)徑流深之間相對(duì)誤差;4)RU為模型修正后徑流深;5)ΔR2為模型修正徑流深與實(shí)測(cè)徑流深之間相對(duì)誤差;6)SO為實(shí)測(cè)產(chǎn)沙量;7)SC為模型原始計(jì)算產(chǎn)沙量;8)ΔS1為模型原始計(jì)算產(chǎn)沙量與實(shí)測(cè)產(chǎn)沙量之間相對(duì)誤差;9)SU為模型修正后產(chǎn)沙量;10)ΔS2為模型修正產(chǎn)沙量與實(shí)測(cè)產(chǎn)沙量之間相對(duì)誤差.
從表2可以看出,應(yīng)用降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線后,理想案例水流部分徑流深誤差由8.29%降低到2.67%,納什效率系數(shù)由0.989提高到0.997;泥沙部分產(chǎn)沙量相對(duì)誤差由6.96%降低到4.77%,納什效率系數(shù)由0.990提高到0.996. 從圖3可以看出,給定降雨誤差與修正得到的降雨誤差之間NS系數(shù)為0.821,相關(guān)系數(shù)為0.908,二者非常接近. 從圖4可以看出,相比較有誤差的降雨系列PC,修正后的降雨系列PU更接近于降雨基準(zhǔn)數(shù)據(jù)PO系列,這說明本文中所采用的降雨誤差估計(jì)方法能夠準(zhǔn)確地估計(jì)出給定的降雨誤差,精確反演降雨誤差. 從圖5可以看出,修正后的流量和輸沙率更接近流量真值和輸沙率真值,流量和輸沙率擬合過程都有所改善. 這說明降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法能夠同時(shí)提高水沙模型的水流和泥沙模擬精度,可以進(jìn)一步用于實(shí)際流域的應(yīng)用檢驗(yàn).
圖3 理想案例給定雨量誤差ΔP及修正雨量誤差ΔP′值(圖中1個(gè)時(shí)刻為模型計(jì)算時(shí)間尺度0.5 h)Fig.3 Given rainfall error ΔP and corrected rainfall error ΔP′ of ideal case
圖4 理想案例修正前后降雨值(圖中1個(gè)時(shí)刻為模型計(jì)算時(shí)間尺度0.5 h)Fig.4 Rainfall values before and after revision of ideal case
圖5 理想案例修正效果(圖中1個(gè)時(shí)刻為模型計(jì)算時(shí)間尺度0.5 h)Fig.5 Correction effect of ideal case
本文選取黃土地區(qū)曹坪流域作為實(shí)驗(yàn)流域. 曹坪流域位于黃河大理河流域中下游地區(qū),流域大部分為黃土丘陵溝壑區(qū)域,植被稀疏,水土流失嚴(yán)重. 是典型的大陸性季風(fēng)氣候區(qū),夏季炎熱潮濕多暴雨,冬季寒冷干燥. 降雨年內(nèi)變化大,降雨多集中在7-9月. 流域多年平均降雨量447 mm,多年平均氣溫7.8~9.6℃,流域泥沙侵蝕模數(shù)為2.2萬t/(km2·a).
研究流域曹坪(岔巴溝)為干旱半干旱氣候,產(chǎn)流模式符合混合產(chǎn)流模式,可以運(yùn)用水沙概念模擬模型進(jìn)行計(jì)算. 用于計(jì)算修正的資料系列為1961-1989年的13場(chǎng)洪水資料,模型計(jì)算時(shí)間尺度為0.5 h,流域水系及站網(wǎng)布設(shè)見圖6.
圖6 曹坪流域水系及站網(wǎng)布設(shè)Fig.6 Layout of water system and station network in Caoping Basin
水沙模型經(jīng)函數(shù)曲面參數(shù)率定方法[25]進(jìn)行參數(shù)率定后,對(duì)降雨進(jìn)行動(dòng)態(tài)系統(tǒng)響應(yīng)修正,檢驗(yàn)?zāi)P偷男拚Ч? 模型參數(shù)率定情況見表3. 13場(chǎng)洪水水流和泥沙修正結(jié)果分別見表4和表5.
表3 水沙模型曹坪流域參數(shù)率定
表4 曹坪流域水流修正效果
表5 曹坪流域泥沙修正效果
圖7 曹坪流域13場(chǎng)洪水納什效率系數(shù)修正前后對(duì)比Fig.7 Contrast of NS coefficient before and after correction of 13 floods in Caoping Basin
圖8 第8場(chǎng)次洪水修正效果(圖中1個(gè)時(shí)刻為模型計(jì)算時(shí)間尺度0.5 h)Fig.8 Correction effect of flood 8
從表4洪水修正效果來看,應(yīng)用降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線后,修正后的徑流深相對(duì)誤差與修正前相比都有所減小,納什效率系數(shù)均有所提高. 平均徑流深相對(duì)誤差由2.48%降低到0.62%,納什效率系數(shù)由0.689提高到0.810,納什效率系數(shù)提高幅度為17.56%. 從表5泥沙修正效果來看,產(chǎn)沙量經(jīng)修正后相對(duì)誤差全部降低,納什效率系數(shù)全部提高. 平均產(chǎn)沙量相對(duì)誤差從1.71%降低到1.28%,納什效率系數(shù)從0.643提高到0.745,納什效率系數(shù)提高幅度為15.86%. 從圖7洪水的納什效率系數(shù)修正效果來看,修正后的納什效率系數(shù)均比修正前有顯著提升. 選取第8場(chǎng)次洪水的流量和輸沙率修正結(jié)果查看修正效果,從圖8可以看出,洪水修正后的流量和輸沙率更接近流量真值和輸沙率真值,流量和輸沙率擬合過程也有所改善. 這說明降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法應(yīng)用于水沙模型后,能夠減小模型徑流深和產(chǎn)沙量的相對(duì)誤差,提高納什效率系數(shù),最終提高模型的模擬精度.
在現(xiàn)今應(yīng)用的黃河流域土壤侵蝕產(chǎn)沙模型中,流域的復(fù)雜性和資料的制約,造成模型模擬結(jié)果精度欠佳的不足[26]. 流域面平均雨量的誤差,影響了模型的模擬精度,使模型初始模擬效果沒有想象中理想. 用于降雨動(dòng)態(tài)系統(tǒng)響應(yīng)修正的13場(chǎng)洪水降雨大多分為兩種情況:一是由于原始降雨資料存在降雨時(shí)段均化、缺測(cè)、漏測(cè)等問題,部分降雨數(shù)據(jù)由插補(bǔ)展延而來導(dǎo)致降雨誤差. 這種類型降雨輸入模型進(jìn)行計(jì)算后,相比實(shí)際流量和輸沙率過程,計(jì)算洪水流量和輸沙率過程相對(duì)呈現(xiàn)“矮胖型”,洪水緩漲緩落,洪峰、沙峰偏低,模型模擬精度不高. 二是由于原始資料存在誤測(cè)等問題,某些時(shí)刻降雨數(shù)據(jù)明顯偏大或者偏小,造成降雨數(shù)據(jù)誤差. 這部分降雨數(shù)據(jù)輸入模型后也會(huì)造成計(jì)算誤差,使計(jì)算流量和輸沙率過程與實(shí)際過程不相符合,模型模擬精度較低. 降雨是模型計(jì)算的重要輸入項(xiàng),也是模型誤差的重要來源之一. 因此,在使用水沙模型模擬黃土區(qū)域水沙時(shí),對(duì)于有誤差的降雨資料,尤其針對(duì)以上兩種降雨分布情況的洪水場(chǎng)次,可以應(yīng)用降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法,通過修正面平均雨量,以提高模型模擬精度.
水沙模擬概念模型物理概念清楚,實(shí)用性強(qiáng),適用于黃土區(qū)域干旱半干旱地區(qū). 應(yīng)用于黃土地區(qū)流域,合理地模擬了流域水沙變化規(guī)律. 降雨是水沙模型的重要輸入項(xiàng),降雨誤差能夠直接造成模型誤差,從而影響模型的模擬精度. 因此,本文將降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法與水沙模型相結(jié)合,通過修正面平均雨量,將修正后降雨輸入模型進(jìn)行計(jì)算重新得到新的流量和輸沙率過程,以達(dá)到減小模型誤差的目的. 結(jié)果表明,在對(duì)降雨進(jìn)行修正后,水沙模型水流徑流深更接近實(shí)際值,流量過程納什效率系數(shù)提高;產(chǎn)沙量誤差降低,輸沙率過程納什效率系數(shù)提高. 說明水沙模型應(yīng)用降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法后,能夠提高模型的模擬精度,最終水流和泥沙模擬平均提高幅度分別是17.56%和15.86%,其修正效果是顯著的.
降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法能夠很好地提高水沙模型的模擬精度,但是由于水文系統(tǒng)畢竟是非線性系統(tǒng),影響因素比較復(fù)雜,降雨誤差因素并非模型誤差的唯一來源. 考慮到本文僅修正面平均雨量一個(gè)變量,因此,對(duì)于如何細(xì)化降雨誤差,以及綜合考慮降雨和其他影響因素之間量化關(guān)系有待進(jìn)一步研究. 根據(jù)流域降雨分布情況,選擇降雨動(dòng)態(tài)系統(tǒng)響應(yīng)曲線修正方法應(yīng)用于水沙模型模擬中,可以為研究黃河水沙變化規(guī)律提供一定的技術(shù)指導(dǎo)和理論支撐.