系統(tǒng)微分響應(yīng)誤差修正方法在VIC 模型中的應(yīng)用與研究

蔣語珣，瞿思敏，蔣思軍，嵇海祥，司 偉，石 朋，萬 皓

（1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2. 浙江省臺州市長潭水庫事務(wù)中心，浙江 臺州 318024；3. 南京水利自動化研究所 江蘇 南京 210008； 4. 水利部南京水利水文自動化研究所 江蘇 南京 210012）

0 引 言

水文模型是探索和認(rèn)識流域水循環(huán)和水文過程的重要手段，在進(jìn)行水文規(guī)律研究、水文預(yù)報、水資源規(guī)劃與管理、水文分析與計算和其他的生產(chǎn)實際問題中起著重要的作用［1，2］。水文模型按照對流域水文過程描述的離散程度分為集總式模型、分布式模型和半分布式模型，其中分布式模型是根據(jù)流域各處的氣候信息和下墊面特性的差異性，劃分流域為多個小單元，充分考慮了流域內(nèi)部的各項氣候地理條件的差異性。在考慮流域各處氣候信息和下墊面特性的不同的基礎(chǔ)上，將流域劃分成多個基本單元，具有從機(jī)理上考慮降雨和下墊面條件空間分布不均勻?qū)α饔蚪涤陱搅餍纬捎绊懙墓δ埽?］。VIC（Variable Infiltration Capacity）模型作為分布式模型代表，已被廣泛應(yīng)用于不同尺度的流域水文模擬中，并都取得了不錯的模擬效果，說明VIC模型在流域水文模擬中適用性較好［3-5］。

水文預(yù)報誤差因素很多，包括輸入數(shù)據(jù)的誤差，模型結(jié)構(gòu)的誤差和模型參數(shù)的誤差等，因此，對預(yù)報結(jié)果進(jìn)行誤差修正就顯得尤為重要。常用的誤差修正方法有：自回歸模型（AR 模型）［6］、卡爾曼濾波修正技術(shù)［7］、抗差修正技術(shù)［8］和基于人工經(jīng)驗的修正方法［9］等。這些修正方法大多數(shù)是基于流量誤差的統(tǒng)計分析，并未將實測流量中所包含的信息進(jìn)行分類提取利用，導(dǎo)致研究要素的誤差與模型本身分離，修正效果有時并不理想。河海大學(xué)包為民教授在2015 年提出系統(tǒng)微分響應(yīng)理論［9-12］，基于該理論的修正方法就是利用自變量誤差與計算結(jié)果誤差之間的關(guān)系，根據(jù)模型計算誤差量來估計自變量或影響因素誤差量，再據(jù)估計誤差量來修正模型計算結(jié)果，此方法結(jié)構(gòu)簡單，且不增加模型參數(shù)。并已在很多流域的水文模型［13］中得到了驗證。在VIC 模型中，系統(tǒng)微分響應(yīng)理論僅用于模型的參數(shù)率定，尚未用于模型的誤差修正中。因此，本文以閩江建陽流域為研究區(qū)域，構(gòu)建區(qū)域VIC 模型，分析模擬結(jié)果的誤差，應(yīng)用系統(tǒng)微分響應(yīng)誤差修正方法通過對降雨量進(jìn)行誤差修正來修正模擬結(jié)果，并與AR 模型修正方法進(jìn)行對比分析，探究系統(tǒng)響應(yīng)誤差修正方法在VIC模型中的適用效果。

1 研究方法

1.1 VIC模型

VIC（Variable Infiltration Capacity）模型，也可稱為“可變下滲容量模型”，屬于流域大尺度分布式水文模型。模型易于氣候模式的嵌套，且具有很強(qiáng)的物理基礎(chǔ)，能靈活地處理各種復(fù)雜的應(yīng)用條件。VIC模型主要通過考慮大氣、植被、土壤之間的物理化學(xué)生物等轉(zhuǎn)化過程來反映陸面各部分之間的水熱動態(tài)變化及傳輸狀況，可同時開啟水量平衡和能量平衡模式，也可只進(jìn)行水量平衡的計算［14］。

