汪 偉, 方秀琴, 杜曉彤, 吳陶櫻, 吳小君

(河海大學(xué) 地球科學(xué)與工程學(xué)院, 南京 211100)

水文模型是為了模擬水循環(huán)過程而構(gòu)建的。水文模型通常分為集總式水文模型[1]和分布式水文模型[2-4]，前者的代表性模型有NAM(丹麥語(yǔ)，“降雨徑流模型”的縮寫)模型和新安江模型，后者的代表模型則有TOPMODEL(Topographic Hydrologic Model )模型[5]、SWAT(Soil and Water Assessment Tool)模型[6]、DHSVM(Distributed Hydrology Soil Vegetation Model)[7]模型以及VIC(Variable Infiltration Capacity)模型等。得益于計(jì)算機(jī)技術(shù)的不斷發(fā)展進(jìn)步，將下墊面劃分為大量網(wǎng)格的分布式水文模型具有充分考慮流域內(nèi)降雨和土地覆蓋在空間上分布的不均勻性特點(diǎn)，因此相較于用單純物理參數(shù)來(lái)描述水文過程變化的集總式水文模型而言，分布式水文模型可以充分描述流域內(nèi)降雨和下墊面要素的空間變化對(duì)徑流產(chǎn)生的影響，從而在水文模擬過程中更具優(yōu)勢(shì)。除此以外，分布式水文模型還更容易與GIS(Geographic Information Science)、RS(Remote Sensing)相結(jié)合。隨著GIS與RS技術(shù)的不斷進(jìn)步，分布式水文模型的基礎(chǔ)數(shù)據(jù)來(lái)源更加的廣泛也更加可靠，從而能更好的描述客觀世界水文物理過程。

柳江是珠江水系西江的第二大支流，汛期流域內(nèi)洪水災(zāi)害發(fā)生較為頻繁。但由于研究區(qū)研究資料的限制以及其他原因，目前針對(duì)柳江流域的水文水資源方面的研究相對(duì)較少。因此針對(duì)柳江流域的水文物理研究對(duì)于柳江流域的水資源管理、洪水災(zāi)害防治均具有十分重要的意義[7-9]。李澤鋒等利用收集到的水文、氣象資料對(duì)柳江流域不同時(shí)間尺度的降雨徑流利用支持向量機(jī)模型進(jìn)行了模擬分析[10]。劉歡等基于Green-Ampt模型量化空氣阻力的影響，同時(shí)結(jié)合分布式水文模型WEP-L模型對(duì)柳江流域徑流進(jìn)行了模擬[11]，結(jié)果表明：通過量化空氣阻力的影響能夠改善柳江流域的水文模擬效果，但改善效果不夠顯著。本研究結(jié)合現(xiàn)有觀測(cè)資料，利用分布式水文模型VIC模型對(duì)柳江流域水文物理過程進(jìn)行模擬，探索分布式水文模型(VIC模型)在柳江流域的適用性，為柳江流域水資源管理、洪水災(zāi)害的防治提供參考依據(jù)。

1 研究區(qū)概況

柳江流域地處東經(jīng)107°30′—111°15′，北緯24°25′—27°10′，流域面積約5.7萬(wàn)km2，橫跨湘、桂、黔三省。柳江發(fā)源于貴州省獨(dú)山縣，河流全長(zhǎng)約77 km，其主要支流有寨嵩河、古宜河、金城江、龍江、都柳江等。柳江水系屬于樹枝型流域，上游河道灘多且河水湍急，中下游水勢(shì)較上游相對(duì)平緩。

柳江流域地勢(shì)整體上為呈現(xiàn)西北高而東南低(圖1)，上游除河源地區(qū)大多屬于變質(zhì)巖組成的峽谷地形為主，海拔多在1 000 m以上；中下游以喀斯特地貌為主，山地、丘陵的占比較大，平均海拔高度為500～1 000 m，耕地集中，人口稠密。柳江流域年平均氣溫為17～20℃，年降雨量在1 200～1 800 mm，降水時(shí)空分布極不均勻：在時(shí)間上主要集中在4—9月，占全年的70%以上；空間上，流域降水呈現(xiàn)出東南高而西北低的趨勢(shì)。

