VIC模型在西江流域的水文模擬及其應用

杜婷婷, 郭夢京, 張晉梅, 田世野

(1.綿陽職業技術學院, 四川 綿陽621000; 2.西安理工大學 水利水電學院, 西安 710077;3.江西省水利科學研究院, 南昌330029; 4.江西省鄱陽湖流域農業資源與生態重點實驗室, 南昌 330045)

水文模型[1]是用數學公式和物理性質來描述流域自然水循環和能量交換過程，隨著計算機科學的不斷發展，水文模型被廣泛的應用在水文模擬預測、旱澇災害預警、流域水資源管理等各個方面。得益于不斷進步的地理信息系統(GIS)[2]技術與遙感技術(RS)[3]，結合了GIS與RS技術的分布式水文模型不但考慮了傳統水文模型所需要的降雨數據還重點考慮到了下墊面以及氣候的變化對于流域水文生態的影響，使得分布式水文模型在流域自然水循環模擬過程中更具優勢。

珠江流域主要支流有三條，分別是東江、北江以及西江，其中西江是珠海、澳門等地區城市用水最重要的上游來水源[4]。縱觀前人的研究，對于珠江流域的水文研究主要集中在流域面積較小的北江和東江，李建慶，羅顯剛等人基于不同土地利用情景利用SWAT模型對北江流域進行了水文模擬，結果表明：城市化以及退耕還林還草等土地利用變化都會對流域自然水文過程造成較大影響[5]，綦昕瑤，劉貴花等利用IHACRES水文模型對東江流域徑流變化的原因進行分析時發現：氣候變化是引起東江流域徑流變化的主要因素，人類活動是次要因素[6]。而在西江流域水文方面的研究由于研究資料的匱乏以及其他原因相對較少，廖衛紅，雷曉輝等利用分布式水文模型EasyDHM對西江流域徑流進行模擬，結果表明該模型在西江流域具有良好的適用性[7]；曾凌，熊立華等利用DDRM(DEM-based distributed rainfall-runoff model)模型模擬了西江流域的土壤濕度的時空分布，結合ASCAT衛星遙感反演土壤濕度產品作為對比，結果表明DDRM模型和衛星遙感土壤濕度指數具有良好的時間一致性和空間一致性[8]；趙胤懋，廖衛紅等以西江為研究區利用地面雨量站的數據對比分析了CMORPH降水產品精度，分析結果表明：CMORPH降水產品具有較高的產品精度，能夠替代地面雨量站的數據[9]；吳志勇，林青霞等選取1951—2010年水文站徑流數據，利用SRDI干旱指數研究西江流域水文干旱的空間變化特征，結果表明：西江流域輕旱發生范圍有擴大趨勢，極旱發生范圍有減小趨勢[10]。綜合來看針對西江的水文研究大都集中在某一點或者某一產品的數據精度上，對于西江流域整體的系統性水文研究則相對較少。本研究利用遙感技術獲取水文分析的基礎數據、利用GIS空間分析技術對水文數據進行轉換與處理，以西江作為目標研究流域，探究分布式水文模型VIC模型在西江流域的適用性情況，以期為西江流域水資源科學管理及洪水災害防治提供一定的科學依據。

1 研究區概況

西江[11]地處北緯21°31′—26°49′、東經102°14′—114°48′，流域橫跨云南、廣西以及貴州，集水面積約35萬km2，占整個珠江流域總面積超過七成，年均徑流量超過2 300億m3，主要支流包括左江、右江、黔江、紅水河、郁江、桂江以及柳江等。

西江流域在氣候上的垂直差異和水平差異相當明顯；不同海拔，不同經緯度地區之間差異明顯。流域氣溫差異明顯，年平均氣溫常年保持在14～23℃。流域年均降水無論在空間上還是時間上分布都極不均勻。時間上：汛期(4—9月)水量占到全年總徑流量的70%；空間上則呈現出東南多西北少的趨勢，年均降雨量在1 000～2 200 mm。