VIC模型在水平方向上將研究區(qū)域劃分為大小相同的正方形網(wǎng)格，在各個網(wǎng)格上遵從水量平衡和能量平衡原理模擬產(chǎn)匯流過程。在模型初建時將土壤分為上下兩層，后改進(jìn)為三層，增加頂薄層描述表層土壤水的動態(tài)變化，能顯著提高模擬表層濕度的精度［15］。

VIC模型蒸散發(fā)計算包括冠層蒸發(fā)、植被蒸騰和裸地蒸發(fā)，每個部分單獨計算后按照面積權(quán)重相加即可得到單個網(wǎng)格的總蒸發(fā)蒸騰量。

產(chǎn)流計算中，考慮次網(wǎng)格尺度的土壤空間分布不均勻性，將超滲產(chǎn)流和蓄滿產(chǎn)流機(jī)制相結(jié)合［16，17］，運用到每個計算網(wǎng)格內(nèi)，模擬計算地表徑流量。上層土壤（包括第一層和第二層）產(chǎn)生直接徑流Qd，計算公式如式（1），下層土壤（包括第三層）產(chǎn)生基流Qb，計算公式如式（2）。

式中：Δt為計算時間；P為降雨量；W c1為上層土壤的最大土壤含水量；W-1為上層土壤初始土壤含水量；i0為初始下滲能力；im為最大下滲能力；b為可變下滲能力曲線形狀參數(shù)。

式中：Dm為最大基流；Ds為非線性基流出現(xiàn)時占Dm的比例；W-2為下層土壤初始土壤含水量；W c2為下層土壤的最大土壤含水量；Ws為非線性基流出現(xiàn)時占W c2的比例，且有Ds≤Ws。

匯流模型采用“先演后合”的線性匯流方案，河道匯流采用線性圣維南方程，坡面匯流采用單位線法［18］。

1.2 系統(tǒng)微分響應(yīng)誤差修正方法

將模型計算概化成一個系統(tǒng)，將出口斷面流量誤差作為信息來源，計算時段降雨引起的模型微分響應(yīng)來反演計算降雨量誤差，再修正流域出口斷面的流量過程。圖1 中方框內(nèi)G（t）為系統(tǒng)輸入項，Q（t）為系統(tǒng)輸出項，VIC 表示模型計算，RDR 表示降雨微分響應(yīng)，CP表示由RDR計算出的降雨量修正量。

圖1 誤差修正方法示意圖Fig.1 Schematic diagram of error correction method

根據(jù)上述系統(tǒng)計算過程，可用以下公式來表示：

式中：in(t)為模型的輸入變量，如降雨，溫度，風(fēng)速等；m(t)表示模型的中間狀態(tài)變量；θ表示模型的參數(shù)；t表示時間。

本文只考慮由于模型輸入降雨量P的變化引起的系統(tǒng)微分響應(yīng)，因此，將式（3）簡化為：

假設(shè)樣本系列長度為L，式中P(t) =[P1，P2，P3，…，PL]T表示時段降雨量系列。本文采用微分線性化的思想將其線性化，建立降雨量增量與響應(yīng)函數(shù)f之間的微分響應(yīng)關(guān)系。

對式（4）兩邊同時求微分得系統(tǒng)輸出響應(yīng)過程，如式（5）所示：

式中：dQ(t)表示系統(tǒng)輸出增量響應(yīng)表示系統(tǒng)微分響應(yīng)函數(shù)；dP表示降雨量增量。

式（5）用矩陣表示為：

式中：CQ=[ΔQ1，ΔQ2，ΔQ3，…，ΔQL]T，ΔQi為i時刻實測與計算流量誤差；CP為求解的降雨修正量，假設(shè)需要修正的降雨量長度為h，Cp=[ΔP1，ΔP2，ΔP3，…，ΔPh]T；C為流量觀測隨機(jī)誤差項，一般認(rèn)為服從零均值分布，為白噪聲向量，C=[e1，e2，e3，…，eL]T。

S矩陣的表達(dá)式為：

