吳夢(mèng)琪，廖衛(wèi)紅，廖春梅，劉 攀

(1. 武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072；2.中國水利水電科學(xué)院，北京 100038；3.國家電網(wǎng)公司西北分部，西安 710048)

0 引 言

增江是東江的一級(jí)支流，縱貫增城市中心城區(qū)，增江流域地區(qū)為自北向南地勢(shì)由高逐漸變低，中南部丘陵山地較多，然而由于降雨時(shí)空分布不均，造成流域內(nèi)徑流量的時(shí)空分布不均。同時(shí)，隨著經(jīng)濟(jì)的不斷發(fā)展和人民生活水平的不斷提高，地區(qū)用水量急劇增加、水質(zhì)日益惡化，致使廣州市增城區(qū)出現(xiàn)水質(zhì)性缺水[1]。

對(duì)增江流域的徑流量變化、流域內(nèi)降水的時(shí)空分布等方面研究較多，但基于水文模型對(duì)增江流域水文循環(huán)過程進(jìn)行模擬的研究還相對(duì)較少，其中毛維新等[2]研究了數(shù)字高程模型和新安江模型相結(jié)合的半分布式水文模型的建模思路，并證明在該流域內(nèi)切實(shí)可行；劉戰(zhàn)友等[3]針對(duì)增江流域水循環(huán)和水文資料特點(diǎn)，建立了三水源新安江模型與GIS結(jié)合的半分布式水文模型。本文考慮到增江流域地區(qū)為丘陵地帶，下墊面情況復(fù)雜，因此建立增江流域分布式水文模型，通過對(duì)水文站的徑流進(jìn)行模擬，為下一步實(shí)現(xiàn)流域洪水預(yù)報(bào)及水資源評(píng)價(jià)工作奠定基礎(chǔ)。

1 分布式水文模型EasyDHM

本文采用的是中國水利水電科學(xué)研究院自主研發(fā)的分布式水文模型EasyDHM模型。目前，分布式水文模型的難點(diǎn)在于模型的建立，而EasyDHM模型擁有一套自己的建模工具，由EasyDHM模型代碼和MWEasyDHM系統(tǒng)組成。EasyDHM模型代碼包含產(chǎn)流模塊、匯流模塊、地下水模塊、水質(zhì)模塊、土壤侵蝕模塊、人工用水模塊、水庫調(diào)度模塊、參數(shù)識(shí)別模塊等8大模塊[4]。EasyDHM模型系統(tǒng)--MWEasyDHM以開源GIS軟件MapWindow作為基礎(chǔ)，由以下兩大模塊組成：模型前處理，其中包括水文分析、子流域劃分、模型參數(shù)推求、氣象數(shù)據(jù)插值等功能，為EasyDHM模型提供了一個(gè)快速建模的軟件環(huán)境；模型后處理，其中包括模型實(shí)時(shí)計(jì)算、參數(shù)自動(dòng)/手動(dòng)識(shí)別模塊、模型結(jié)果分析、模型結(jié)果時(shí)空展示等，為EasyDHM模型提供了一個(gè)快速計(jì)算、快速率定的軟件環(huán)境[5]。

該模型基于流域DEM、河流分布、土壤分布、土地利用分布等下墊面數(shù)據(jù)對(duì)研究區(qū)域進(jìn)行單元?jiǎng)澐?。首先，根?jù)水文站和水利工程的空間分布將研究區(qū)域劃分為若干個(gè)參數(shù)分區(qū)，可在各個(gè)參數(shù)分區(qū)中分別進(jìn)行參數(shù)率定工作；為更好地體現(xiàn)流域的空間差異性，同時(shí)能解決采用小網(wǎng)格單元帶來的計(jì)算負(fù)擔(dān)過重的問題，基于下墊面屬性，每個(gè)參數(shù)分區(qū)又被細(xì)化分為多個(gè)子流域和內(nèi)部空間單元。在研究流域被進(jìn)行空間劃分之后，運(yùn)用插值方法將降水、風(fēng)速、溫度、濕度、太陽輻射等[6]氣象信息也進(jìn)行空間離散。在對(duì)水循環(huán)過程進(jìn)行模擬時(shí)，該模型不僅支持自主開發(fā)的EasyDHM產(chǎn)流模型，同時(shí)也支持新安江模型、WetSpa模型、TOPMODEL等模型。匯流模擬可選用圣維南方法、馬斯京根法進(jìn)行計(jì)算。

