GLUE框架下似然函數(shù)對水文模型不確定性的影響

姜倩妮, 李占玲, 張永勇

(1.中國地質(zhì)大學(xué)(北京) 水資源與環(huán)境學(xué)院， 北京 100083； 2.中國科學(xué)院地理科學(xué)與資源研究所 陸地水循環(huán)及地表過程重點實驗室， 北京 100101)

1 研究背景

水文模型是探索自然界水循環(huán)過程的重要手段，如今已經(jīng)廣泛應(yīng)用在水資源管理、水文預(yù)報等方面。隨著對水文循環(huán)過程的不斷了解和計算機性能的不斷提高，水文過程的模擬已經(jīng)越來越精確。但由于更為復(fù)雜水文過程的引入，模型對水文物理過程描述的逐漸增加、模型結(jié)構(gòu)呈現(xiàn)復(fù)雜化趨勢[1]，由此也導(dǎo)致模型參數(shù)的不斷增加。謹慎地率定參數(shù)和分析參數(shù)的不確定性，對于提高水文預(yù)報精度具有重要意義。

普適似然不確定性估計法(GLUE：Generalized Likelihood Uncertainty Estimation)目前已廣泛應(yīng)用在參數(shù)的率定和不確定性分析中[2-4]，該方法是Beven和Binley在1992年基于“異參同效”現(xiàn)象提出的[5]。與其它不確定性分析方法相比，GLUE方法操作簡單并可通過靈活定義似然函數(shù)來分離可行和不可行參數(shù)組[6]。似然函數(shù)是GLUE方法的重要組成部分，對參數(shù)的率定和不確定性分析產(chǎn)生一定的影響。目前單目標GLUE方法中，Nash-Sutcliffe系數(shù)(NSCE)是最常用的似然函數(shù)[7-8]，但NSCE僅能描述一部分的水文過程，如劉艷麗等[9]將NSCE單個似然判據(jù)擴展為洪峰誤差、洪量、峰現(xiàn)時間、NSCE4個目標的多準則似然判據(jù)時，能更好地反映模型的實際不確定性情況。考慮多個似然函數(shù)的GLUE方法雖然綜合描述了多個水文過程，但權(quán)重的主觀分配進一步增加了參數(shù)的不確定性[10-11]。

目前基于單目標的GLUE方法評價參數(shù)不確定性的研究已經(jīng)有很多，但不同似然函數(shù)對參數(shù)敏感性和不確定性的影響評價研究還比較少。本文以淮河上游大坡嶺控制流域為例，探索了GLUE框架下NSCE和水量平衡系數(shù)作為似然函數(shù)對流域水循環(huán)模型參數(shù)敏感性和不確定性的影響，旨在為GLUE方法中似然函數(shù)的選取以及為多目標GLUE方法權(quán)重的分配提供參考依據(jù)。

2 水循環(huán)模型與GLUE方法

2.1 流域水循環(huán)模型

模型采用ZhangYongyong等[12]于2012-2015年期間自主開發(fā)的流域水循環(huán)系統(tǒng)模型(Hydrological,EcologicalandwaterQualityModel:HEQM)。該模型以水和營養(yǎng)物(碳、氮和磷)循環(huán)作為聯(lián)系流域內(nèi)與水相關(guān)各過程的紐帶，綜合考慮水文過程和營養(yǎng)物循環(huán)在陸面、土壤、植被和河流水體中的相互作用關(guān)系及閘壩調(diào)控、排污等人類活動對水和營養(yǎng)物循環(huán)的影響[12]。HEQM包含八大模塊，本文主要采用水文循環(huán)模塊。其關(guān)鍵的水量平衡方程為：

Pi+SWi=SWi+1+Rsi+Eai+Rssi+Rbsi+Ini

(1)式中：P為降水，mm；SW為土壤濕度，mm；Ea為實際蒸散發(fā)，mm，包括土壤蒸發(fā)(Es, mm)和植被蒸散發(fā)(Ep,mm)；Rs、Rss和Rbs分別為地表徑流、壤中流和基流，mm；In為植被截留，mm；i為時間步長，d。