式（7）中，S表示時段降雨量Pi所對應(yīng)模型的微分響應(yīng)，其每一項都可用以下差分近似求解：

式（8）中，t=1，…，L，當(dāng)t＜i時：

本文研究選擇一個單位的降雨量，即令式（8）中的ΔPi= 1，即C 矩陣就代表降雨量的單位變化量引起的模型降雨微分響應(yīng)。

降雨量的系統(tǒng)微分響應(yīng)修正步驟如下：

（1）計算實測流量與計算流量誤差CQ=QO-Qm，QO為實測流量，Qm為計算流量。

（2）以日為單位確定需要修正的降雨量時段數(shù)n，在每日的降雨量Pi上（i=1…n）依次增加一個單位的降雨量，其他時段降雨量不變，得到新的降雨量系列。

（3）用新的降雨量系列通過VIC模型計算后得到流量過程。

（4）用步驟（3）計算得到的流量過程減去用原降雨量系列P計算得到的流量過程線，即為降雨量Pi的微分響應(yīng)，表示為Si。

S矩陣中的每一列均用相同的方法求得。根據(jù)最小二乘估計原理，得到的降雨修正量的計算公式是：

修正后的降雨量系列為：

P'為修正后的降雨量，將其重新代入VIC模型中進(jìn)行計算，即可得到修正后的流域出口斷面的計算流量過程。

2 研究流域概況及數(shù)據(jù)來源

2.1 研究流域概況

選取閩江上游建陽站以上流域為研究區(qū)域。建陽站所在的河流屬于閩江上游建溪水系，流域范圍在東經(jīng)117°31'～118°19'，北緯27°10'～28°5' 之間，建陽站以上流域控制面積為4 848 km2，屬于亞熱帶季風(fēng)氣候，氣候溫和，年平均氣溫為18.1 ℃，多年平均降雨量1 741 mm，年日照平均1 802 h，年徑流系數(shù)為0.56。選取的研究區(qū)內(nèi)有邵武、浦城、武夷山、建陽、建甌、玉山、衢州、麗水、東溪、黃坑、黎源、興田、書坊共13 個氣象站和1個水文站（建陽水文站），水系及站點分布圖見圖2。

圖2 閩江建陽流域水系及站點分布圖Fig.2 Distribution map of water system and stations in Jianyang Watershed of Minjiang River

2.2 資料與數(shù)據(jù)

在建立VIC 模型時，需要將流域劃分成多個10 km×10 km分辨率的網(wǎng)格單元并生成水系，模型輸入的文件包括植被參數(shù)文件、土壤參數(shù)文件和氣象驅(qū)動數(shù)據(jù)。制作模型輸入文件所需的數(shù)據(jù)來源如表1。

表1 原始數(shù)據(jù)說明Tab.1 Raw data description

3 建陽流域VIC模型構(gòu)建與參數(shù)率定

3.1 建陽流域VIC模型構(gòu)建

對原始高程數(shù)據(jù)進(jìn)行剪裁和重采樣，將流域劃分成66 個10 km×10 km 分辨率的網(wǎng)格，按照從上到下、從左到右的順序依次為其編號，并生成閩江建陽流域水系，如圖3（a）。根據(jù)Maryland University 全球1 km 分辨率土地覆蓋類型數(shù)據(jù)集對閩江建陽流域的下墊面覆蓋情況進(jìn)行研究，如圖3（b），閩江建陽流域植被類型以常綠針葉林為主，占比達(dá)到55.89%，其次是林地，占比達(dá)到20.33%。制備的VIC 模型的頂薄層及上層土壤類型根據(jù)HWSD 數(shù)據(jù)庫上層土壤分布確定，下層土壤類型同樣根據(jù)HWSD 數(shù)據(jù)庫下層土壤分布確定，如圖3（c），（d），上層土壤分布以砂質(zhì)壤土（Sandy loam）和砂質(zhì)黏土（Sandy clay）為主，分別占比達(dá)到51.73%和16.47%；下層土壤分布以砂質(zhì)壤土（Sandy loam）和粉土（Silt）為主，分別占比達(dá)到51.73%和24.69%。本文氣象數(shù)據(jù)采用1989-2000年中國地面國際交換站氣候資料日值數(shù)據(jù)集，該數(shù)據(jù)集共包含氣象站756個，選取閩江建陽流域內(nèi)以及周邊13個氣象站日尺度數(shù)據(jù)作為原始數(shù)據(jù)（氣象站分布圖見圖2），將氣象站數(shù)據(jù)根據(jù)反距離加權(quán)法（IDW）插值到各網(wǎng)格的計算單元中心點，制備各網(wǎng)格的氣象驅(qū)動文件。