2 VIC模型概述、數據來源及研究方法

2.1 VIC模型簡介

最初版本的VIC[12-14]模型又叫VIC-2L模型，其中的2L代表的是對土壤的分層，即將土壤分成了上下兩層。隨著越來越多的水文工作者對VIC模型研究水平不斷提高，VIC版本隨之不斷更迭，將第二層土壤進行了分割，劃出頂薄層，新增了對于土壤表層的動態描述，提高了土壤濕度的模擬精度，最終形成了如今被廣大水文研究工作者廣泛采用的VIC-3L模型。VIC-3L模型考慮了裸土以及不同植被覆蓋類型對于降雨的動態響應。

VIC模型在產流方式上參考新安江模型同時考慮了兩種產流方式(蓄滿產流和超滲產流[15-17])，在模擬歷時任意一個Δt內，該段時間內的降雨量會被VIC模型產流過程拆分為蓄滿產流R1、超滲產流R2以及下滲到土壤的降雨ΔW3個部分。Wt是t時刻的土壤含水量，而不同性質土壤的空間分布的不均勻性導致降雨的入滲產流的方式也隨著空間變化而變化，從蓄滿產流變成超滲產流，網格內降雨在土壤中的入滲能力也隨土壤性質的變化而變化，具體如下式：

f=fm[1-(1-C)1/B]

(1)

式中:f表征不同土壤的降雨入滲能力在空間上的變化;fm代表流域內網格降雨量的最大入滲能力大小;C表征的是流域內降雨的入滲能力值小于或者等于f的網格面積占流域總網格面積的比例；B的意義與上式類似,是表征流域內土壤入滲能力參數。

P=R1(y)+R2(y)+ΔW(y)

(2)

式中：P為在歷時內的流域降雨量總和；R1為蓄滿產流量；R2為超滲產流量；ΔW為入滲到土壤的雨量，且y=R1(y)+ΔW(y)。其中蓄滿產流和超滲產流加上入滲到土壤的總水量ΔW可以用y式來進行表示,具體如下:

(3)

(4)

(5)

式中:fmm為流域土壤平均入滲能力,其他參數的意義與前面公式說明相同,經過上式的計算后就可以得到地下基流、地表徑流以及三層土壤的含水量。

VIC模型的一大特點是通過全局參數文件的參數設置可以自由決定是否開啟能量平衡模擬從而在研究流域內同時進行能量平衡模擬以及水量平衡模擬。此外VIC模型將流域內氣候、土壤屬性、地形地貌以及植被的綜合作用以數學公式的方式進行集成處理，使得VIC模型具備扎實的數學意義和較高的物理基礎。最后VIC模型源碼開源，方便后續研究人員做模型的改進與發展。

2.2 數據來源與處理

VIC模型在應用過程中，模型的預處理是一個必不可少的環節，借助Arcmap平臺，基于STRM 90 m高程數據對西江流域進行劃分，以2 km×2 km模型分辨率在Arcmap軟件中利用西江流邊界將西江流域劃分為723個網格單元。

本研究所使用的部分VIC模型輸入數據來源如下：(1) 地表覆蓋數據來源于馬大1 km分辨率全球地表覆蓋數據集。該數據集的處理方法為：在流域所劃分好的所有網格中，統計網格內所有的地表覆蓋類型，同時將所有類型的地表覆蓋屬性信息賦值到對應網格上。(2) 土壤質地數據采集自世界土壤數據庫；土壤質地數據集的處理原則與地表覆蓋數據大體類似，區別之處在于網格內占比最大的土壤質地的理化性質數據代表整個網格的理化性質屬性數據。(3) 氣象驅動數據來源中國地面累年日值數據集。該數據集包含了全國699個氣象站點的70 a日值氣象數據。本研究針對西江流域提取了流域周圍共88個氣象站點的氣象數據進行研究，其數據處理原則為：根據流域周圍共88個氣象站點的氣象要素數據，從中提取降水、最高氣溫、最低氣溫、風速4個子要素組成新的氣象要素數據集，作為運行VIC模型的氣象驅動數據。(4) 用于VIC模型模擬結果的率定和驗證所需要的實際徑流量數據來自于武宣水文站1981—1989年逐日徑流量數據。

