基于SWAT的徑流模擬參數優化方案及不確定性分析

陳 長 征，甘 容，楊 峰，左 其 亭

(1.鄭州大學 水利科學與工程學院，河南 鄭州 450001； 2.河南省地下水污染防治與修復重點實驗室，河南 鄭州 450001； 3.河南省出山店水庫建設管理局，河南 信陽 464043)

0 引 言

水文模型是對水流運動過程及其邏輯機理的概化，隨著計算機技術的高速發展，不同種類的水文模型陸續出現，逐步被應用到水旱災害防治、水文預報、水資源管理、水環境保護和水文變化成因分析等領域[1]。

參數率定是水文模型優化的重要內容，但模型中通常涉及大量參數，難以同時保障全部參數的精確性，因此，評估參數敏感性、遴選關鍵參數對模型優化具有重要意義[2-3]。順序不確定擬合算法(Sequential Uncertainty Fitting Version2，SUFI-2)是SWAT(Soil and Water Assessment Tool)模型率定的常用工具，劉偉等[4]比較了SUFI-2和廣義似然不確定性估計算法GLUE(Generalized Likelihood Uncertainty Estimation)在參數率定和不確定性評估等方面的性能，結果表明SUFI-2算法的率定效率更高。周帥等[5]利用SUFI-2識別敏感參數及其置信區間，進而量化了參數交互作用不確定性對徑流模擬的敏感性。Abbaspour等利用SUFI-2篩選出敏感參數后，通過超拉丁立方體抽樣構造參數組合，將目標函數如納什效率NS(Nash-Sutcliffe)的最優值對應參數組合的模擬結果選為最佳模擬，計算其決定系數R2(Coefficient of Determination)、百分比偏差PBIAS(Percentage Bias)、克林-古普塔效率KGE(Kling-Gupta Efficiency)等指標，依據各指標評估模擬性能[6]。應用SUFI-2的研究中，多數將NS作為目標函數確立最佳模擬[7-8]，但是除NS外的評價指標取值不確定性大，會出現NS評價的模型性能優秀，其他指標評價的性能差的現象，需要反復調參才能避免，降低了模型優化速率。

本文以淮河出山店水庫為例，基于SWAT-CUP軟件下廣泛應用的SUFI-2算法，設置NS、R2和PBIAS3種目標函數，比較不同目標函數下的參數敏感性；結合模型性能評價準則，確立模型最佳模擬參數，建立了適用于該水庫流域的分布式SWAT模型，可為流域未來水資源形勢的科學把握和水庫的合理調度提供基礎支撐。

1 材料與方法

1.1 研究區概況

出山店水庫是淮河防洪系統的骨干工程，總庫容12.51億m3，控制流域面積約2 900 km2，于2019年投入使用，是淮河干流上首個控制性工程。水庫以上的淮河干流長110 km，近幾十年流域內下墊面變化較小，總體上人類活動對河道徑流影響較低，對水文數值模擬研究有利。流域內水系發育，支流眾多，南側支流坡陡水急，北側支流彎多水淺；氣候為亞熱和暖溫過渡帶季風氣候，多年平均氣溫15.1 ℃。多年平均降水量1 130.9 mm，汛期集中在6～9月份，時空分布不均，水旱災害頻生。研究區水文氣象站位置如圖1所示，其中大坡嶺水文站控制面積1 640 km2，其來水信息對出山店水庫的水情預報及洪水調度起著非常重要的作用。

圖1 研究區概況Fig.1 Overview of the study basin

1.2 數據介紹

模型構建需要的基礎空間數據包括：數字高程模型DEM(Digital Elevation Model)、土地利用(見圖2)和土壤類型(見圖3)。土地利用采用的是資源環境科學與數據中心發布的《2010年中國土地利用現狀遙感監測數據》。由DEM可知：研究區西北平均高程700～1 200 m，地勢較高，山巒眾多，中東平均高程100～300 m，地勢平起伏小。由土地利用數據可知：耕地和林地占研究區總面積比例超過90%，其中耕地多種植水稻和小麥，水稻田約占70%耕地。由土壤類型數據可知：研究區中東分布累積性紅土和不飽和始成土，分別占26%和23%，西北分布石灰性粗骨土和飽和黏性土，約占14%和20%。空間數據來源及分辨率如表1所列。