圖1 柳江流域及其主要支流分布

2 模型和數(shù)據(jù)

2.1 VIC模型

VIC模型是基于物理機(jī)制的分布式水文模型，由Wood[12]等在1992年提出，經(jīng)梁旭等人[13-14]改進(jìn)。由VIC-2L模型逐漸發(fā)展到如今普遍采用的VIC-3L模型。相較于其他水文模型，VIC模型的主要特點(diǎn)[15]有：(1) 能夠同時(shí)進(jìn)行陸—?dú)忾g能量平衡以及水量平衡的模擬。通過全局參數(shù)文件的參數(shù)設(shè)置可以自由決定是否開啟能量平衡模擬，從而彌補(bǔ)了傳統(tǒng)水文模型對(duì)能量交換過程描述的不足；(2) 利用土壤下滲能力參數(shù)B來(lái)描述次網(wǎng)格內(nèi)土壤的不均勻性對(duì)產(chǎn)流的影響；(3) 產(chǎn)流過程的計(jì)算方式與新安江模型相同，同時(shí)考慮兩種產(chǎn)流機(jī)制(蓄滿產(chǎn)流和超滲產(chǎn)流)。除此以外，VIC模型代碼開源，易于做模型的改進(jìn)與發(fā)展。

VIC模型由陸面過程模型和匯流模型組成。陸面過程模型將研究區(qū)劃分為若干大小相同的網(wǎng)格(本研究將柳江流域劃分為676個(gè)網(wǎng)格)，將每個(gè)網(wǎng)格對(duì)應(yīng)的理化性質(zhì)數(shù)據(jù)，作為驅(qū)動(dòng)VIC模型的輸入數(shù)據(jù)。VIC模型采用新安江模型的產(chǎn)流計(jì)算方式：假設(shè)上層土壤蓄水能力在模型網(wǎng)格單元內(nèi)或者研究區(qū)域內(nèi)變化，土壤蓄水能力的空間變化主要是土層的深度和土壤特性的各向異性產(chǎn)生的，土壤蓄水能力的空間變化特性表示如下[16-17]：

i=im[1-(1-C)1/b]

(1)

式中：i為網(wǎng)格的土壤蓄水能力；im為最大土壤蓄水能力；C為土壤蓄水能力小于或者等于i部分的面積比例；b為土壤蓄水能力形狀特征參數(shù)，是土壤蓄水能力空間變化的表征，可以表示土壤特性的空間變異性。

VIC模型對(duì)基流的描述采用的是ARNO模型[18]，默認(rèn)基流是發(fā)生在最下層土壤中。各土層間的水汽通量服從達(dá)西(Darcy)定律從而計(jì)算出基流。其計(jì)算公式為：

(2)

匯流模型[19]則是以陸面過程的輸出數(shù)據(jù)作為輸入數(shù)據(jù)，使用單位線文件進(jìn)行坡面匯流，使用線性圣維南方程來(lái)進(jìn)行河網(wǎng)匯流。其匯流方案是先演后合，即假設(shè)每個(gè)網(wǎng)格的產(chǎn)流量對(duì)流域出口徑流量的貢獻(xiàn)相互獨(dú)立，先計(jì)算每個(gè)網(wǎng)格到流域出口斷面的流量過程，然后疊加，最終得到流域的模擬徑流。

2.2 數(shù) 據(jù)

VIC模型所需的原始輸入數(shù)據(jù)及其描述見表1，所有原始數(shù)據(jù)均通過ArcMap平臺(tái)重采樣為應(yīng)用在柳江流域的模型分辨率9 km×9 km。

表1 VIC模型輸入數(shù)據(jù)