根據已經收集的數據可知，西江流域土壤主要以沙壤土、黏土以及沙質黏土為主，其中又以沙壤土占據最主要部分。土地覆蓋類型主要草地、耕地和林地3種基本覆蓋類型

2.3 研究方法

M-K[18-19](Mann-Kendell)趨勢檢驗方法是分析水文氣象要素在時間上變化特點常用的一種分析方法，該方法的優勢在于用于檢驗的時間要素序列樣本無須遵循某一特定的分布規律，在進行變化特點分析時極少會受到異常值的影響，而且計算過程簡單且易于分析。M-K趨勢檢驗方法的原理是假設有一個時間序列Xn，將每一個時間序列的子元素進行對比，根據其數值的大小來判斷該時間序列是否具有某種變化趨勢，具體見表1。

表1 M-K趨勢檢驗方法公式及其描述

表1中Z值的正負表示該時間要素序列變化趨勢為增加趨勢還是減小趨勢，當值大于零時為增加趨勢，小于零表示該要素序列表現出減小趨勢，Z的值的大小用來定量表征該時間序列變化趨勢是否顯著，參考他人研究[20]當Z>|±1.64|，則該要素序列具有某一顯著的變化趨勢，正值具有顯著增加趨勢，反之則有減小趨勢。

3 結果與分析

3.1 參數率定

VIC模型土壤參數眾多，絕大部分的參數有其固定物理意義[21]，通過不同土壤屬性數據的獲取，得到流域不同網格的土壤理化性質參數，而參與率定的土壤參數[22]共6個，其中控制研究流域地下基流量的土壤參數變量包括：Dmax，Ds，Ws，表征流域不同區域土壤蓄水能力大小的形狀變量b以及由VIC模型所劃分出來的土壤第二層厚度D2和第三層土壤厚度D3。本研究借助Matlab科學計算語言將待率定的VIC模型土壤參數設置成土壤參數變量，每次用不同的土壤參數變量的運行VIC模型，得到對應的模型模擬日值數據，將該次模擬數據與水文站實測徑流數據進行對比，對比完成后，更改土壤參數變量的值，再次運行VIC模型，重復上述過程，直到所有的土壤變量組合方式都運行完畢。

根據現有的武宣水文站實測日值徑流數據以及模型模擬日值徑流數據通過評價指標納什效率系數(Nash-Sutcliffe，NS)[23]、相關性系數(r)[24]以及相對偏差(BIAS)[25]對VIC模型的在對應土壤參數變量條件下的率定結果進行評估，評價指標具體數學公式見表2。

表2 模擬效果評價指標及其公式

由于西江流域面積廣闊其數據量相應的也會更大，導致每次模擬過程所需時間較長(30 min以上)，若將每種可能的組合方式均運行一次，在時間成本上是不可接受的，因此本研究參考Rosenbrock法[26]通過確定參數變化步長的方法來減少模型運行所需時間成本，即對VIC模型土壤參數變量的變化范圍進行劃分，每次只使用一個變化步長范圍內的土壤參數變量，輸入到VIC模型中，運行VIC模型，得到一組模擬徑流數據，重復上述過程直到所有的土壤參數變量組合都運行完完畢，最后根據評價指標的最優值找出最優土壤參數變量組合。

本研究將整個VIC模型模擬過程劃分為率定期(1981—1986年)和驗證期(1987—1989年)兩個階段。土壤參數變量的物理意義及率定工作完成后得到的土壤參數變量最優值結果見表3。