方程(1)中地表徑流Rs的計算采用時變增益產(chǎn)流方程(Time Variant Gain Model: TVGM)[13]；壤中流Rss和基流Rbs分別與上下層土壤濕度呈線性關(guān)系；潛在蒸散發(fā)的計算采用Hargreaves方法[14]；實際蒸散發(fā)包括植被蒸騰和土壤蒸發(fā)，植被蒸騰通過潛在蒸散發(fā)量和葉面積指數(shù)獲得，土壤蒸發(fā)通過潛在蒸發(fā)量和表層土壤殘余度獲得；上下層水流之間的運動采用土壤蓄水演算法；河道匯流演算采用馬斯京根法或動力波方程。

2.2 GLUE方法

作為分析參數(shù)不確定性的重要方法，GLUE方法的重要觀點是，模型模擬結(jié)果的好壞不是由模型中的某個參數(shù)所決定，而是由一組模型參數(shù)來決定。一般的參數(shù)率定只選出一組最優(yōu)的參數(shù)組作為模型的參數(shù)，但是因為“異參同效”現(xiàn)象，“最優(yōu)參數(shù)”很難獲得，GLUE方法通過選取似然函數(shù)，設(shè)置閾值，篩選“可行參數(shù)組”，分析模擬結(jié)果的不確定性[5]，可以避免“異參同效”現(xiàn)象帶來的影響。

GLUE方法參數(shù)敏感性和模擬結(jié)果的不確定性分析步驟如下：

(1)選擇似然函數(shù)。似然函數(shù)是評價水文過程模擬好壞的重要指標。本文選取NSCE和水量平衡系數(shù)作為似然函數(shù)，分別從時間和水量兩個方面描述水文過程，其計算公式如下：

Nash-Sutcliffe系數(shù)：

(2)

水量平衡系數(shù)：

(3)

(2)模型參數(shù)范圍的確定和先驗分布的選取。對于模型參數(shù)的選取和取值范圍的確定，參考相關(guān)文獻[12]并根據(jù)模型參數(shù)的物理特征，確定參數(shù)的初始范圍。由于先驗信息不足，無法準確描述先驗分布，通常選擇均勻分布作為先驗分布。

(3)參數(shù)抽樣及其優(yōu)化。參數(shù)優(yōu)化采用SCE-UA單目標自動優(yōu)化算法。SCE-UA優(yōu)化算法可以解決復(fù)雜水文模型參數(shù)的優(yōu)選問題，由DuanQingyun等[15]結(jié)合單純性法、隨機搜索和生物競爭進化方法的優(yōu)點而研發(fā)的。

(4)確定閾值。閾值的選取帶有一定的主觀性，本文將NSCE的閾值設(shè)定為0.7；水量平衡系數(shù)在-0.05～0.05之間。如果似然值滿足閾值條件，其對應(yīng)的參數(shù)集稱為可行參數(shù)集，反之則為不可行參數(shù)集。

(5)參數(shù)的不確定性分析。根據(jù)參數(shù)組模擬情況確定可行參數(shù)組，繪制模型參數(shù)與似然值的散點圖、可行參數(shù)組頻率直方圖，定性分析模型參數(shù)的不確定性。在散點圖中，若參數(shù)似然分布接近于原始均勻分布，則參數(shù)對目標函數(shù)的影響不顯著，為不敏感性參數(shù)；反之，則為敏感性參數(shù)。在頻率直方圖中，各參數(shù)分布的形狀表明估計的不確定性的頻率，參數(shù)在某一值處的頻率越低說明參數(shù)具有越多的不確定性。

(6)模擬結(jié)果不確定性評估。給每個可行參數(shù)組一個權(quán)重，模型輸出不確定性用權(quán)重化后的參數(shù)組決定的輸出分位數(shù)表示，即基于95%置信區(qū)間對應(yīng)的分位數(shù)組成的預(yù)測區(qū)間。同時還要對模擬結(jié)果的不確定性進行評估。