圖2 土地利用類型(2010年)Fig.2 Type of land use (2010)

圖3 土壤類型Fig.3 Type of soil

表1 空間數據介紹

氣象數據源于中國氣象數據網和地方水文年鑒，研究區境內氣象站有桐柏國家站和8個地方站，時間序列為1958～2018年，包括逐日降水、最高氣溫、最低氣溫、風速、濕度和太陽輻射等；水文數據為水庫上游大坡嶺水文站1960～2018年逐月徑流資料。

1.3 研究方法

研究包括3個主要部分：① 參數敏感性分析，② 參數范圍優選和最佳模擬的確定，③ 模型驗證，研究流程如圖4所示。

圖4 研究流程Fig.4 Flow chart of the study

1.3.1SWAT模型和率定參數

SWAT是基于物理的分布式水文模型，它將流域劃為多個子流域SUBs(Subbasins)，再依據土地利用、土壤和坡度特征將SUB劃為水文響應單元HRUs(Hydrological Response Units)。模型中流域水文循環分為兩個部分：第一部分為陸地階段，控制進入河道的水、泥沙和營養物等；第二部分為河道演算階段，計算水、泥沙等沿著河道運動至流域出口的過程。基于水量平衡方程計算陸地產流，如公式(1)，河網中的水經過主河道演算和水庫演算，最終匯流至流域出口[9]。

(1)

式中：SWt為土壤最終含水量，mm；SW0為土壤前期含水量，mm；t為時間步長，d；Rd為第i天降水量，mm；Qsur為第i天地表徑流量，mm；Ea為第i天蒸發量，mm；Ws為第i天土壤剖面底層的滲流和側流量，mm；Qg為第i天地下水出流量，mm。

選用14種水文相關參數參與模型率定，如表2所列，其中數值空間分布具有異質性的參數(如CN2)，將其調參方式設為基于原始值進行縮放，目的是保留其空間分布特性，使流域建模的特征和實際更加貼近[10-12]。

表2 模型參數

1.3.2目標函數

目標函數是SUFI-2的關鍵輸出，建立目標函數與模擬變量的聯系可分析參數敏感性。選取納什效率NS，確定系數R2和百分比偏差PBIAS作為本次模擬的目標函數，分別由公式(2)～(4) 表示。其中，NS可體現模擬徑流和實測徑流的總體貼合情況，但會賦予洪峰段更高的計算權重，容易忽視平水期或枯水期的擬合情況；R2是模擬值與實測值的擬合優度，弊端是R2對模擬值整體偏高或偏低的偏差響應不明顯；百分比偏差PBIAS可體現模擬值和實測值的累積偏差，當模擬水文過程與實際趨勢貼合情況良好時，PBIAS可更加精確地評估模型總水量平衡的效果[13]。

(2)

(3)

(4)

式中：Qm為實測流量；Qs為模擬流量；i和n為樣本編號和總樣本數量。

1.3.3參數敏感性分析

敏感性評估有助于識別關鍵參數，降低參數不確定性影響，進而提升模型優化效率，采用SUFI-2提供的一次性OAT (One-At-a-Time) 敏感性分析和全局敏感性分析方法評估參數敏感性。其中，OAT屬于局部敏感性分析方法，做法是保持其他參數不變，調動單個參數，觀察目標函數變化。OAT操作簡單，可快速估計參數敏感程度；全局敏感性分析是目標函數在參數同時變化的情況下，對某個參數的平均變化率，方法如下[6]：

(1) 計算目標函數g(b)的雅可比行列式J：

(5)

(2) 通過Gauss-Newton法和克萊默定理估計參數下限的協方差矩陣C：

(6)

(3) 通過雅可比矩陣的中列平均值計算參數全局敏感度S：

(7)

1.3.4最佳模擬

基于95%置信水平上的不確定性區95PPU(95% prediction uncertainty)的評價準則，優選合適的參數范圍。95PPU有兩個特征指標：P-factor和R-factor。其中，P-factor是觀測數據被95PPU包絡的百分比，R-factor是95PPU的平均厚度，如式(8)和(9)所示。P-factor越大且R-factor越小，表明更窄的95PPU可包含更多的觀測值，即95PPU不確定性越小。理論上P-factor=1且R-factor=0時，模擬值與實測值完全吻合，此時模型參數完全和真實流域一致。實際建模中認為P-factor>0.7且R-factor<1時，模擬的不確定性可接受[6]。

