VIC模型系統(tǒng)微分響應(yīng)參數(shù)率定方法在大坡嶺流域的應(yīng)用

孫苗苗，石 朋，瞿思敏，陳學(xué)秋，徐 瑤，游 洋，包沐曦，宋蘭蘭

（1.河海大學(xué)水文水資源學(xué)院，南京 210098；2.水利部珠江水利委員會(huì)水文局，廣州 510611；3.江蘇省水利科教中心，南京 210003；4.江蘇省水利廳，南京 210029；5.南京海事局，南京 210009）

0 引 言

流域水文模型是研究降雨徑流形成機(jī)理及時(shí)空變化規(guī)律的有效工具［1］。分布式水文模型能考慮流域特征空間變異性且具有物理基礎(chǔ)，在水文模擬中具有優(yōu)越性，是當(dāng)前流域水文模擬的重要發(fā)展趨勢(shì)之一［2］。在眾多分布式水文模型中，VIC（Variable Infiltration Capacity）模型因其操作原理簡(jiǎn)單、模擬精度較好的特點(diǎn)，被廣泛成功應(yīng)用于國(guó)內(nèi)外不同尺度的流域，研究區(qū)域覆蓋廣，研究?jī)?nèi)容涉及豐富，包括徑流模擬、模型改進(jìn)及氣候和土地覆被變化對(duì)水資源影響等諸多方面［3-9］。

水文模型的模擬效果很大程度上依賴于模型參數(shù)的選擇，因此，選擇合理、高效的參數(shù)率定方法極其重要［10］。目前，大多數(shù)的參數(shù)率定方法采用誤差平方和目標(biāo)函數(shù)來(lái)進(jìn)行評(píng)估，但誤差平方和的累加過(guò)程損失了有效信息，其提供的信息存在諸多不合理甚至是錯(cuò)誤的。近年來(lái)，河海大學(xué)的包為民教授提出了系統(tǒng)微分響應(yīng)參數(shù)率定方法［11，12］，該方法能有效避免現(xiàn)有參數(shù)率定方法的絕大多數(shù)問(wèn)題，并在新安江模型、SAC 模型等的參數(shù)優(yōu)化中均取得了較好的應(yīng)用效果［10，13］，但該方法在大尺度分布式水文模型上的應(yīng)用研究相對(duì)較少?；谏鲜鲈?，本文選取淮河上游大坡嶺子流域?yàn)槔趯?duì)模型參數(shù)敏感性分析的基礎(chǔ)上，將系統(tǒng)微分響應(yīng)參數(shù)率定方法與VIC模型相結(jié)合，探究驗(yàn)證其在大尺度分布式水文模型中的適用性。

1 研究方法

1.1 VIC模型

VIC 模型［14］基于SVATS（Soil Vegetation Atmospheric Transfer Schemes）思想的空間分布網(wǎng)格化的大尺度分布式水文模型。模型主要通過(guò)考慮大氣、植被、土壤之間的物理化學(xué)生物等轉(zhuǎn)化過(guò)程來(lái)反映陸面各部分間的水熱狀態(tài)變化和傳輸狀況，可同時(shí)模擬水循環(huán)過(guò)程的水量平衡和能量平衡，在實(shí)際應(yīng)用中，也可只計(jì)算水量平衡［15］。

VIC 模型規(guī)定對(duì)研究區(qū)域網(wǎng)格化，在每個(gè)網(wǎng)格內(nèi)分別考慮了裸地及多種植被覆蓋類型。在蒸散發(fā)過(guò)程中，同時(shí)考慮了植被冠層蒸發(fā)、植被蒸騰及裸土蒸發(fā)。在產(chǎn)流計(jì)算方面，VIC-3L模型將土壤劃分為三層，頂薄層和上層土壤產(chǎn)生直接徑流，下層土壤產(chǎn)生基流［15］，其中，直接徑流的計(jì)算考慮了蓄滿產(chǎn)流和超滲產(chǎn)流兩種機(jī)制，基流則根據(jù)Arno 概念模型［16］計(jì)算。VIC 模型無(wú)匯流模型，本文采用了Dag Lohmann 等人開(kāi)發(fā)的匯流模型與之耦合。

