系統微分響應誤差修正方法在VIC 模型中的應用與研究

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

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

0 引 言

水文模型是探索和認識流域水循環和水文過程的重要手段，在進行水文規律研究、水文預報、水資源規劃與管理、水文分析與計算和其他的生產實際問題中起著重要的作用［1，2］。水文模型按照對流域水文過程描述的離散程度分為集總式模型、分布式模型和半分布式模型，其中分布式模型是根據流域各處的氣候信息和下墊面特性的差異性，劃分流域為多個小單元，充分考慮了流域內部的各項氣候地理條件的差異性。在考慮流域各處氣候信息和下墊面特性的不同的基礎上，將流域劃分成多個基本單元，具有從機理上考慮降雨和下墊面條件空間分布不均勻對流域降雨徑流形成影響的功能［1］。VIC（Variable Infiltration Capacity）模型作為分布式模型代表，已被廣泛應用于不同尺度的流域水文模擬中，并都取得了不錯的模擬效果，說明VIC模型在流域水文模擬中適用性較好［3-5］。

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

1 研究方法

1.1 VIC模型

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

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

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

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

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

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

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

1.2 系統微分響應誤差修正方法

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

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

根據上述系統計算過程，可用以下公式來表示：

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

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

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

對式（4）兩邊同時求微分得系統輸出響應過程，如式（5）所示：

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

式（5）用矩陣表示為：

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

S矩陣的表達式為：

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

式（8）中，t=1，…，L，當t＜i時：

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

降雨量的系統微分響應修正步驟如下：

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

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

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

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

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

修正后的降雨量系列為：

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

2 研究流域概況及數據來源

2.1 研究流域概況

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

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

2.2 資料與數據

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

表1 原始數據說明Tab.1 Raw data description

3 建陽流域VIC模型構建與參數率定

3.1 建陽流域VIC模型構建

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

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

3.2 VIC模型參數敏感性分析

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

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

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

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

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

表3 VIC模型參數敏感性分析結果Tab.3 Parameters sensitivity results of the VIC model

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

參數C在這4 種指標下都不敏感，說明該參數對模型模擬結果的影響不大，可認為是不敏感參數，在模型運行前根據經驗方法取固定值即可。

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

3.3 VIC模型參數率定與模型檢驗

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

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

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

4 誤差修正結果與分析

4.1 評價指標

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

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

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

4.2 誤差修正結果

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

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

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

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

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

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

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

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

5 結 論

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

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

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