EasyDHM中還嵌入了多套參數(shù)敏感性分析方法和參數(shù)優(yōu)化方法?，F(xiàn)有參數(shù)敏感性分析方法一般可分為局部分析法和全局分析法。一次二階矩法、隨機(jī)OAT法等均屬于局部分析法，主要是選取局部參數(shù)變化，分析對(duì)模型輸出造成的影響。全局分析法側(cè)重分析所有參數(shù)對(duì)模型模擬結(jié)果的影響，包括Monte Carlo方法、Latin-Hypercube模擬法、基于方差的分析法、LH-OAT方法等。在分布式水文模型的參數(shù)優(yōu)化中，優(yōu)化算法一般分為局部優(yōu)化算法和全局優(yōu)化算法。全面優(yōu)化的結(jié)果較局部優(yōu)化更有說服力，模型中可選用遺傳算法(GA)、多目標(biāo)粒子群算法(MOPSO)、動(dòng)態(tài)維度搜索算法(DDS)、SCE-UA算法等進(jìn)行優(yōu)化計(jì)算。對(duì)于目標(biāo)函數(shù)，模型中也有多種類型可供選擇。在此次增江流域的EasyDHM模型參數(shù)率定中，選擇DDS算法，以殘差平方和SSQ為目標(biāo)函數(shù)進(jìn)行參數(shù)率定。

2 模型建模

建立分布式水文模型需要的數(shù)據(jù)包括增江流域的DEM數(shù)據(jù)，數(shù)字河網(wǎng)，土地利用數(shù)據(jù)，土壤類型數(shù)據(jù)及流域的水文數(shù)據(jù)。本研究采用的DEM地表高程原始數(shù)據(jù)來自美國聯(lián)邦地質(zhì)調(diào)查局的HYDRO1K；土地利用源信息來自于《中國資源環(huán)境遙感宏觀調(diào)查與動(dòng)態(tài)研究》課題的研究成果數(shù)據(jù)的全國分縣土壤覆蓋矢量數(shù)據(jù)，采用增江流域2000年的土地利用類型數(shù)據(jù)；土壤基礎(chǔ)信息采用來自第二次全國土壤普查的匯總資料。

本研究中選取增江干支流上的2個(gè)水文控制站點(diǎn)和上游1個(gè)水庫參與流域模型構(gòu)建。在EasyDHM中，為了精細(xì)刻畫模擬流域水文過程，采用了子流域劃分、氣象數(shù)據(jù)插值、根據(jù)流域的下墊面屬性初步確定產(chǎn)匯流參數(shù)，再根據(jù)參數(shù)的敏感性分析結(jié)果對(duì)敏感性較高的參數(shù)進(jìn)行參數(shù)的自動(dòng)率定等步驟。

首先，對(duì)原始DEM數(shù)據(jù)進(jìn)行降低實(shí)際河網(wǎng)位置對(duì)應(yīng)的柵格的修正處理，通過給定最小水道集水面積閾值50 km2，將集水面積閾值超過給定閾值的柵格定為河網(wǎng)的單元，從而可得到數(shù)字河網(wǎng)水系，如圖1所示。EasyDHM模型可根據(jù)不同的土地利用類型，推求得到模型所需的部分產(chǎn)流、匯流參數(shù)初始值，而土壤質(zhì)地特征直接影響降雨的產(chǎn)流過程[7-9]。

圖1 增江流域DEM圖及水文氣象站點(diǎn)分布圖

在生成的數(shù)字河網(wǎng)的基礎(chǔ)之上，利用通用復(fù)雜流域、區(qū)域子流域劃分算法進(jìn)行子流域的劃分[11,12]。由于麒麟咀水文站之后只有水位站的數(shù)據(jù)而沒有水文站及雨量站的資料，缺乏實(shí)測(cè)流量記錄，無法進(jìn)行參數(shù)率定，同時(shí)麒麟咀以下河段受到潮位影響，流量會(huì)有周期性變化，不能單純地用水文模型來進(jìn)行預(yù)報(bào)。綜合以上因素，此次建模選擇了麒麟咀水文站作為流域出口。研究區(qū)域共劃被分為25個(gè)子流域。在本模型研究中，根據(jù)2個(gè)水文站、1個(gè)水庫，將增江流域劃分為3個(gè)參數(shù)分區(qū)，如圖2。天堂山水庫、渡頭水文站、麒麟咀(二)水文站為各區(qū)控制點(diǎn)，各個(gè)分區(qū)分別進(jìn)行參數(shù)率定工作。

圖2 增江流域的子流域及參數(shù)分區(qū)劃分示意圖