1.2 系統(tǒng)微分響應(yīng)參數(shù)率定方法

系統(tǒng)微分響應(yīng)參數(shù)率定方法（System Response Parameter Calibration Method，SRPCM）是一種基于參數(shù)函數(shù)曲面的參數(shù)估計(jì)尋優(yōu)方法［17，18］。不同于傳統(tǒng)優(yōu)化在目標(biāo)函數(shù)曲面尋找參數(shù)優(yōu)值的想法，該方法在參數(shù)函數(shù)曲面上尋優(yōu)，成功解決了目標(biāo)函數(shù)帶來(lái)的問(wèn)題［10］。

水文模型常可看作為沒(méi)有明確表達(dá)式的函數(shù)Yt=f(P，Xt)，其中，模型的輸入是自變量Xt，模型的輸出為因變量Yt，模型參數(shù)為函數(shù)參數(shù)P，當(dāng)Xt確定不變后，Yt隨P而改變，這便形成了定義在參數(shù)空間域上的函數(shù)曲面。

給定一個(gè)實(shí)測(cè)樣本(Xt，Qt)，將自變量代入函數(shù)或模型可得：

式中：Xt取實(shí)測(cè)值；Yt隨參數(shù)P而改變，把該式定義的函數(shù)稱為參數(shù)函數(shù)。

將實(shí)測(cè)樣本(Xt，Qt)代入函數(shù)或模型得：

上式是與特定實(shí)測(cè)樣本對(duì)應(yīng)的參數(shù)域函數(shù)，將式（2）定義的函數(shù)稱為樣本函數(shù)。

利用系統(tǒng)微分理論，對(duì)于一個(gè)多參數(shù)函數(shù)Yt=f(Pi，Xt)，對(duì)自變量{P1，P2，…，Pn}進(jìn)行微分可得：

式中：函數(shù)增量dYt隨參數(shù)增量{dP1，dP2，…，dPn}而改變，稱為參數(shù)增量的全微分系統(tǒng)響應(yīng)。當(dāng)只改變某一參數(shù)Pi其他參數(shù)不變時(shí)，函數(shù)增量dYt是參數(shù)增量dPi的偏微分系統(tǒng)響應(yīng)。

因水文模型常常沒(méi)有明確數(shù)學(xué)表達(dá)式，無(wú)法進(jìn)行偏微分的計(jì)算，因此可用差分格式。對(duì)式（3）中向前差分可得：

式中：Pji表示參數(shù)Pi的第j次率定；Pji+1=Pi+ ΔPi表示參數(shù)Pi的第j+ 1 次率定；其中ΔPi為參數(shù)Pi增加的一個(gè)微小值。將f ij表示為代入Pji的函數(shù)值，f ij+1表示為代入Pji+1的函數(shù)值，則式（4）可表示為：

則考慮所有參數(shù)的差分表達(dá)式為：

假設(shè)有N組實(shí)測(cè)樣本(X1，Q1)，(X2，Q2)，…，(XN，QN)，將其代入式（6），可得：

用向量形式表示式（7）如下：

式中：e為實(shí)測(cè)值Q與函數(shù)計(jì)算值f j+1的偏差；s是參數(shù)靈敏度矩陣，表達(dá)式如下：

對(duì)新的估計(jì)參數(shù)Pj+1應(yīng)使式（8）的誤差平方和（SLS）最小，即minP?Rn{SLSj+1=(Q-f j+1)T(Q-f j+1)}。

當(dāng)minP?RnSLSj+1= 0時(shí)，結(jié)合式（8）可得Pj+1的最佳估計(jì)：

系統(tǒng)微分響應(yīng)參數(shù)率定方法的具體計(jì)算步驟如下：