圖3 閩江建陽流域VIC模型構(gòu)建圖Fig.3 Building map of VIC model in Jianyang Watershed of Minjiang River

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

VIC 模型的部分土壤參數(shù)是根據(jù)概念性產(chǎn)流模型確定，缺少實測數(shù)據(jù)，需要經(jīng)過經(jīng)驗估計或參數(shù)率定來確定取值。本文對這部分參數(shù)首先進(jìn)行敏感性分析，定量評估其敏感程度，對于較為敏感的參數(shù)采取優(yōu)化率定的方法確定其取值，不敏感的參數(shù)則可以根據(jù)經(jīng)驗估算取值，敏感性分析所涉及的部分土壤參數(shù)含義及范圍見表2。

表2 部分土壤參數(shù)及取值范圍Tab.2 Partial soil parameters and value range

采用LH-OAT敏感性分析方法對這8個參數(shù)進(jìn)行敏感性分析，選用相對誤差RE、納什效率系數(shù)NSE、高水量誤差系數(shù)（fg）和低水量誤差系數(shù)（fd）共4 個誤差計算方法評估各個參數(shù)對日流量過程模擬的敏感程度。計算公式如下：

式中：T為計算的時間長度；為實測的流量序列；為模擬的流量序列為實測流量的平均值。

參數(shù)敏感性分析計算結(jié)果如表3 所示，根據(jù)各參數(shù)總敏感指標(biāo)對參數(shù)敏感性進(jìn)行分級，1 級表示最敏感，2～7 級表示較為敏感，8 級表示最不敏感［19，20］。表中最后一列為各參數(shù)敏感性等級最小值，以此作為該參數(shù)的全局敏感性等級。

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

由表3 可以看出，參數(shù)d1在4 種指標(biāo)下都是最敏感參數(shù)，上層土壤主要與表層產(chǎn)流有關(guān)，上層土壤厚度d1的大小會影響上層土壤含水量，厚度增加會導(dǎo)致蒸發(fā)增大，最終影響產(chǎn)流結(jié)果。參數(shù)B和d2的全局敏感性等級為2 級，參數(shù)B通過影響可變下滲曲線形狀來改變流域下滲率的空間分布，參數(shù)d2表示下層土壤厚度，但敏感度小于d1，說明閩江建陽流域下層土壤對產(chǎn)流過程的影響小于上層土壤。此外，DS、Dm和WS與計算基流過程有關(guān)，而d0影響頂薄層濕度和裸土蒸發(fā)，對徑流模擬結(jié)果也有一定的影響。

參數(shù)C在這4 種指標(biāo)下都不敏感，說明該參數(shù)對模型模擬結(jié)果的影響不大，可認(rèn)為是不敏感參數(shù)，在模型運行前根據(jù)經(jīng)驗方法取固定值即可。

研究結(jié)果與孫苗苗、張續(xù)軍和宋星原［21-23］對VIC 模型參數(shù)敏感性的分析結(jié)果相近，VIC 模型中共有7 個參數(shù)比較敏感，分別是B、DS、Dm、WS、d0、d1、d2，在參數(shù)率定中主要對這7個參數(shù)進(jìn)行率定，參數(shù)C為不敏感參數(shù)，取固定值2即可。

3.3 VIC模型參數(shù)率定與模型檢驗

選取閩江建陽流域1988-2000 年共13 年的數(shù)據(jù)進(jìn)行模擬計算，其中1988年為預(yù)熱期，率定期選擇的是1989-1996年共8年的數(shù)據(jù)，驗證期選擇的是1997-2000 年共4 年的數(shù)據(jù)。閩江建陽水文站1989-2000年的日流量過程為實測流量數(shù)據(jù)。