(8)

(9)

式中：n*為95PPU包含的觀測樣本個數，n為觀測樣本總數；Qu和Ql分別為模擬值累積分布的97.5%和2.5%分位點；σ為觀測流量序列的標準差。

利用優選的參數最終范圍模擬徑流，計算全部模擬結果的目標函數，結合模型性能評價的分級準則(見表3)[14]，篩選出目標函數性能評價處于當前最高級的所有模擬結果；從這個范圍中選取最佳模擬，便可避免使用單個目標函數確立最佳模擬的弊端。

表3 模型性能的評價準則

2 結果和討論

2.1 參數敏感性分析

2.1.1OAT結果分析

模擬流域1960～1979年逐月徑流過程，圖5為OAT方法下NS，R2和PBIAS的變化特征，由圖5可知：參數CN2、SOL_BD，SOL_AWC，CANMX，GWQMN，RCHRG_DP，ESCO和GW_REVAP改變時，目標函數的變化相對明顯，參數的敏感性較強。其中，反映總體水量平衡的PBIAS對參數的變化反應最為靈敏，NS和R2對參數變化反應相對遲鈍。

圖5 調動單個參數時目標函數的變化特征Fig.5 Variation characteristics of the objective functions when changing single parameter

2.1.2全局敏感性分析

表4顯示了基于不同目標函數計算的參數全局敏感度，其中T檢驗的統計值t-Stat(Statistical value of the T-test)是標準化后的取值，其絕對值越大代表參數越敏感；圖6 為t-Stat的顯著性檢驗結果，其中P-Value(Probability of rejecting the original assumption)是T檢驗拒絕原假設(參數不敏感)的概率，當P-Value≤0.05時，拒絕原假設即參數敏感。由圖6知，參數CN2，CANMX，SOL_BD和GW_DELAY基于NS，R2，PBIAS計算的P-Value全小于0.05，即NS、R2、PBIAS都將其判定為敏感參數；基于R2判定為敏感的參數有CN2，CANMX，SOL_BD，GW_DELAY，ALPHA_BF和GWQMN；基于PBIAS判定為敏感的參數有CN2，CANMX，SOL_BD，GW_DELAY，GWQMN，GW_REVAP，ESCO，RCHRG_DP和SOL_AWC。

圖6 全局敏感度顯著性檢驗P-ValueFig.6 P-value of significance test for the global sensitivity

綜上可發現：參數CN2，CANMX，SOL_BD和GW_DELAY的敏感性最強，是參數率定時需要被重點關注的對象；參數ALPHA_BF和GWQMN被R2判定為敏感，被NS判定為不敏感，表明ALPHA_BF和GWQMN對平水期和枯水期的水文過程影響更強；參數GW_REVAP，ESCO，RCHRG_DP和SOL_AWC只被PBIAS判定為敏感，原因可能是參數引起的豐枯水季流量偏差的作用相當，在徑流總體走勢上體現不明顯。基于PBIAS計算的P-Value基本處于3個目標函數中的最低水平，基于R2和NS計算的P-Value互有高低，表明PBIAS對參數變化的響應度最強，與OAT結果一致。

全局敏感性分析確定了11個敏感參數，比OAT多3個，原因是OAT的前提是參數互相獨立，而實際中參數敏感性會隨其他參數改變而發生變化。因此，全局敏感性分析得出的參數敏感性結果更加可靠，但面對參數極多情況，OAT仍是初步掌握參數敏感度的實用方法，必要時可采用OAT迅速實現參數降維。

2.2 確定最佳模擬

采用PBIAS判斷出的敏感參數參與模型率定，基于P-factor>0.7且R-factor<1的準則降低模擬不確定性，直到效果滿意，參數最終范圍如表4所列。圖7是參數最終范圍上500次隨機模擬結果的目標函數的空間散點分布和坐標面投影。其中，NS-R2相關系數為0.67，NS-PBIAS相關系數為0.52，R2-PBIAS相關系數為0.07。這表明利用3個目標函數評價模型性能時，NS的綜合能力最強，兼具了R2和PBIAS的部分特性。由圖7可發現NS即將取到最優值時，R2和PBIAS的取值范圍較大，NS最優值對應的PBIAS取值在優秀等級范圍外。PBIAS評價的模型性能包括優秀、良好、一般和較差4個等級，NS評價的模型性能包括優秀和良好2個等級，R2評價的模型性能只有優秀等級，表明PBIAS確實對參數變化響應最靈敏，變化幅度最大。綜上，采用限制評定等級的方法遴選最佳模擬，首先，限定NS，R2，PBIAS都處于優秀等級范圍內，然后采用性能評價綜合能力最高的NS作為尋優標準，以其最高值對應模擬為最佳模擬。