對于不確定性結(jié)果的區(qū)間估計，本文選取比較常用的2個指標，即覆蓋率CS和平均對稱度S[16]。CS表示預(yù)測區(qū)間包含實測流量的比例，通常用百分比表示，CS越接近100%則說明模擬值接近實測值的可能性越大，模擬效果越好。在絕對對稱的情形及理想情況下，S值為0；若S<0.5，則說明在平均水平上流量過程線在預(yù)測不確定性區(qū)間內(nèi)；若S≥0.5，則說明平均水平上流量過程線在預(yù)測不確定性區(qū)間之外，而且S越大則不確定性區(qū)間的性能越差。其公式如下：

(4)

式中：nQin為實測值落在置信區(qū)間的個數(shù)；n為可行參數(shù)的個數(shù)。

3 實例分析

3.1 研究區(qū)概況與水循環(huán)模型構(gòu)建

本文研究區(qū)為淮河上游的大坡嶺流域(E113°15′～114°46′，N31°31′～32°43′)，流域面積為1 573.3km2，見圖1。流域內(nèi)分布有回龍寺、黃崗、毛集、潘莊、固縣等13個雨量站，多年平均降雨量為918mm，其中50%集中在汛期( 6-9月)，而多年平均氣溫11～16℃。流域地形呈現(xiàn)山脈較多，坡度較大的特點。流域出口大坡嶺水文站歷年最高水位為104. 86m，最大洪峰流量為4 200m3/s[17]。

研究搜集了數(shù)字高程模型(DEM: 90m×90m)、 土地利用(1km×1km)等GIS數(shù)據(jù)，包括流域內(nèi)13個雨量站點1994-2002年日降水序列、以及大坡嶺水文站日徑流序列。利用ArcGIS平臺提取HEQM所需的流域?qū)傩詳?shù)據(jù)；并利用泰森多邊形方法求得該流域的面雨量序列，驅(qū)動已構(gòu)建的模型，進而模擬大坡嶺斷面日徑流過程。選取1994年為模型預(yù)熱期，1995-1999年為模型率定期，2000-2002年為模型驗證期。HEQM模型水文相關(guān)的參數(shù)主要有9個，其中與土壤有關(guān)的參數(shù)5個：土壤穩(wěn)定下滲率(fc)、上層飽和土壤含水量(WM)、土壤水出流系數(shù)(Kr) 、田間持水量(WMc)，土壤蒸散發(fā)補充系數(shù)(ketp)；與地表水有關(guān)的參數(shù)2個：產(chǎn)流參數(shù)g1和g2；地下水有關(guān)參數(shù)2個：地下水出流系數(shù)(Krg)，水滯留時間(Tg)[12]。具體的參數(shù)范圍見表1。

圖1 大坡嶺流域圖

類別參數(shù)初始取值范圍NSCE閾值為0.7參數(shù)取值范圍水量平衡系數(shù)滿足條件的為在±0.05之間土壤穩(wěn)定下滲率(fc)/(mm·h-1)0～120.0038.00～120.000～120.00飽和土壤含水量(WM)0.45～0.750.45～0.500.45～0.70土壤土壤水出流系數(shù)(Kr)0～1.000.36～0.720～1.00土壤田間持水量(WMc)0.20～0.450.26～0.360.20～0.45土壤蒸散發(fā)補充系數(shù)(ketp)0～2.000.50～0.810～1.90地表水g1調(diào)整系數(shù)0～3.001.60～3.000～3.00g2調(diào)整系數(shù)0～3.001.20～2.100～3.00地下水地下水出流系數(shù)(Krg)0～1.000.54～1.000～1.00地下水滯留時間(Tg)/d1.00～100.001.00～3.201.00～97.00

3.2 以效率系數(shù)為似然函數(shù)的參數(shù)不確定性分析