表3 VIC模型土壤參數變量物理意義及其最優值

3.2 模型適用性評估

以武宣站站率定期(1981—1986年)和驗證期(1987—1989年)實測日值徑流作為參考，VIC模型模擬結果見圖1—2。

圖1 武宣站率定期(1981-1986年)實測與模擬日流量過程對比

通過圖1率定期與圖2驗證期實測與模擬日流量過程對比以及表4西江流域VIC模型的評價指標在率定期和驗證期的結果可以看出在率定期：NS系數在日尺度上的值為0.71，在月尺度上的值為0.84。在相關性系數r的值上，日尺度的值為0.84，與NS系數的表現類似，月尺度的值較日尺度有所上升為0.95，表現優異。在相對偏差BIAS上，月尺度較日尺度也有所下降，分別為9%以及8%。在驗證期內：3個評價指標總體上較率定期均有不同程度的下降，其中NS系數在日尺度和月尺度上的值分別為0.65，0.85。相關性系數r的變化情況與NS系數的變化情況類似，在月尺度上的表現要優于日尺度，分別為0.81，0.94，在相對偏差上，日尺度的值為13%而月尺度的值為9%較率定期內的表現均有不同程度的增大。

表4 西江流域模擬結果

圖2 武宣站驗證期(1987—1989年)實測與模擬日流量過程對比

為了進一步分析VIC模型在西江流域年內豐水期(4—9月)和枯水期(10—12月、1—3月)的模擬效果是否存在差別，本研究將整個模擬歷時進行拆分，將每年的豐水期以及枯水期分別提取出來單獨組成一個徑流時間序列，分析VIC模型在年內不同時期的模擬效果是否存在差別。

圖3和圖4為率定期豐水期以及枯水期的水文過程線，由圖3可以看出在1981—1986年的豐水期中雖然水文過程曲線的變化趨勢保持一致，但在多個年份出現了洪水峰值模擬不足。但VIC模型在枯水期(圖4)的表現上，其模擬效果無論是從變化趨勢上看，還是對于峰值的模擬，其效果都要好于豐水期的模擬效果，部分情況下實測徑流量曲線甚至與模擬徑流量曲線處于重合狀態。

圖3 率定期西江流域豐水期水文過程線

圖4 率定期西江流域枯水期水文過程線

由表5中可知：日尺度下，率定期VIC模型在豐水期的NS的值為0.60，BIAS的值為13%，月尺度NS的值為0.63，BIAS的值為11%，而枯水期的NS的值為0.89，BIAS的值為4%。在月尺度下，枯水期的NS的值為0.91，BIAS的值略差于日尺度為6%。驗證期的模擬結果與率定期類似，具體為在日尺度驗證期的模擬結果中；豐水期NS的值為0.55，略優于率定期日尺度NS的值，BIAS的值為17%。在月尺度上NS的值為0.58，BIAS的值為13%。在枯水期的模擬中，日尺度上NS的值為0.78，BIAS的值為6%。月尺度上NS的值為0.87，BIAS的值為8%。因此在不同時間尺度條件下，VIC模型在枯水期的模擬效果無論是在率定期還是驗證期都要優于豐水期的模擬效果。

表5 西江流域模擬結果

由率定期和驗證期的水文過程線可以看出，VIC模型在年內不同時期的模擬精度上枯水期的模擬精度都要優于豐水期的模擬精度，但從總體上年際模擬效果上看，對比廖衛紅等人利用EasyDHM模型對西江流域的徑流模擬結果，雖然VIC模型的輸入數據種類更多、數據要求也更高，但該模型也有著更好的模擬效果，同時也能夠對流域熱量的傳輸過程進行模擬，這是EasyDHM模型所不具備的。綜合來說VIC模型在西江流域表現良好，可以為西江流域水資源合理分配、旱澇災害預警與防治提供一定的科學依據。