雨量數(shù)據(jù)來自17個(gè)雨量站的實(shí)測(cè)數(shù)據(jù)，氣象數(shù)據(jù)來自該流域內(nèi)及周邊氣象站的觀測(cè)數(shù)據(jù)，包括溫度、太陽輻射、濕度、風(fēng)速等，所有氣象及雨量數(shù)據(jù)均為點(diǎn)數(shù)據(jù)，需要進(jìn)行空間離散展布為面數(shù)據(jù)。增江流域模型中采用泰森多邊形法進(jìn)行空間展布[10]，根據(jù)距離柵格最小站點(diǎn)的數(shù)據(jù)作為此柵格的數(shù)據(jù)，再統(tǒng)計(jì)到子流域中，得到子流域的面平均氣象數(shù)據(jù)。

3 模型參數(shù)敏感性分析

EasyDHM模型參數(shù)包括產(chǎn)流參數(shù)和匯流參數(shù)[10,11]。本文采用LH-OAT全局分析法[16]對(duì)殘差平方和SSQ及平均流量兩個(gè)指標(biāo)對(duì)參數(shù)進(jìn)行多次敏感性分析，得出各參數(shù)的相對(duì)敏感度大小，從而可分析出各因子對(duì)模型影響程度。兩目標(biāo)函數(shù)分別如式(1)、(2)計(jì)算:

(1)

(2)

式中：Xsim,i為模型在時(shí)段i模擬出的水文站出口流量，m3/s；Xabs,i為時(shí)段i對(duì)應(yīng)的水文站實(shí)測(cè)流量，m3/s。

(3)

算法中采用公式(3)計(jì)算參數(shù)的敏感度，式中Si,j為參數(shù)ei第j個(gè)LH抽樣集合的敏感度，fi為參數(shù)ei的變化比例，M(e1,…,ei,…,ep)代表第j個(gè)LH抽樣集合的指標(biāo)函數(shù)，M[e1,…,ei(1+fi),…,ep]代表第i個(gè)LH抽樣集合上改變ei為ei(1+fi)后的指標(biāo)函數(shù)。根據(jù)相對(duì)敏感度大小，可分為4類：極高敏感(≥1.0)、高敏感(0.2～1.0)、中敏感(0.05～0.2)和低敏感(<0.05)。

模擬結(jié)果以Nash-Sutcliffe系數(shù)(以下簡稱Nash系數(shù))來進(jìn)行評(píng)定。Nash系數(shù)越接近1，則模擬效果越好。Nash系數(shù)計(jì)算公式為：

(4)

參數(shù)敏感性分析結(jié)果以渡頭站為例，如表1所示，圖3為渡頭站的不同目標(biāo)函數(shù)各參數(shù)敏感度分布圖。

表1 渡頭站參數(shù)敏感性分析結(jié)果

圖3 不同目標(biāo)函數(shù)各參數(shù)敏感度分布圖(渡頭站)

由表1及圖3可知，在以平均徑流量為目標(biāo)函數(shù)時(shí)，最大葉面積指數(shù)修正系數(shù)、最大冠層截留能力修正系數(shù)的敏感性較高，土壤凍結(jié)凍融參數(shù)、融雪積雪參數(shù)等表現(xiàn)出完全不敏感，這與增江流域上游區(qū)域內(nèi)的植被較多且地處亞熱帶氣候，冬季嚴(yán)寒氣候較少相符合。對(duì)于渡頭分區(qū)，我們可以發(fā)現(xiàn)對(duì)比以模擬SSQ、平均徑流量作為目標(biāo)函數(shù)進(jìn)行敏感性分析時(shí)，以SSQ為目標(biāo)函數(shù)的參數(shù)敏感性更高，說明了模型模擬效果不僅僅體現(xiàn)在模擬結(jié)果，還應(yīng)關(guān)注模擬值與實(shí)測(cè)值的相對(duì)或絕對(duì)誤差。同時(shí)，結(jié)合其他站點(diǎn)的參數(shù)敏感性分析結(jié)果，可知各個(gè)站點(diǎn)參數(shù)敏感值最大的均為土壤田間持水率修正系數(shù)、地表填洼能力修正系數(shù)也較高，說明在增江流域分布式水文模型中土地利用類型及土壤類型的影響較大；而從總徑流構(gòu)成上來看，土壤參數(shù)敏感度及地表徑流參數(shù)與地下水模擬參數(shù)的敏感度相比都較高，說明在增江流域中，地下水的影響較小，地下徑流穩(wěn)定，徑流量的變化主要體現(xiàn)在地表徑流及壤中流中；觀察到子流域主河道長度修正系數(shù)均為極端高敏感系數(shù)，說明，在增江流域中除了產(chǎn)流模塊影響較大之外，匯流模塊也對(duì)模擬結(jié)果有著重要影響。

4 模型參數(shù)率定