3.2.1 參數(shù)敏感性分析 利用SCE-UA算法，隨機均勻抽取的10 000組參數(shù)分別運行HEQM模型，模擬大坡嶺流域日徑流量的過程，最終獲得NSCE大于0.7的參數(shù)似然函數(shù)散點圖、可行參數(shù)頻率直方圖，見圖2，可行參數(shù)組的取值范圍見表1。篩選的9個參數(shù)形成的散點圖均有明顯的峰值，可作為模型的敏感性參數(shù)。其中可行參數(shù)的取值范圍與初始范圍變化最為明顯的是與土壤有關(guān)的參數(shù)fc、ketp和與地下水有關(guān)的參數(shù)Tg，分別由原來的0～120.00mm/h、0～2.00和1～100.00d，縮小到38.00～120.00mm/h、0.50～0.81和1～3.20d。因此，以NSCE為似然函數(shù)的GLUE方法可有效縮小部分參數(shù)的取值范圍；另外，這9個參數(shù)均具有良好的單峰分布，其中fc的取值主要集中在80～100mm/h，ketp的取值在0.50～0.75，其他7個參數(shù)的取值均在一定范圍內(nèi)比較穩(wěn)定。通過對各組可行參數(shù)的頻率分析得：以NSCE為似然函數(shù)選取的可行參數(shù)值出現(xiàn)較高頻次的范圍比較集中，這進一步減小了參數(shù)的不確定性。由于篇幅有限，本文只展示土壤穩(wěn)定下滲率和蒸散發(fā)補充系數(shù)相關(guān)結(jié)果，見圖2。

3.2.2 模擬結(jié)果的不確定性分析 根據(jù)篩選的可行參數(shù)組模擬每組參數(shù)對應(yīng)的徑流過程，計算獲得95%置信度下NSCE作為似然函數(shù)時大坡嶺日徑流模擬的不確定性區(qū)間，見圖3。結(jié)合不確定性評價指標，可以看出：在率定期和驗證期模擬的日徑流序列區(qū)間涵蓋了26.2%和15.8%的實測值；當觀測值很大或者很小時，均不能將其涵蓋，特別是在驗證期。這可能是由于模型結(jié)構(gòu)的不確定性或者人為選取閾值的原因。

圖2 參數(shù)-NSCE散點圖與可行參數(shù)組頻率直方圖

圖3 95%置信度下NSCE系數(shù)作為似然函數(shù)的率定期和驗證期的不確定性范圍

3.3 以水量平衡系數(shù)為似然函數(shù)的參數(shù)不確定性分析

3.3.1 參數(shù)的敏感性分析 參數(shù)-水量平衡系數(shù)散點圖與可行參數(shù)組頻直方圖見圖4。從水量平衡系數(shù)在±0.05之間的可行參數(shù)組的取值范圍可以看出，絕大部分滿足似然判據(jù)的參數(shù)取值范圍與原始取值范圍相同，只有ketp和Tg從原始的0～2.00、1～100.00d縮小至0～1.90和1～97.00d。從圖4(a)可以看出，僅fc為敏感性參數(shù)；而從圖4(b)可以看出，fc取值范圍在80.00～100.00 mm/h之間時，ketp的取值主要集中在0.50～0.75，其余7個參數(shù)的取值范圍均在一定的范圍內(nèi)，在此范圍內(nèi)，參數(shù)的不確定性較小，在此范圍外參數(shù)的不確定性增加；以水量平衡系數(shù)作為似然函數(shù)時，只有fc為敏感性參數(shù)，增加了模型的不確性。由于篇幅有限，本文也只展示土壤穩(wěn)定下滲率和蒸散發(fā)補充系數(shù)相關(guān)結(jié)果(圖4)。

3.3.2 模擬結(jié)果的不確定性分析 從不確定性區(qū)間來看(圖5)：在率定期和驗證期模擬的徑流序列區(qū)間涵蓋了44.5%和32.1%的實測值；并不能涵蓋峰值較大的實測徑流量，特別是在驗證期。這可能是由于模型結(jié)構(gòu)的不確定性或者人為選取閾值的原因。

圖4 參數(shù)-水量平衡系數(shù)散點圖與可行參數(shù)組頻率直方圖

圖5 95%置信度下水量平衡系數(shù)作為似然函數(shù)的率定期和驗證期不確定區(qū)間

3.4 不同似然函數(shù)的結(jié)果比較