表4 參數在NS、R2、PBIAS為目標函數時的全局敏感性統計

2.3 模型應用

將最終參數范圍和最佳模擬應用到流域1980～2018年的逐月徑流模擬，分1980～1999年和2000～2018年兩個驗證期，評估最佳參數范圍和最佳模擬的性能。

由率定期和驗證期的降水分布可知(見圖8)，率定期的逐月降水高于驗證期01和02；由大坡嶺水文站模擬逐月徑流過程可知(見圖9)：率定期和驗證期01和02的P-factor分別為0.86，0.84，0.86，均大于0.7；R-factor分別為0.65，0.77，0.78，均小于1.0。表明模型輸出的不確定性較小，可接受。率定期、驗證期01和02的P-factor基本一致，表明率定出的參數范圍在率定期和驗證期01和02捕捉流域特征的能力相當；率定期的R-factor明顯小于兩個驗證期，可能是由降水分布差異導致的。降水愈多，降水對徑流的影響越強，徑流和降水的關系越密切，參數對SWAT的影響相對越低，SWAT模擬性能更加優秀[15-16]。率定期與驗證期的參數范圍相同，但率定期降水較多，其輸出的95PPU整體性能更優，表現為模擬徑流過程向實測徑流過程貼近，使95PPU的平均寬度更小。

圖8 逐月降水分布箱形圖Fig.8 The box chart of monthly precipitation

圖9 大坡嶺站月徑流模擬過程Fig.9 Monthly runoff simulation process of the Dapoling Hydrological Station

表5是最佳模擬在率定期、驗證期01和驗證期02的性能評價，由表可知不同目標函數的評價標準下，率定期和驗證期的性能評價都為優秀，表明最佳模擬在長時間序列中可以取得穩定良好的性能，基于NS、R2和PBIAS3種目標函數的評價分級準則遴選最佳模擬的方法可行，這種方法可避免評價指標性能差別大時參數的反復迭代，在一定程度上提高了模型優化的效率。

表5 不同時期最佳模擬的性能評價

3 結 論

本文基于設立NS，R2和PBIAS3個目標函數的SUFI-2算法確定了出山店水庫流域的敏感水文參數，率定出參數最終范圍和最佳模擬，結論如下：

(1) 基于不同的目標函數，確立的參數敏感性有所差異。NS，R2和PBIAS3種目標函數中，PBIAS對參數變化響應度最高，基于PBIAS確定的敏感參數最多，包括CN2，CANMX，SOL_BD，GW_DELAY，GWQMN， GW_REVAP，ESCO，RCHRG_DP和SOL_AWC；同類研究中，建議使用PBIAS作為目標函數篩選敏感參數，但當率定參數的數目較大時，建議使用NS作為目標函數，因為NS能夠更快速地實現參數降維，犧牲部分模型性能提高優化效率。

(2) SWAT的徑流模擬不確定性和降水分布相關，降水愈豐，降水對徑流模擬的影響權重愈大，模擬的擬合效果越好；因此，固定參數范圍后，率定期和驗證期的P-factor變化不大，但降水越多，95PPU越向著性能更高的方向集中，使平均寬度R-factor變窄，從而表現為豐水期的模擬不確定性更小。

(3) 利用NS，R2和PBIAS評價參數最終范圍內的全部模擬結果時，PBIAS評價的性能等級的不確定性最大，NS評價性能的綜合能力最高。因此，建議通過限制PBIAS縮小選擇范圍，然后由高到低挑選NS的方法確定最佳模擬，可避免反復調參，提高效率。

(4) 出山店水庫正處于試運行階段，本文的研究結果對于水庫流域開展水資源評價具有重要意義，特別是對流域未來氣候下水資源變化形勢的分析工作有一定參考價值。