3.3 模型應用

利用經過率定與驗證得出的最優土壤參數組合，對VIC模型進行更深層次的應用。基于已經建立好的西江流域VIC模型框架，單獨制作了1971—2010年共四十年的氣象數據作為VIC模型的氣象驅動數據以驅動VIC模型運行，獲取1971—2010年的VIC模型產流過程輸出的西江流域多個水文要素數據。

VIC模型產流過程的輸出文件參數眾多，內容涵蓋了經過率定過程得出的最優水量性質參數以及最優能量性質參數。本研究提取了水量性質參數中的P(多年平均降水)、R(地表徑流)、E(蒸散發)、B(地下基流)4個水文要素，以3個月為一個時間步長，結合M-K趨勢檢驗方法基于網格探究西江流域P，R，E，B這4個水文要素的近40 a的變化規律。

本文以基于網格的形式，對VIC模型所劃分出的西江流域723個網格進行逐網格變化趨勢分析，其原理是根據VIC模型產流過程所生成西江流域723個網格1971—2010年共40 a的水文要素，每一個網格單獨進行M-K變化趨勢檢驗分析，計算各個網格的Z值，在所有網格的完成M-K變化趨勢分析后，將所有的Z值進行匯總，制作成ASCII格式，導入ArcMap中進行可視化處理。

研究表明：在這40 a中降雨在西江流域約一半的地區有增加趨勢，而Z>1.64的網格共213個占西江流域總網格的29%，主要分布在西江流域的西部云南和貴州境內以及東部廣西境內，降雨表現出下降趨勢的網格主要集中在西江流域中部，但其Z均沒有大于1.64，即增加趨勢不明顯。蒸散發具有增加趨勢的網格集中在西江流域的中部貴州境內，而Z>1.64的網格較降水有一定的增加，共266個網格，占流域總網格數的37%；與降雨類似，蒸散發也具有部分網格具有下降趨勢，但其Z值均沒有大于1.64。西江流域地表經流要素具有增加趨勢的網格最多，其Z>1.64的網格共476個，在流域總網格數中占比超過65%，主要集中在廣西境內流域部分。地下基流相較于其他3個水文要素變化最不明顯，具有顯著增加趨勢的網格共10個，主要集中在廣西境內，占比不到2%。

4 結 論

流域自然水循環過程的模擬需要綜合考慮多方面的因素，尤其是大尺度的流域水文物理過程尤其復雜。本研究依托RS技術和GIS空間分析技術完成了分布式水文模型VIC模型的資料收集和數據預處理工作，在西江流域構建了VIC模型數據庫，利用3個評價指標(NS，r，BIAS)評估VIC模型在西江流域徑的適用性，結果表明：VIC模型在不同的時間尺度的徑流模擬過程中均表現出了良好的適用性。文章最后對VIC模型的輸出文件做了進一步的應用，從中提取了P，E，R，B共4個水文要素，利用M-K變化趨勢分析方法，對這西江流域這4個水文要素在1971—2010年共四十年內的變化趨勢進行了逐網格分析，分析結果表明：西江流域有超過65%的流域面積的地表徑流R具有顯著增加趨勢，年均蒸散發E和年均降雨P具有顯著增加趨勢的面積占流域總面積約35%。地下基流B則相對穩定，具有顯著變化趨勢的面積占流域總面積不到2%。

盡管大尺度VIC模型在流域自然水文物理過程的模擬研究中具有眾多優點，但是由于研究區數據資料的限制以及模型參數的不確定性等方面的影響，會對實際的模擬過程產生不良后果。如本研究在西江流域所采用的VIC模型分辨率為2 km，較大的模型分辨率會實測流域內很多的細節被忽略、概化，從而影響模型的模擬精度，因此在未來的研究中可以考慮采用更加精細的數據精度以提高模型的模擬精度從而更好的了解流域自然水循環過程的變化特點。