選用SCE-UA 算法對VIC 模型敏感性參數(shù)進(jìn)行參數(shù)率定，考慮到模型將d0設(shè)置的缺省值為0.1 m，且在后續(xù)計算中發(fā)現(xiàn)如果頂薄層厚度超過上層土壤厚度則會終止計算，因此本文將d0設(shè)置為默認(rèn)值，不參與模型的參數(shù)率定過程。最終率定的參數(shù)結(jié)果如表4所示。

表4 閩江建陽流域參數(shù)率定結(jié)果Tab.4 Parameter calibration results of Jianyang watershed of Minjiang River

4 誤差修正結(jié)果與分析

4.1 評價指標(biāo)

選用納什效率系數(shù)NSE、徑流相對誤差RE作為模型精度的評價指標(biāo)：

（1） 納什效率系數(shù)NSE，表示模擬和實測徑流過程的擬合度，NSE值越大，表示擬合度越高，如式（13）。

（2） 徑流相對誤差RE，表示徑流模擬總量和實測總量的偏差，從另一角度評估模型的模擬效果，RE值越小，偏差越小，如式（12）。

4.2 誤差修正結(jié)果

降雨過程是一個隨機(jī)的過程［24］，受很多因素的影響，比如：氣候條件、流域地形、植被和人類活動等，這些因素在現(xiàn)實中很難精準(zhǔn)把握和控制。水文遙測系統(tǒng)存在著信號接收和機(jī)械等方面的多重誤差影響，造成降雨觀測數(shù)據(jù)產(chǎn)生誤差［25］。因此，雨量數(shù)據(jù)的誤差，是水文模型誤差的重要來源。VIC 模型中需要研究流域日降雨數(shù)據(jù)制作氣象數(shù)據(jù)文件，降雨量的誤差會影響到產(chǎn)流量，最終影響模擬的流量值。

采用閩江建陽流域1989-2000 年共4 383 d 的數(shù)據(jù)進(jìn)行日徑流量的誤差修正，將率定好的參數(shù)代入模型中，對VIC模型中66個網(wǎng)格在日降雨量上增加一個單位，其他時段降雨量不變得到新的降雨量系列，按照系統(tǒng)微分響應(yīng)誤差修正的步驟，得到修正后的每日降雨量，將其重新代入VIC模型中進(jìn)行計算，最終得到修正后的閩江建陽流域的流量過程。并將此方法與常用的AR 模型誤差修正方法進(jìn)行比較。閩江建陽流域年均誤差修正結(jié)果如表5 所示。本文選用1989-2000 年共12 年的數(shù)據(jù)，從表5 中可以看出，系統(tǒng)微分響應(yīng)修正結(jié)果的納什效率系數(shù)均在0.85 以上，而AR 模型的納什效率系數(shù)均在0.85 以下，說明系統(tǒng)微分響應(yīng)誤差修正結(jié)果的納什效率系數(shù)高于AR 模型的值，此兩種修正方法的納什效率系數(shù)均高于未修正的模擬結(jié)果。系統(tǒng)微分響應(yīng)修正結(jié)果的徑流相對誤差也均低于AR 模型的結(jié)果，由此可以說明，按照年均誤差結(jié)果來看，用系統(tǒng)微分響應(yīng)修正方法的結(jié)果優(yōu)于傳統(tǒng)的AR 模型修正結(jié)果，能較好地模擬閩江建陽流域的徑流過程，具有較好的適用性。

表5 閩江建陽流域年均誤差修正結(jié)果Tab.5 Average annual error correction results of Jianyang basin of Minjiang River