通過比較NSCE和水量平衡系數(shù)兩個似然函數(shù)對模型參數(shù)敏感性的影響，可以看出：似然函數(shù)選取對模型參數(shù)的敏感性分析影響很大。以NSCE作為似然函數(shù)時，9個與水文有關(guān)的參數(shù)均較敏感，且參數(shù)的最終取值范圍縮小；但以水量平衡系數(shù)作為似然函數(shù)時，只有fc是敏感性參數(shù)。說明參數(shù)的敏感性與似然函數(shù)的選取有很大關(guān)系，以“時間”屬性的NSCE和“水量”屬性的水量平衡系數(shù)由于考慮了不同水文過程，衡量的標準也不同。但是從參數(shù)集中出現(xiàn)的范圍來看，這兩種似然函數(shù)參數(shù)主要集中的范圍相似。值得注意的是，雖然fc在兩種似然函數(shù)判別下均為敏感性參數(shù)，但是其散點圖形狀卻有很大不同，在NSCE作為似然函數(shù)時，fc有非常明顯的波峰和波谷，但是在水量平衡系數(shù)作為似然函數(shù)時，其數(shù)值范圍明顯縮小，呈現(xiàn)類似于“柱狀”的圖形。

通過圖3和5可以看出，無論是在率定期還是驗證期，水量平衡系數(shù)作為似然函數(shù)的模擬結(jié)果要好于NSCE作為似然函數(shù)時的模擬結(jié)果，特別是在率定期流量較大的時候。這可能與水量平衡系數(shù)作為似然函數(shù)時，各參數(shù)取得的可行參數(shù)集取值范圍較大有關(guān)，參數(shù)集能夠充分體現(xiàn)模型模擬的不確定性區(qū)間。

根據(jù)實測值在置信區(qū)間內(nèi)數(shù)量評價模擬結(jié)果的不確定性(表2)。可以看出以NSCE作為似然函數(shù)時，驗證期和率定期的覆蓋率較低，分別為26.2%和15.1%；以水量平衡系數(shù)作為似然函數(shù)時，驗證期和率定期的覆蓋率分別為44.5%和32.1%，說明水量平衡系數(shù)作為似然函數(shù)，對不確定區(qū)間評價的結(jié)果要優(yōu)于NSCE；此外，從平均對稱度S看，NSCE作為似然函數(shù)時，驗證期和率定期的對稱度分別為0.26和0.25；以水量平衡系數(shù)作為似然函數(shù)時，驗證期和率定期的對稱度分別為0.23和0.25，也可證明水量平衡系數(shù)作為似然函數(shù)，選取的不確定性區(qū)間要優(yōu)于NSCE作為似然函數(shù)選取的不確定性區(qū)間。

表2 兩種似然函數(shù)的不確定性區(qū)間評價結(jié)果

4 結(jié) 論

本文以淮河流域上游大坡嶺控制流域為研究區(qū)，分析了GLUE框架下不同似然函數(shù)對模型參數(shù)選取和模擬結(jié)果的影響，主要的結(jié)論如下：

(1)在不同的似然函數(shù)條件下，HEQM模型參數(shù)的取值范圍明顯不同。以NSCE作為似然函數(shù)時，9個參數(shù)的取值范圍縮小較大，其中fc、ketp和Tg3個參數(shù)的取值范圍縮小最大；以水量平衡系數(shù)作為似然函數(shù)時，參數(shù)范圍沒有明顯變化；在不同的似然函數(shù)條件下，參數(shù)對模型的敏感性也不同，以NSCE作為似然函數(shù)時，9個參數(shù)的散點圖均有明顯的波動，因此9個參數(shù)均為敏感性參數(shù)；而以水量平衡系數(shù)作為似然函數(shù)時，只有fc是敏感性參數(shù)。

(2)NSCE和水量平衡系數(shù)均可作為GLUE方法的似然函數(shù)，但從參數(shù)不確定性區(qū)間評價結(jié)果來看，水量平衡系數(shù)作為似然函數(shù)時，覆蓋率和對稱度均好于NSCE作為似然函數(shù)。特別是對于極端徑流，水量平衡系數(shù)作為似然函數(shù)的GLUE方法對參數(shù)的不確定性分析更為準確。

參數(shù)的敏感性和不確定性與似然函數(shù)的選取密切相關(guān)，在今后的研究中可以選取更多具有代表性的似然函數(shù)，來分析似然函數(shù)與參數(shù)之間的關(guān)系，從而為多目標似然函數(shù)的權(quán)重分配提供理論依據(jù)。