（1）設(shè)置參數(shù)初值P0，輸入N組實(shí)測(cè)值(X1，Q1)，(X2，Q2)，…，(XN，QN)，令j= 0，Pj=P0；

（2）將模型輸入X1，X2，…，XN和參數(shù)Pj代入模型計(jì)算得到一組輸出f j=(f1，f2，…，fN)，并根據(jù)式（9）計(jì)算參數(shù)靈敏度矩陣s；

（3）根據(jù)式（10）確定新參數(shù)Pj+1；

（4）判斷參數(shù)及模擬結(jié)果是否滿足收斂條件，若是則退出尋優(yōu)，否則令j=j+ 1，跳轉(zhuǎn)到第（2）步繼續(xù)循環(huán)。

2 區(qū)域概況及數(shù)據(jù)

2.1 研究區(qū)域概況

淮河流域位于我國(guó)東部，介于黃河、長(zhǎng)江兩流域之間，流域經(jīng)緯度范圍分別為：111°55′～121°25′E，30°55′～36°36′N。本文選取的研究區(qū)域是位于淮河干流上游的大坡嶺水文站以上流域（圖1），該站控制流域面積為1 640 km2，干流長(zhǎng)度為73 km。大坡嶺流域?qū)倥瘻貛О霛駶?rùn)季風(fēng)氣候，氣候溫和，流域多年平均降雨量為918 mm。流域內(nèi)地形以山區(qū)和丘陵區(qū)為主，水利工程不多。流域地處亞熱帶常綠闊葉林帶，植被資源豐富，主要的土地利用類型為農(nóng)林業(yè)，主要作物為水稻［19］。

2.2 研究區(qū)數(shù)據(jù)

構(gòu)建VIC模型需要的輸入數(shù)據(jù)主要分為以下4類：

（1）地形數(shù)據(jù)。采用了SRTMDEM（90M）分辨率原始高程數(shù)據(jù)，使用ARCGIS 將研究區(qū)劃分成61 個(gè)0.05°×0.05°分辨率的網(wǎng)格單元并生成水系等，見(jiàn)圖2（a）。

（2）氣象強(qiáng)迫數(shù)據(jù)。來(lái)自中國(guó)地面氣候資料日值數(shù)據(jù)集（V3.0），選取了大坡嶺流域及臨近氣象站2007年1月1日-2014年12月31日的氣象資料，采用IDW 法插值獲取降雨量、日最高氣溫、日最低氣溫及風(fēng)速四個(gè)氣象要素在流域內(nèi)的分布。

（3）植被數(shù)據(jù)。植被分類采用的是馬里蘭大學(xué)（UMD）發(fā)布的全球1 km 分辨率土地覆蓋類型數(shù)據(jù)集，見(jiàn)圖2（b），大坡嶺流域內(nèi)草原和耕地所占的面積比例較大，說(shuō)明流域內(nèi)以人工種植的植被為主，各植被類型對(duì)應(yīng)的特定參數(shù)主要依據(jù)LDAS 提供的值來(lái)確定。

（4）土壤數(shù)據(jù)。選取了1 km分辨率的基于HWSD的中國(guó)土壤數(shù)據(jù)集及FAO 世界糧農(nóng)組織的5′分辨率的土壤質(zhì)地分布數(shù)據(jù)集，VIC模型上（下）層土壤類型根據(jù)HWSD數(shù)據(jù)庫(kù)確定，如圖2（c）、（d）所示，可見(jiàn)大坡嶺流域上層土壤以黏壤土和砂質(zhì)黏壤土為主，下層土壤以黏壤土為主。

3 流域應(yīng)用

3.1 參數(shù)敏感性分析