比較分析修正前后的評價指標(biāo)，結(jié)果見表6。從表6 中可以看出，系統(tǒng)微分響應(yīng)修正后的納什效率系數(shù)（NSE）從修正前的0.752 增加到修正后的0.895，徑流相對誤差（RE）從修正前的7.89%降低到修正后的5.71%，AR 模型計算出的納什效率系數(shù)為0.807，徑流相對誤差為6.77%，系統(tǒng)微分響應(yīng)誤差修正結(jié)果的納什效率系數(shù)和徑流相對誤差均優(yōu)于AR 模型的結(jié)果。由此可以看出，采用系統(tǒng)微分響應(yīng)誤差修正方法對VIC 模型的降雨量進(jìn)行修正，修正后的閩江建陽流域的模擬精度得到顯著提升，說明系統(tǒng)微分響應(yīng)誤差修正方法對閩江建陽流域的日徑流過程模擬具有較好的適用性。

表6 閩江建陽流域誤差修正效果Tab.6 Error correction effect in Jianyang basin of Minjiang River

兩種誤差修正方法的實測流量與模擬流量過程如圖4和圖5，圖4 中為率定期的模擬流量結(jié)果過程線圖，圖5 為驗證期的模擬流量結(jié)果過程線圖。對比表5 的年均誤差修正結(jié)果和圖4和圖5中修正前后的流量過程線可以看出，1995和1998年的年實測徑流量大，日徑流峰值段量大且持續(xù)時間長，誤差修正后的峰值比修正前總體偏大，所以誤差修正后的徑流相對誤差為正值，分析其原因，可能與模型本身參數(shù)和結(jié)構(gòu)的局限性等因素有關(guān)。圖4 和圖5 中對比兩種誤差修正結(jié)果可以看出，系統(tǒng)微分響應(yīng)誤差修正的總體擬合度比AR 自回歸模型誤差修正方法好。

圖4 率定期（1989-1996年）系統(tǒng)微分響應(yīng)和AR模型的誤差修正結(jié)果過程線圖Fig.4 Error correction result hydrograph of system differential response and AR model of periodic systems（1989-1996）

圖5 驗證期（1997-2000年）系統(tǒng)微分響應(yīng)和AR模型的誤差修正結(jié)果過程線圖Fig.5 Error correction result hydrograph of system differential response and AR model of Validation period（1997-2000）

5 結(jié) 論

系統(tǒng)微分響應(yīng)誤差修正方法是一種結(jié)構(gòu)簡單、具有物理成因機(jī)理的預(yù)報誤差修正方法，該方法的理論基礎(chǔ)扎實，且不損失預(yù)見期，本文基于此，提出從輸入數(shù)據(jù)中的降雨量來進(jìn)行系統(tǒng)響應(yīng)誤差修正，將該方法應(yīng)用于VIC模型中，驗證系統(tǒng)微分響應(yīng)誤差修正方法在VIC模型中的適用性。得出如下結(jié)論：

（1）在建立閩江建陽流域的VIC模型后，應(yīng)用SCE-UA 算法先對VIC 模型中的敏感性參數(shù)進(jìn)行參數(shù)率定，再根據(jù)系統(tǒng)微分響應(yīng)誤差修正的原理對日徑流模擬結(jié)果進(jìn)行修正，修正后的NSE從0.752 提高至0.895，徑流相對誤差從7.89% 降低至5.71%。并與傳統(tǒng)的AR 自回歸模型誤差修正方法進(jìn)行對比，AR 模型修正結(jié)果的納什效率系數(shù)為0.807，徑流相對誤差為6.77%，計算結(jié)果表明，該方法對于VIC 模型具有較好的應(yīng)用效果，結(jié)果優(yōu)于傳統(tǒng)的AR 自回歸模型修正結(jié)果，閩江建陽流域的日徑流模擬精度得到顯著提升。

（2）在系統(tǒng)響應(yīng)誤差修正方法中，不僅有對降雨量的修正，還有對產(chǎn)流量和土壤初始含水量的修正，在后續(xù)的研究中可繼續(xù)研究對這些輸入量或者中間變量的修正效果，同時，系統(tǒng)響應(yīng)誤差修正方法也可應(yīng)用于洪水預(yù)報中，可深入研究VIC 模型中系統(tǒng)微分響應(yīng)方法對多場洪水的修正，探討系統(tǒng)微分響應(yīng)誤差修正方法與VIC模型的聯(lián)系。