陳 長(zhǎng) 征，甘 容，楊 峰，左 其 亭

(1.鄭州大學(xué) 水利科學(xué)與工程學(xué)院，河南 鄭州 450001； 2.河南省地下水污染防治與修復(fù)重點(diǎn)實(shí)驗(yàn)室，河南 鄭州 450001； 3.河南省出山店水庫(kù)建設(shè)管理局，河南 信陽(yáng) 464043)

0 引 言

水文模型是對(duì)水流運(yùn)動(dòng)過(guò)程及其邏輯機(jī)理的概化，隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，不同種類(lèi)的水文模型陸續(xù)出現(xiàn)，逐步被應(yīng)用到水旱災(zāi)害防治、水文預(yù)報(bào)、水資源管理、水環(huán)境保護(hù)和水文變化成因分析等領(lǐng)域[1]。

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

本文以淮河出山店水庫(kù)為例，基于SWAT-CUP軟件下廣泛應(yīng)用的SUFI-2算法，設(shè)置NS、R2和PBIAS3種目標(biāo)函數(shù)，比較不同目標(biāo)函數(shù)下的參數(shù)敏感性；結(jié)合模型性能評(píng)價(jià)準(zhǔn)則，確立模型最佳模擬參數(shù)，建立了適用于該水庫(kù)流域的分布式SWAT模型，可為流域未來(lái)水資源形勢(shì)的科學(xué)把握和水庫(kù)的合理調(diào)度提供基礎(chǔ)支撐。

1 材料與方法

1.1 研究區(qū)概況

出山店水庫(kù)是淮河防洪系統(tǒng)的骨干工程，總庫(kù)容12.51億m3，控制流域面積約2 900 km2，于2019年投入使用，是淮河干流上首個(gè)控制性工程。水庫(kù)以上的淮河干流長(zhǎng)110 km，近幾十年流域內(nèi)下墊面變化較小，總體上人類(lèi)活動(dòng)對(duì)河道徑流影響較低，對(duì)水文數(shù)值模擬研究有利。流域內(nèi)水系發(fā)育，支流眾多，南側(cè)支流坡陡水急，北側(cè)支流彎多水淺；氣候?yàn)閬啛岷团瘻剡^(guò)渡帶季風(fēng)氣候，多年平均氣溫15.1 ℃。多年平均降水量1 130.9 mm，汛期集中在6～9月份，時(shí)空分布不均，水旱災(zāi)害頻生。研究區(qū)水文氣象站位置如圖1所示，其中大坡嶺水文站控制面積1 640 km2，其來(lái)水信息對(duì)出山店水庫(kù)的水情預(yù)報(bào)及洪水調(diào)度起著非常重要的作用。

圖1 研究區(qū)概況Fig.1 Overview of the study basin

1.2 數(shù)據(jù)介紹

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

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

圖3 土壤類(lèi)型Fig.3 Type of soil

表1 空間數(shù)據(jù)介紹

氣象數(shù)據(jù)源于中國(guó)氣象數(shù)據(jù)網(wǎng)和地方水文年鑒，研究區(qū)境內(nèi)氣象站有桐柏國(guó)家站和8個(gè)地方站，時(shí)間序列為1958～2018年，包括逐日降水、最高氣溫、最低氣溫、風(fēng)速、濕度和太陽(yáng)輻射等；水文數(shù)據(jù)為水庫(kù)上游大坡嶺水文站1960～2018年逐月徑流資料。

1.3 研究方法

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

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

1.3.1SWAT模型和率定參數(shù)

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

(1)

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

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

表2 模型參數(shù)

1.3.2目標(biāo)函數(shù)

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

(2)

(3)

(4)

式中：Qm為實(shí)測(cè)流量；Qs為模擬流量；i和n為樣本編號(hào)和總樣本數(shù)量。

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

敏感性評(píng)估有助于識(shí)別關(guān)鍵參數(shù)，降低參數(shù)不確定性影響，進(jìn)而提升模型優(yōu)化效率，采用SUFI-2提供的一次性O(shè)AT (One-At-a-Time) 敏感性分析和全局敏感性分析方法評(píng)估參數(shù)敏感性。其中，OAT屬于局部敏感性分析方法，做法是保持其他參數(shù)不變，調(diào)動(dòng)單個(gè)參數(shù)，觀(guān)察目標(biāo)函數(shù)變化。OAT操作簡(jiǎn)單，可快速估計(jì)參數(shù)敏感程度；全局敏感性分析是目標(biāo)函數(shù)在參數(shù)同時(shí)變化的情況下，對(duì)某個(gè)參數(shù)的平均變化率，方法如下[6]：

(1) 計(jì)算目標(biāo)函數(shù)g(b)的雅可比行列式J：

(5)

(2) 通過(guò)Gauss-Newton法和克萊默定理估計(jì)參數(shù)下限的協(xié)方差矩陣C：

(6)

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

(7)

1.3.4最佳模擬

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

(8)

(9)