在制備土壤參數(shù)時(shí)發(fā)現(xiàn)，VIC模型的部分水文參數(shù)（表1）依托于概念性產(chǎn)流模型，缺乏實(shí)測(cè)數(shù)據(jù)，需要通過(guò)經(jīng)驗(yàn)估計(jì)或參數(shù)率定來(lái)確定取值。在眾多有關(guān)VIC 模型的研究應(yīng)用中，對(duì)模型參數(shù)的敏感性分析研究相對(duì)很少［20］。鑒于不同參數(shù)對(duì)模擬效果有不同程度的影響，探究模型參數(shù)的敏感性，僅選擇敏感參數(shù)進(jìn)行率定，對(duì)提高模型應(yīng)用效率具有重要意義。

表1 部分土壤參數(shù)及其取值范圍Tab.1 Some soil parameters and their value ranges

本文使用LH-OAT 敏感性分析方法對(duì)這8個(gè)參數(shù)進(jìn)行敏感性分析，設(shè)置m=5（即將參數(shù)的空間分布等分為5 層），fi= 0.1（即每個(gè)參數(shù)的微小變化設(shè)置為0.1）。選用納什系數(shù)（NSE）、相對(duì)誤差（RE）、高水量誤差系數(shù)（fg）及低水量誤差系數(shù)（fd）共4個(gè)評(píng)估指標(biāo)來(lái)評(píng)估各個(gè)參數(shù)對(duì)日徑流過(guò)程的敏感性。4個(gè)評(píng)估指標(biāo)的計(jì)算公式如下：

式中：Qc(i)為計(jì)算流量；Qob(i)為實(shí)際流量；N為模擬的序列長(zhǎng)度。

本文針對(duì)不同參數(shù)，對(duì)其敏感性按從大到小進(jìn)行了排序（表2），并將敏感度劃分等級(jí)，分別是：排序?yàn)? 的是極敏感；排序?yàn)?～7 的是較敏感；排序?yàn)? 的是不敏感［21］。最后綜合統(tǒng)計(jì)各參數(shù)在不同指標(biāo)下的最小排序值作為其全局敏感性的排序結(jié)果。

表2 VIC模型參數(shù)敏感性分析結(jié)果Tab.2 Parameters sensitivity results of the VIC model

由表2可以看出，參數(shù)d1（上層土壤厚度）在4 個(gè)指標(biāo)下均為最高敏感度，即d1是最敏感的參數(shù)。在VIC 模型中，上層土壤厚度直接影響上層土壤含水量與裸土蒸發(fā)，從而影響直接徑流的產(chǎn)生，因而d1對(duì)模型的敏感性較高。

對(duì)于關(guān)注高流量值的評(píng)估指標(biāo)fg來(lái)說(shuō)，蓄水容量曲線參數(shù)B 的敏感性較高，為較敏感參數(shù)。參數(shù)B反映了流域下滲能力空間分布的不均勻性，對(duì)產(chǎn)流過(guò)程影響巨大。大坡嶺流域面積雖不大，下墊面分布不均勻性卻較為明顯，所以B對(duì)產(chǎn)流較敏感。

參數(shù)Ds（非線性基流發(fā)生時(shí)占最大基流比例）、Dm（最大基流量）、Ws（非線性基流發(fā)生時(shí)占最大土壤含水量比例）、d0（頂薄層厚度）和d2（下層土壤厚度）對(duì)產(chǎn)流起間接影響，全局敏感性均為較敏感。其中，參數(shù)Ds、Dm和Ws與模型基流計(jì)算有關(guān)，Dm通過(guò)影響基流量大小從而影響徑流，Ds和Ws與非線性基流有關(guān)，這3個(gè)參數(shù)在低水量徑流計(jì)算中敏感性較高，說(shuō)明結(jié)果與實(shí)際水文過(guò)程較為相符。d0影響頂薄層濕度和裸土蒸發(fā)，一般在模型運(yùn)算時(shí)設(shè)置為缺省值0.1 m，但從結(jié)果來(lái)看，其對(duì)徑流模擬效果有一定的影響。

參數(shù)C在4 種指標(biāo)下均表現(xiàn)為不敏感，說(shuō)明該參數(shù)對(duì)模型模擬結(jié)果的影響不大，在模型運(yùn)行開(kāi)始時(shí)設(shè)置經(jīng)驗(yàn)值即可。