針對(duì)不同的數(shù)據(jù)類型有不同的數(shù)據(jù)處理原則，分別為：

(1) 氣象驅(qū)動(dòng)數(shù)據(jù)：選取在柳江流域周圍20個(gè)站點(diǎn)(圖1)的降水、最高氣溫、最低氣溫、風(fēng)速4個(gè)氣象要素，同時(shí)使用IDW(Inverse Distance Weighted)插值方法得到柳江流域各個(gè)網(wǎng)格的氣象數(shù)據(jù)，所得出數(shù)據(jù)作為運(yùn)行VIC模型的氣象輸入數(shù)據(jù)。

(2) 土地覆蓋數(shù)據(jù)：統(tǒng)計(jì)網(wǎng)格內(nèi)所有的土地覆蓋類型，計(jì)算網(wǎng)格內(nèi)每一種土地覆蓋類型占比，同時(shí)附上該網(wǎng)格內(nèi)每一種類型土地覆蓋的理化性質(zhì)數(shù)據(jù)。由圖2A可知柳江流域的主要土地覆蓋類型為草地、林地和耕地。

(3) 土壤性質(zhì)數(shù)據(jù)：其數(shù)據(jù)處理的原則與土地覆蓋數(shù)據(jù)大致相同。主要區(qū)別是根據(jù)流域網(wǎng)格的劃分，在流域內(nèi)任意一個(gè)網(wǎng)格內(nèi)若該類土壤占比最大，則使用該類土壤的理化性質(zhì)代表該網(wǎng)格的所有土壤理化性質(zhì)。由圖2B可知柳江流域的土壤類型以黏土、壤土和粉質(zhì)壤土為主。

圖2 柳江流域土地覆蓋和土壤性質(zhì)數(shù)據(jù)

3 結(jié)果與分析

3.1 參數(shù)率定

本研究將柳江流域VIC模型的模擬過程劃分為率定期和驗(yàn)證期兩個(gè)階段。其中1981—1989年設(shè)為模型的率定期，2007—2011年為模型驗(yàn)證期。選取3個(gè)目標(biāo)函數(shù)[20]作為評(píng)價(jià)VIC模型在柳江流域徑流模擬效果的評(píng)價(jià)指標(biāo)，分別是Nash-Sutcliffe效率系數(shù)(NS)[21-22]、相關(guān)性系數(shù)(r)以及相對(duì)偏差(BIAS)，目標(biāo)函數(shù)的具體計(jì)算公式見表2。

其中NS通過對(duì)比在日尺度和月尺度上模型模擬徑流量與實(shí)際徑流量的差異來(lái)反映本次模擬過程的效果，其值越接近1表明本次模擬擬合程度越好。r用來(lái)評(píng)價(jià)模擬徑流量與實(shí)測(cè)徑流量之間的相關(guān)性，其值越接近于1表示模型模擬的數(shù)據(jù)越接近實(shí)測(cè)徑流數(shù)據(jù)。BIAS用來(lái)描述模型模擬徑流量與實(shí)際徑流量在總量上的相似程度，其值越接近于0，表明模型模擬的徑流量越接近于實(shí)際徑流量。以上3個(gè)評(píng)價(jià)指標(biāo)分別從不同的角度對(duì)模型模擬效果進(jìn)行評(píng)價(jià)。

表2 模擬效果評(píng)價(jià)指標(biāo)及其公式