由于分布式水文模型參數(shù)數(shù)目較多，對(duì)所有參數(shù)同時(shí)優(yōu)化效率太低，因而僅優(yōu)化經(jīng)參數(shù)敏感性分析所選出的極端高敏感度參數(shù)和高敏感度參數(shù)(統(tǒng)一稱為敏感參數(shù))，利用DDS算法[7]，先進(jìn)行全局搜索，隨著迭代次數(shù)的增加逐漸轉(zhuǎn)向近似局部的搜索，以SSQ為目標(biāo)函數(shù)對(duì)敏感度高的參數(shù)進(jìn)行優(yōu)化率定，得出優(yōu)化結(jié)果。

由于天堂山水庫站資料的缺乏，本文僅以渡頭站、麒麟咀站為例來說明模擬參數(shù)率定的效果。以2008年至2010年為率定期，選取三年內(nèi)多場(chǎng)次洪水采用日時(shí)間步長水進(jìn)行組合率定，得出參數(shù)最優(yōu)值，如表2及表3。麒麟咀水文站選取2012及2013兩年作為驗(yàn)證期，而由于渡頭水文站缺失2012年的氣象數(shù)據(jù)，因此選定2007年及2013年驗(yàn)證進(jìn)一步驗(yàn)證各參數(shù)分區(qū)率定結(jié)果的可靠性。表4為率定期及驗(yàn)證期的部分洪水場(chǎng)次模擬結(jié)果。

表2 渡頭站產(chǎn)匯流參數(shù)最優(yōu)值

表3 麒麟咀站產(chǎn)匯流參數(shù)最優(yōu)值

從表4模擬結(jié)果可以看出，水文站的率定結(jié)果都較好，挑選的洪水場(chǎng)次的Nash系數(shù)均在0.5以上。麒麟咀水文站的洪水模擬效果略差，從2008年的麒麟咀水文站的模擬結(jié)果圖可以看出，模擬結(jié)果和實(shí)測(cè)值存在滯后效果，實(shí)測(cè)值的徑流流量相對(duì)模擬值滯后將近一天的時(shí)間，麒麟咀水文站處于流域的中下游地區(qū)，在增江流域中建有多個(gè)攔水壩，而在徑流過程中攔水壩影響了匯流過程，位于模擬區(qū)域的出口處的麒麟咀水文站所測(cè)得的徑流流量，導(dǎo)致實(shí)測(cè)值滯后于模擬值。且上游建有多個(gè)中小型水庫，徑流流量受水庫調(diào)蓄影響較大。

表4 各站場(chǎng)次洪水模擬結(jié)果的Nash效率系數(shù)

續(xù)表4 各站場(chǎng)次洪水模擬結(jié)果的Nash效率系數(shù)

為了繼續(xù)驗(yàn)證模型中率定的參數(shù)對(duì)此流域的適應(yīng)性，根據(jù)各個(gè)站點(diǎn)的水文數(shù)據(jù)，以2013年為例，驗(yàn)證期各參數(shù)分區(qū)的結(jié)果見表5。挑選的洪水場(chǎng)次的模擬見表4。從表4中2013年麒麟咀站的模擬結(jié)果和驗(yàn)證期2008年時(shí)的結(jié)果相似，存在滯后效果，證明猜想正確，流域中的攔水壩影響著匯流過程，造成了模擬結(jié)果的偏差。同時(shí)從各參數(shù)分區(qū)的驗(yàn)證結(jié)果可以看出，驗(yàn)證期間的Nash效率系數(shù)均達(dá)到0.7以上，說明模擬效果較好，分布式水文模型EasyDHM在增江流域具有良好的實(shí)用性，因此，可進(jìn)一步用于增江流域?qū)嶋H問題的研究。

表5 模型各參數(shù)分區(qū)模擬結(jié)果

5 結(jié) 語

本文針對(duì)增江流域降水較多，下墊面情況復(fù)雜的特點(diǎn)，在對(duì)以收集的水文、氣象、地形等數(shù)據(jù)的整理編輯之后，建立了增江流域EasyDHM分布式水文模型。在使用分布式水文模型的通過敏感性分析，識(shí)別出影響該流域的主要參數(shù)因子為：土壤添加持水率修正系數(shù)，子流域主河道長度修正系數(shù)，不透水面積修正系數(shù)等?？傮w上，各參數(shù)分區(qū)模擬效果較為理想，說明分布式水文模型EasyDHM在增江流域應(yīng)用良好，可利用該模型實(shí)現(xiàn)對(duì)流域的水雨情分析，為水資源的管理提供支持，使得水資源利用更加合理，從而緩解地區(qū)資源性缺水的問題。經(jīng)比較，麒麟咀分區(qū)由于流域內(nèi)的攔水壩影響，導(dǎo)致實(shí)測(cè)值相對(duì)模擬結(jié)果存在滯后現(xiàn)象，將在之后的應(yīng)用中重點(diǎn)展開對(duì)此影響因子的研究。