本研究的結(jié)果與張續(xù)軍［15］和宋星原［22］的結(jié)果較為接近。張續(xù)軍認(rèn)為d1對(duì)多年徑流相對(duì)誤差較為敏感，B 對(duì)日徑流過(guò)程擬合程度較為敏感。宋星原認(rèn)為d1、Ws和B 對(duì)徑流相對(duì)誤差敏感度相對(duì)較高，對(duì)于納什系數(shù)，d1的敏感度遠(yuǎn)高于其他參數(shù)。

綜上，VIC 模型土壤參數(shù)中的B、Ds、Dm、Ws、d0、d1及d2這7個(gè)參數(shù)較為敏感，在后續(xù)研究中應(yīng)注重考慮率定優(yōu)化，參數(shù)C為不敏感參數(shù)，本文便不對(duì)其進(jìn)行率定，設(shè)置為固定值2。

3.2 理想模型應(yīng)用

為檢驗(yàn)系統(tǒng)微分響應(yīng)參數(shù)率定方法（SRPCM）在VIC 模型中的適用性，構(gòu)建VIC 理想模型進(jìn)行單參數(shù)率定檢驗(yàn)。本文隨機(jī)設(shè)置一套參數(shù)真值（表3），將實(shí)測(cè)輸入數(shù)據(jù)和參數(shù)真值代入到VIC模型計(jì)算得到的流量過(guò)程即為“實(shí)測(cè)徑流”。文中通過(guò)率定單參數(shù)來(lái)驗(yàn)證SRPCM 的應(yīng)用效果，根據(jù)上節(jié)敏感性分析結(jié)果，選擇土壤參數(shù)中B、Ds、Dm、Ws、d1及d2這6 個(gè)參數(shù)進(jìn)行參數(shù)優(yōu)化（因考慮到模型將d0設(shè)置了缺省值0.1 m，且在后續(xù)運(yùn)行時(shí)發(fā)現(xiàn)當(dāng)d0＞d1時(shí)會(huì)終止計(jì)算，故本文將此值設(shè)置成默認(rèn)值，不再進(jìn)行率定）。

表3 理想模型參數(shù)真值Tab.3 The real value of ideal model parameters

圖3為分別對(duì)6 個(gè)參數(shù)率定時(shí)，各參數(shù)初始值和最優(yōu)值的分布，其中，平均迭代次數(shù)最多的為參數(shù)d2，共迭代4次，單組迭代次數(shù)最多的是d2第五組，共迭代7次。從圖3可以看出，對(duì)于隨機(jī)給定的不同初值，各參數(shù)經(jīng)過(guò)優(yōu)化均能收斂到表3的真值?？梢?jiàn)，SRPCM 在VIC 理想模型單參數(shù)率定中具有一定的適用性，且具有高效、穩(wěn)定、結(jié)果合理的特性，說(shuō)明將此方法應(yīng)用于VIC模型是可行的。

3.3 實(shí)際流域應(yīng)用

參數(shù)率定及驗(yàn)證數(shù)據(jù)資料選用了大坡嶺流域2006-2016年共11年的日徑流數(shù)據(jù)。根據(jù)VIC模型參數(shù)率定原則，選擇2006年為預(yù)熱期，模擬結(jié)果不計(jì)入目標(biāo)函數(shù)，2007-2014年為率定期，2015、2016年為驗(yàn)證期。實(shí)測(cè)流量數(shù)據(jù)為大坡嶺水文站日流量過(guò)程。

本文使用系統(tǒng)微分響應(yīng)參數(shù)率定方法進(jìn)行參數(shù)率定，為防止迭代過(guò)程陷入死循環(huán)，設(shè)置最大迭代次數(shù)為500次，以納什系數(shù)（NSE）作為目標(biāo)函數(shù)，率定結(jié)果見(jiàn)表4。