本研究通過編寫matlab腳本的方式完成VIC模型的參數(shù)率定和驗(yàn)證工作。主要思路如下：將VIC模型的6個(gè)參數(shù)(B，Ds，Dm，Ws，d2，d3)設(shè)置成待率定的變量，確定好待率定參數(shù)所有可能的組合方式，作為輸入?yún)?shù)按照每種組合方式生成VIC模型土壤參數(shù)文件。土壤參數(shù)文件生成后調(diào)用VIC模型源程序生成陸面過程結(jié)果文件。將該次結(jié)果文件作為輸入文件運(yùn)行匯流模型源程序，最終產(chǎn)生該次參數(shù)組合方式的徑流量模擬結(jié)果，將該次徑流結(jié)果數(shù)據(jù)與實(shí)際徑流數(shù)據(jù)通過目標(biāo)函數(shù)進(jìn)行比較，記錄每次模擬的目標(biāo)函數(shù)結(jié)果，尋找目標(biāo)函數(shù)最優(yōu)值，直到所有的參數(shù)組合方式都運(yùn)行完畢。本研究采用Rosenbrock[23]法來(lái)進(jìn)行模擬，即每次只劃定一個(gè)參數(shù)的范圍，同時(shí)確定其參數(shù)變化步長(zhǎng)，剩余5個(gè)參數(shù)用唯一值代替，模擬過程不變，直到所有的參數(shù)都率定完畢。柳江流域VIC模型參數(shù)率定結(jié)果見表3。

表3 VIC模型率定的參數(shù)物理意義及其取值

3.2 模型模擬結(jié)果

3.2.1 基于流域斷面的模擬效果分析 表4給出了柳江流域VIC模型在日、月兩個(gè)時(shí)間尺度上的率定期和驗(yàn)證期模擬結(jié)果。在日尺度上，率定期內(nèi)的NS，r分別為0.74，0.86，均超過了0.7。在徑流量上的相對(duì)偏差BIAS為0.10。驗(yàn)證期內(nèi)的NS，r分別為0.70，0.80，均達(dá)到或者超過了0.70，相對(duì)偏差BIAS為0.06。雖然較率定期的模擬效果均有不同程度的下降，但在徑流量上的模擬甚至要優(yōu)于率定期的BIAS的值，且達(dá)到了NS>0.5同時(shí)BIAS的絕對(duì)值小于0.15的要求，說(shuō)明柳江流域VIC模型的模擬效果符合預(yù)期結(jié)果。在月尺度上無(wú)論是率定期還是驗(yàn)證期，對(duì)比日尺度的模擬結(jié)果，其模擬效果均要優(yōu)于日尺度。其中率定期的NS，r分別為0.86，0.93，而在徑流量上則和日尺度的模擬效果相差不大，BIAS的值為0.10。而在驗(yàn)證期內(nèi)月尺度的模擬效果與率定期相差不多，NS出現(xiàn)了小幅上漲，為0.89，而相關(guān)性系數(shù)和相對(duì)偏差出現(xiàn)了小幅下降趨勢(shì)，其值分別為0.93，0.07。

表4 柳江流域模擬結(jié)果

圖3和圖4是在日尺度和月尺度上柳江流域模擬和實(shí)測(cè)流量散點(diǎn)圖，從這四張圖中可以看出不管是率定期還是驗(yàn)證期，絕大部分的散點(diǎn)均比較均勻分布在對(duì)稱線的兩側(cè)，但是在日尺度流量模擬過程(圖3)中，不管是率定期還是驗(yàn)證期，都有零星幾個(gè)散點(diǎn)位于偏左上的位置，說(shuō)明在日尺度流量模擬過程中VIC模型對(duì)于洪峰的峰值模擬偏小，但總體上還處于比較均勻的分布狀態(tài)。而在月尺度流量模擬過程(圖4)中，率定期散點(diǎn)分布較為集中且總體處于對(duì)稱線偏下一點(diǎn)位置說(shuō)明月尺度模擬徑流量稍稍偏大。驗(yàn)證期的散點(diǎn)分布稍顯分散但總體均勻處在對(duì)稱線的左右兩側(cè)，這也與圖3中的數(shù)據(jù)結(jié)果相一致。

圖3 柳州站實(shí)測(cè)與模擬日流量散點(diǎn)圖

圖4 柳州站實(shí)測(cè)與模擬月流量散點(diǎn)圖