[1]SHENZhenyao，CHENLei，CHENTao.AnalysisofparameteruncertaintyinhydrologicalandsedimentmodelingusingGLUEmethod:acasestudyofSWATmodelappliedtoThreeGorgesReservoirRegionChina[J].HydrologyandEarthSystemSciences，2012，16(1)：121-132.

[2] 張質(zhì)明，王曉燕，于 洋，等.基于GLUE法的多指標水質(zhì)模型參數(shù)率定方法[J].環(huán)境科學(xué)學(xué)報，2014，34(7)：1853-1861.

[3] 舒 暢，劉蘇峽，莫興國，等.新安江模型參數(shù)的不確定性分析[J].地理研究，2008，27(2)：343-352.

[4]MCCONNELLP.Hydrologicalforecastinguncertaintyassessment：IncoherenceoftheGLUEmethodology[J].JournalofHydrology，2006,330(1-2)：368-381.

[5]KEITHB，BINLEYA.Thefutureofdistributedmodels：modelcalibrationanduncertaintyprediction[J].HydrologicalProcesses，1992，6(3)：279-298.

[6]BLASONERS，MADSENH，DANR.UncertaintyassessmentofintegrateddistributedhydrologicalmodelsusingGLUEwithMarkovchainMonteCarlosampling[J].JournalofHydrology，2008，353(1-2)：18-32.[7]SELLAMIH，BENABDALLAHS，JEUNIESSEIL，etal.ClimatemodelsandhydrologicalparameteruncertaintiesinclimatechangeimpactsonmonthlyrunoffanddailyflowdurationcurveofaMediterraneancatchment[J].InternationalAssociationofScientificHydrologyBulletin，2016，61(8)：1415-1429.[8]ARABIM，RAOSG，HANTUSHMM.Aprobabilisticapproachforanalysisofuncertaintyintheevaluationofwatershedmanagementpractices[J].JournalofHydrology，2007，333(2-4)：459-471.

[9] 劉艷麗，梁國華，周惠成.水文模型不確定性分析的多準則似然判據(jù)GLUE方法[J].四川大學(xué)學(xué)報(工程科學(xué)版)，2009,41(4)：89-96.

[10]SUNMei，ZHANGXiaoling，HUOZailing，etal.Uncertaintyandsensitivityassessmentsofanagricultural-hydrologicalmodel(RZWQM2)usingtheGLUEmethod[J].JournalofHydrology，2016，534：19-30.

[11] 王綱勝，夏 軍，陳軍鋒.模型多參數(shù)靈敏度與不確定性分析[J].地理研究.2010,29(2)：263-270.[12]ZHANGYongyong，SHAOQuanxi，YEAizhong，etal.Integratedwatersystemsimulationbyconsideringhydrologicalandbiogeochemicalprocesses：modeldevelopment，withparametersensitivityandautocalibration[J].HydrologyandEarthSystemSciencesDiscussions，2016，12(5)：4997-5053.[13]XIAJun.IdentificationofaconstrainednonlinearhydrologicalsystemdescribedbyVolterraFunctionalSeries[J].WaterResourcesResearch，1991，27(9)：2415-2420.[14]HARGREAVESGH，SAMANIZA.Estimatingpotentialevapotranspiration[J].JournaloftheIrrigationandDrainageDivision，1982，108(3):225-230.[15]DUANQingyun，SOROOSHS，VIJAIG.Effectiveandefficientglobaloptimizationforconceptualrainfall‐runoffmodels[J].WaterResourcesResearch，1992，28(4)：1015-1031.[16]XIONGLihua，WANMin，WEIXiaojing，etal.Indicesforassessingthepredictionboundsofhydrologicalmodelsandapplicationbygeneralizedlikelihooduncertaintyestimation[J].InternationalAssociationofScientificHydrologyBulletin，2009，54(5)：852-871.

[17] 周敏敏，瞿思敏，石 朋，等.淮河上游大坡嶺流域土地利用方式變化引起的流域滯時變化[J].河海大學(xué)學(xué)報(自然科學(xué)版)，2015，43(2)：100-106.