表4 大坡嶺流域參數(shù)率定結(jié)果Tab.4 The parameter calibration results of Dapoling catchment

本文僅進(jìn)行了水量平衡模擬，率定期和驗(yàn)證期的納什系數(shù)（NSE）和徑流總量相對(duì)誤差（RE）見(jiàn)表5，從結(jié)果可以看出，VIC模型對(duì)大坡嶺流域的模擬效果較好，率定期和驗(yàn)證期的徑流總量相對(duì)誤差分別為4.64%、7.09%，均控制在8%以內(nèi)，徑流過(guò)程的納什系數(shù)分別為0.600、0.601，均達(dá)到了0.6，驗(yàn)證期的NSE略高，據(jù)此可以看出采用系統(tǒng)微分響應(yīng)參數(shù)率定方法優(yōu)選的VIC模型參數(shù)，對(duì)于大坡嶺流域的日徑流過(guò)程模擬具有較好的適用性。

表5 大坡嶺流域徑流模擬效果Tab.5 Runoff simulation results of Dapoling catchment

圖4點(diǎn)繪了率定期與驗(yàn)證期日流量模擬值與實(shí)測(cè)值的散點(diǎn)圖，可以看出在率定期和驗(yàn)證期，絕大部分散點(diǎn)比較均勻分布在對(duì)稱線的兩側(cè)，但在高流量段，不論是率定期還是驗(yàn)證期，都有幾個(gè)散點(diǎn)位于偏右下的位置，說(shuō)明在日徑流模擬過(guò)程中，VIC模型對(duì)于流量峰值段的模擬結(jié)果偏小。率定期和驗(yàn)證期的相關(guān)系數(shù)分別為0.790 7、0.811 7，線性回歸系數(shù)R2分別為0.625 3、0.658 9，相較于率定期，驗(yàn)證期的模擬效果均有所提高，這表明采用系統(tǒng)微分參數(shù)率定方法得到的最優(yōu)參數(shù)值具有一定可信度，該方法適用于VIC模型參數(shù)率定。

對(duì)比率定期及驗(yàn)證期流量過(guò)程線（圖5、圖6）可知，模擬與實(shí)測(cè)徑流過(guò)程總體擬合度良好，特別是洪峰出現(xiàn)時(shí)間和退水時(shí)間的模擬基本與實(shí)測(cè)徑流過(guò)程相一致。但對(duì)于部分較大的流量峰值（特別是2008年、2014年及2015年等日徑流峰值段），模擬結(jié)果不是特別理想，在流量峰值時(shí)模擬值基本偏小，這與上文圖4散點(diǎn)圖得出的結(jié)論相一致。分析其原因，可能與模型本身參數(shù)和結(jié)構(gòu)的局限性、研究區(qū)域偏小以及收集的資料缺乏代表性等因素有關(guān)。

4 結(jié) 語(yǔ)

本文采用LH-OAT 方法對(duì)VIC 模型參數(shù)進(jìn)行敏感性分析，將系統(tǒng)微分響應(yīng)參數(shù)率定方法與VIC模型結(jié)合并應(yīng)用于大坡嶺流域，得出以下結(jié)論：①VIC模型的敏感參數(shù)以上層土壤厚度d1為主，模型參數(shù)的敏感性分析，對(duì)后續(xù)參數(shù)率定具有一定的參考價(jià)值；②理想模型和實(shí)際流域的應(yīng)用檢驗(yàn)表明，用系統(tǒng)微分響應(yīng)參數(shù)率定方法率定VIC模型參數(shù)是可行的，率定結(jié)果合理、高效；③將該方法率定VIC 模型敏感參數(shù)的結(jié)果應(yīng)用于大坡嶺流域，取得了較好的應(yīng)用效果，但對(duì)徑流過(guò)程中的較大洪峰模擬效果不是十分理想，可能與模型本身的局限性等原因有關(guān)，值得未來(lái)進(jìn)一步的研究與改進(jìn)予以提高水文模擬精度?！?/p>