圖5和圖6為柳州站率定期和驗(yàn)證期模擬與實(shí)測(cè)日值流量過程線對(duì)比圖。從圖中我們可以看出，無(wú)論是率定期還是驗(yàn)證期，VIC模型都能較好的模擬出柳江流域?qū)崪y(cè)徑流過程線，尤其是洪峰出現(xiàn)時(shí)間和退水時(shí)間以及秋冬季節(jié)枯水期水量的模擬，都基本與實(shí)測(cè)徑流過程保持一致。在圖上可以清楚看出“1988.8”大洪水災(zāi)害的時(shí)間變化過程。但是在洪峰流量的模擬上，VIC模型模擬的數(shù)值都出現(xiàn)了不同程度的偏低，這與上文圖3散點(diǎn)圖的日流量模擬過程中VIC模型對(duì)于洪峰的峰值模擬偏小的結(jié)論相一致。這可能和模型本身參數(shù)的局限性以及所收集的氣象資料的局限性有關(guān)。此外無(wú)論是在率定期還是驗(yàn)證期，在模擬時(shí)段剛開始的時(shí)候都出現(xiàn)了模擬徑流量的值要遠(yuǎn)大于實(shí)測(cè)徑流量的值，此后隨著時(shí)間的推移又快速恢復(fù)的情況。這可能是由于模擬過程中出現(xiàn)的“邊緣效應(yīng)”所導(dǎo)致的。

3.2.2 基于網(wǎng)格的流域模擬效果分析 利用VIC模型的輸出參數(shù)文件，本研究提取了率定期模擬過程輸出文件中的4個(gè)要素(包括多年平均降水(P)、蒸散發(fā)(E)、地表徑流(R)、地下基流(B))進(jìn)行可視化處理，如圖7所示。該降水?dāng)?shù)據(jù)和驅(qū)動(dòng)VIC模型的輸入降水?dāng)?shù)據(jù)相同，所以降水分布空間變化比較均勻，流域年均降水量為1 069.3～1 723.6 mm，空間上表現(xiàn)為東南多、西北少。蒸散發(fā)數(shù)據(jù)和地表徑流數(shù)據(jù)的空間變化主要和土壤類型有關(guān)，同時(shí)結(jié)合圖2的土壤類型分布可以得出：柳江流域內(nèi)壤土和壤質(zhì)沙土地區(qū)的蒸發(fā)量要大于黏土地區(qū)的蒸發(fā)量，相應(yīng)的壤土和壤質(zhì)沙土地區(qū)地表徑流普遍要小于黏土地區(qū)的地表徑流。柳江流域的蒸發(fā)量主要在408.4～692.9 mm之間，地表徑流在219～680.1 mm之間變動(dòng)。地下基流的變化則不太明顯，空間分布相對(duì)均勻。其值主要在331.1～545.8 mm之間變動(dòng)。

圖5 柳州站率定期(1981-1989)實(shí)測(cè)與模擬日流量過程對(duì)比

圖6 柳州站驗(yàn)證期(2007-2011)實(shí)測(cè)與模擬日流量過程對(duì)比

圖7 柳江流域主要水文要素空間分布

根據(jù)柳江流域VIC模型的輸出要素(P，E，R，B)，結(jié)合流域水量平衡公式：

ΔW=P-E-R-B

(3)

其中ΔW為時(shí)段內(nèi)流域水量變化，得出流域水量平衡相對(duì)誤差ΔR為：

(4)

根據(jù)該公式可得出基于網(wǎng)格的柳江流域多年平均水量平衡相對(duì)誤差。由圖8我們可以看出柳江流域的水量平衡相對(duì)誤差均在1.2%以下。且在整個(gè)空間上分布較為均勻，說(shuō)明VIC模型對(duì)柳江流域的模擬效果達(dá)到與預(yù)期。

3.3 模型應(yīng)用