式中：n*為95PPU包含的觀(guān)測(cè)樣本個(gè)數(shù)，n為觀(guān)測(cè)樣本總數(shù)；Qu和Ql分別為模擬值累積分布的97.5%和2.5%分位點(diǎn)；σ為觀(guān)測(cè)流量序列的標(biāo)準(zhǔn)差。

利用優(yōu)選的參數(shù)最終范圍模擬徑流，計(jì)算全部模擬結(jié)果的目標(biāo)函數(shù)，結(jié)合模型性能評(píng)價(jià)的分級(jí)準(zhǔn)則(見(jiàn)表3)[14]，篩選出目標(biāo)函數(shù)性能評(píng)價(jià)處于當(dāng)前最高級(jí)的所有模擬結(jié)果；從這個(gè)范圍中選取最佳模擬，便可避免使用單個(gè)目標(biāo)函數(shù)確立最佳模擬的弊端。

表3 模型性能的評(píng)價(jià)準(zhǔn)則

2 結(jié)果和討論

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

2.1.1OAT結(jié)果分析

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

圖5 調(diào)動(dòng)單個(gè)參數(shù)時(shí)目標(biāo)函數(shù)的變化特征Fig.5 Variation characteristics of the objective functions when changing single parameter

2.1.2全局敏感性分析

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

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

綜上可發(fā)現(xiàn)：參數(shù)CN2，CANMX，SOL_BD和GW_DELAY的敏感性最強(qiáng)，是參數(shù)率定時(shí)需要被重點(diǎn)關(guān)注的對(duì)象；參數(shù)ALPHA_BF和GWQMN被R2判定為敏感，被NS判定為不敏感，表明ALPHA_BF和GWQMN對(duì)平水期和枯水期的水文過(guò)程影響更強(qiáng)；參數(shù)GW_REVAP，ESCO，RCHRG_DP和SOL_AWC只被PBIAS判定為敏感，原因可能是參數(shù)引起的豐枯水季流量偏差的作用相當(dāng)，在徑流總體走勢(shì)上體現(xiàn)不明顯?；赑BIAS計(jì)算的P-Value基本處于3個(gè)目標(biāo)函數(shù)中的最低水平，基于R2和NS計(jì)算的P-Value互有高低，表明PBIAS對(duì)參數(shù)變化的響應(yīng)度最強(qiáng)，與OAT結(jié)果一致。

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

2.2 確定最佳模擬

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

表4 參數(shù)在NS、R2、PBIAS為目標(biāo)函數(shù)時(shí)的全局敏感性統(tǒng)計(jì)

2.3 模型應(yīng)用

將最終參數(shù)范圍和最佳模擬應(yīng)用到流域1980～2018年的逐月徑流模擬，分1980～1999年和2000～2018年兩個(gè)驗(yàn)證期，評(píng)估最佳參數(shù)范圍和最佳模擬的性能。

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

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

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

表5是最佳模擬在率定期、驗(yàn)證期01和驗(yàn)證期02的性能評(píng)價(jià)，由表可知不同目標(biāo)函數(shù)的評(píng)價(jià)標(biāo)準(zhǔn)下，率定期和驗(yàn)證期的性能評(píng)價(jià)都為優(yōu)秀，表明最佳模擬在長(zhǎng)時(shí)間序列中可以取得穩(wěn)定良好的性能，基于NS、R2和PBIAS3種目標(biāo)函數(shù)的評(píng)價(jià)分級(jí)準(zhǔn)則遴選最佳模擬的方法可行，這種方法可避免評(píng)價(jià)指標(biāo)性能差別大時(shí)參數(shù)的反復(fù)迭代，在一定程度上提高了模型優(yōu)化的效率。

表5 不同時(shí)期最佳模擬的性能評(píng)價(jià)

3 結(jié) 論

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

(1) 基于不同的目標(biāo)函數(shù)，確立的參數(shù)敏感性有所差異。NS，R2和PBIAS3種目標(biāo)函數(shù)中，PBIAS對(duì)參數(shù)變化響應(yīng)度最高，基于PBIAS確定的敏感參數(shù)最多，包括CN2，CANMX，SOL_BD，GW_DELAY，GWQMN， GW_REVAP，ESCO，RCHRG_DP和SOL_AWC；同類(lèi)研究中，建議使用PBIAS作為目標(biāo)函數(shù)篩選敏感參數(shù)，但當(dāng)率定參數(shù)的數(shù)目較大時(shí)，建議使用NS作為目標(biāo)函數(shù)，因?yàn)镹S能夠更快速地實(shí)現(xiàn)參數(shù)降維，犧牲部分模型性能提高優(yōu)化效率。

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

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

(4) 出山店水庫(kù)正處于試運(yùn)行階段，本文的研究結(jié)果對(duì)于水庫(kù)流域開(kāi)展水資源評(píng)價(jià)具有重要意義，特別是對(duì)流域未來(lái)氣候下水資源變化形勢(shì)的分析工作有一定參考價(jià)值。