利用建立好的柳江流域VIC模型框架，以1981—2011年的氣象數(shù)據(jù)作為驅(qū)動(dòng)模型運(yùn)行的輸入數(shù)據(jù)，得到對(duì)應(yīng)計(jì)算時(shí)段的模型輸出要素，提取其中的P，E，R，B水文要素。以年為單位，結(jié)合Mann-Kendall趨勢(shì)檢驗(yàn)方法[24-28]對(duì)柳江流域這4個(gè)水文要素的變化情況進(jìn)行趨勢(shì)性分析。

圖8 柳江流域相對(duì)誤差空間分布

M-K趨勢(shì)檢驗(yàn)是指對(duì)樣本的相對(duì)大小進(jìn)行比較，并對(duì)該樣本是否有某種變化趨勢(shì)進(jìn)行定量的判斷。對(duì)一個(gè)隨機(jī)變量組X(x1,x2,…，xn)，M-K趨勢(shì)檢驗(yàn)方法的公式介紹見表5。

表5 M-K趨勢(shì)檢驗(yàn)方法公式及其描述

本文參考曹潔萍[29]等的研究，當(dāng)Z的絕對(duì)值在大于1.64時(shí)表示通過了置信度為95%的顯著性檢驗(yàn)。

柳江流域的M-K檢驗(yàn)結(jié)果見圖9：柵格顏色的深淺代表變化趨勢(shì)的強(qiáng)弱，帶矩形框的柵格表示該柵格通過了顯著性檢驗(yàn)，有顯著的增加或減小趨勢(shì)。由圖9可以看出地表徑流通過顯著性檢驗(yàn)的柵格最多，共480個(gè)，占流域總面積的71%，覆蓋了柳江流域除東北方向的絕大部分地區(qū)。降水通過顯著性檢驗(yàn)的柵格次之，共328個(gè)，占流域總面積的48%。蒸發(fā)量通過顯著性檢驗(yàn)的柵格共有162個(gè)，占流域總面積的24%，且主要分布在柳江流域的西部和西北部，地下基流則基本上沒有明顯變化的趨勢(shì)，其Z值大于1.64的柵格僅8個(gè)，占比1%。

圖9 柳江流域柵格M-K檢驗(yàn)結(jié)果

4 結(jié)論與展望

VIC模型對(duì)于柳江流域水文物理過程的模擬效果無(wú)論是在時(shí)間上還是在空間上都表現(xiàn)出了良好的適用性。時(shí)間上：無(wú)論是驗(yàn)證期還是率定期，在日尺度上其效率系數(shù)NS的值都在0.7以上，相對(duì)偏差BIAS的絕對(duì)值均小于0.11，月尺度上VIC模型的模擬效果較日尺度明顯更優(yōu)?？臻g上：柳江流域所有網(wǎng)格水量平衡相對(duì)誤差的值均小于1.2%。二者皆表明了VIC模型在柳江流域良好的模擬效果，可以為柳江流域水資源管理以及洪災(zāi)防治提供一定的依據(jù)。

針對(duì)柳江流域VIC模型的輸出要素(P，E，R，B)進(jìn)行了M-K趨勢(shì)性檢驗(yàn)，結(jié)果表明柳江流域有超過70%的柵格地表徑流有明顯增加趨勢(shì)，而地下基流普遍無(wú)顯著變化趨勢(shì)，降水和蒸散發(fā)具有顯著增加趨勢(shì)的柵格均不到總柵格的50%。

在研究的過程中我們發(fā)現(xiàn)VIC模型在模擬的起步階段會(huì)出現(xiàn)模擬值遠(yuǎn)大于實(shí)測(cè)值的“邊緣效應(yīng)”現(xiàn)象，其原因可能與模型本身的產(chǎn)流原理以及參數(shù)自身的局限性有關(guān)，未來(lái)可在這一方面進(jìn)行進(jìn)一步的探究與改進(jìn)。此外本文僅對(duì)柳江流域水文要素的時(shí)空分布進(jìn)行了模擬研究；其他諸如土地利用變化、水文參數(shù)變化以及其他要素變化對(duì)于模擬效果的影響并未進(jìn)行深入研究。