張如強，高紅凱，劉俊國,，田 勇

(1.北京林業(yè)大學自然保護區(qū)學院，北京 100083; 2.亞利桑那州立大學，亞利桑那州坦佩 AZ85282； 3.南方科技大學環(huán)境科學與工程學院，廣東 深圳 518055)

水文模型作為模擬水文過程和認識水文規(guī)律的有效工具，已廣泛應用于科學研究領(lǐng)域和解決實際問題，如洪水分析與預報、水資源管理等[1,2]。水文模型是對現(xiàn)實世界的簡化和抽象，模型中包含難以直接測量的參數(shù)，所以需要對模型參數(shù)進行率定，間接得到這些參數(shù)值[3,4]。另外，無論選用哪種水文模型，都需要對模型的模擬效果進行評價[5]。模型的率定和模擬效果的評價通常是對比模擬流量和觀測流量，使目標函數(shù)值最大或最小[6,7]。

模型的模擬效果與目標函數(shù)的選取直接相關(guān)[8]。納什系數(shù)(E)是模型率定過程中應用最為廣泛的目標函數(shù)，也是評價模擬效果的重要數(shù)值指標[7,9,10]。但是，納什系數(shù)存在不足之處：計算公式中觀測流量與模擬流量之差采用平方的形式，這導致洪峰的影響較基流大，即納什系數(shù)對洪峰過于敏感[7]。因此產(chǎn)生了納什系數(shù)的變形形式，主要包括以自然數(shù)為底的對數(shù)形式(Eln)、開方形式(Esqrt)、絕對值形式(Eabs)和相對值形式(Erel)[6,7,11,12]。盡管這些變形形式已經(jīng)被學者應用于模型率定和評價[5,6,13]，但這些目標函數(shù)的選取對模擬效果影響的研究相對較少，選擇哪種變形形式更為合理這一問題尚未得到圓滿解答。

本研究擬通過將納什系數(shù)、各個變形形式設(shè)定為目標函數(shù)，以納什系數(shù)、相關(guān)系數(shù)、平均絕對誤差和平均值之差4個數(shù)值指標綜合分析模型模擬效果，探索不同目標函數(shù)的選取對模擬效果的影響。

1 研究區(qū)域概況及數(shù)據(jù)來源

黑河發(fā)源于祁連山，流經(jīng)甘肅、青海、內(nèi)蒙古自治區(qū)，最后流入居延海，干流全長821 km，是我國西北地區(qū)第二大內(nèi)陸河。以鶯落峽和正義峽為界，黑河流域分為上、中、下游。本文的研究區(qū)為黑河上游，位于祁連山西北，面積約1.000 9 萬km2，海拔高度介于1 700～4 900 m；研究區(qū)覆被包括森林(20%)、草地(52%)、裸地(19%)、濕地或洪泛區(qū)(8%)和永久性冰川雪地(0.8%)。研究區(qū)降雨集中，超過80%大降雨分布在5-9月；多年平均降雨和潛在蒸散發(fā)為430和520 mm/a。

數(shù)據(jù)來源于寒區(qū)旱區(qū)科學數(shù)據(jù)中心。研究使用數(shù)據(jù)為日降雨、日平均氣溫和逐日的徑流量。其中，氣象數(shù)據(jù)為托勒站、祁連站、野牛溝站、張掖站四個氣象站1959-1979年和1990-2006年日降雨和日平均溫度，徑流數(shù)據(jù)為鶯落峽水文站相應時間觀測流量。研究區(qū)位置、覆被以及氣象站水文站位置見圖1。

圖1 研究區(qū)位置、黑河上游覆被及氣象、水文站點分布Fig.1 Location and Land cover of Upper Heihe River basin and distribution of the meteorological andhydrologic stations in and around upper Heihe

2 研究方法

2.1 模型簡介

地形與地質(zhì)條件、土壤類型、氣候、地表覆被等密切相關(guān)，而這些正是影響水文過程的主要因素。地形綜合反映了流域下墊面信息，在不同尺度上影響了降雨-徑流過程[14]?；谶@一假設(shè)Gao等[15]開發(fā)了半分布式水文模型FLEX-Topo，并應用于黑河上游，取得了滿意的模擬效果[15]。因此，本文研究中采用FLEX-Topo模型[15]。依據(jù)地形因素將流域分為4種水文景觀類型。不同水文景觀類型的產(chǎn)匯流特點不同，因此，模型中采用了不同的產(chǎn)匯流機制，使模型具有了半分布式的特點。依據(jù)降雨形成徑流的過程，模型包含了5個水箱：積雪水箱(Sw)、植被冠層截留水箱(Si)、非飽和土壤層水箱(Su)、快速退水水箱(Sf)和慢速退水水箱(Ss)。整個模型通過5個水箱模擬降雨經(jīng)過截留(I)、蒸散發(fā)(Ei&Ea)、優(yōu)先流(RfB)、壤中流(Qf)、地下徑流(Qs)等最終流入河網(wǎng)的過程。模型結(jié)構(gòu)見圖2。FLEX-Topo模型包含24個參數(shù)，參數(shù)含義及取值范圍見表1。

圖2 FLEX-Topo 模型結(jié)構(gòu)[15]Fig.2 The structure of the FLEX-Topo model

參數(shù)名稱取值范圍含義FDD[1，8] 度日模型參數(shù)，表征當溫度高于Tt時，每升高1°每日融化形成的水量Tt[-2．5，2，5] 當溫度高于該參數(shù)時，產(chǎn)生截留Pt[5，35] 當有效降雨大于該參數(shù)時，發(fā)生霍頓超滲產(chǎn)流SumaxB[5，500] 裸地景觀土壤最大持水能力PmaxB[0．1，10] 裸地景觀的最大深層滲漏能力DB[0，1] 裸地景觀中將徑流分成流入快速退水水箱和慢速退水水箱的參數(shù)TlagB[0，5] 裸地景觀中徑流進入快速退水水箱的滯后時間KfB[2，50] 裸地景觀中存儲于快速退水水箱中水流入河網(wǎng)的時間ImaxFH[1，10] 林坡景觀最大截留量SumaxFH[100，1000] 林坡景觀土壤最大持水量βFH[0．1，5] 表征景觀的空間異質(zhì)性，用于計算林坡的徑流系數(shù)D[0，1] 林坡和草坡景觀中將徑流分成流入快速退水水箱和慢速退水水箱的參數(shù)TlagFH[0，5] 林坡景觀徑流進入快速退水水箱的滯后時間KfH[1，20] 林坡和草坡景觀中快速退水水箱中的水流入河網(wǎng)的時間

續(xù)表1 FLEX-Topo模型參數(shù)及參數(shù)取值范圍

模型的輸入數(shù)據(jù)為DEM、日降雨、日平均溫度，輸出為逐日的蒸散發(fā)(Ei&Ea)、徑流量(Qf&Qs)、滲流量(RfB&f)。DEM主要用于生成HAND(Height Above the Nearest Drainage)[15-18]和坡度。黑河上游依據(jù)這兩個因素分為4種水文景觀類型，即裸地(bare soil/rock)、林坡(forest hillslope)、草坡(grass hillslope)、濕地或洪泛區(qū)(wetland/flood plain)，其中，裸地的產(chǎn)流機制主要為深層滲漏補給地下水(DP)和超滲產(chǎn)流(HOF)；林坡和草坡的產(chǎn)流機制為壤中快速徑流(SSF)，主要是優(yōu)先流；濕地或洪泛區(qū)的產(chǎn)流機制為蓄滿坡面產(chǎn)流(SOF)。降雨和溫度是模型的驅(qū)動數(shù)據(jù)，模型采用Hamon方程[19]計算潛在蒸散發(fā)。黑河流域上游面積大、海拔高，而研究區(qū)周邊僅有4個氣象站，所以模型采用泰森多邊形法計算平均降雨量，并同時考慮海拔梯度對溫度和降雨的影響。徑流量是模型率定和驗證的主要數(shù)據(jù)。

2.2 模型率定算法

FLEX-Topo模型采用的參數(shù)率定算法為MOSCEM-UA(Multi-Objective Shuffled Complex Evolution Metropolis-University of Arizona)[20]。Vrugt等[20]開發(fā)了該算法(算法流程圖見圖3)，并將其用于參數(shù)數(shù)量、復雜程度不同的3個模型中，證明該算法的有效性和高效性。該算法已被國內(nèi)外學者廣泛使用[21-25]。MOSCEM-UA算法需要設(shè)置3個參數(shù)：迭代次數(shù)、復合型個數(shù)以及隨機樣本數(shù)。根據(jù)Vrugt等[20]的建議本研究設(shè)定的3個參數(shù)值分別為100 000、24、2 304。

圖3 MOSCEM-UA算法流程圖[20]Fig.3 Flowchart of the MOSCEM-UA algorithm

2.3 研究方法

納什系數(shù)是模型率定過程中應用最為廣泛的目標函數(shù)，鑒于納什系數(shù)的不足之處，產(chǎn)生了變形形式，主要包括以自然數(shù)為底的對數(shù)形式(Eln)、開平方形式(Esqrt)、絕對值形式(Eabs)和相對值形式(Erel)4種。4種變形的計算公式見表2。

FLEX-Topo模型及MOSCEM-UA算法的可靠性已經(jīng)得到證明。本研究采用該模型和算法。首先，使用相同的參數(shù)范圍、迭代次數(shù)(100 000)、復合型個數(shù)(24)以及隨機樣本數(shù)(2 304)，設(shè)定不同的目標函數(shù)(納什系數(shù)及變形形式)進行參數(shù)率定(流程見圖3)，獲得5套參數(shù)。然后，將每一套參數(shù)代入模型進行驗證，獲得5套模擬結(jié)果。最后對5套模擬結(jié)果進行分析、比較。本研究中模型參數(shù)率定期設(shè)為1959-1978年，驗證期設(shè)為1990-2006年，兩階段分別取1959和1990年為模型預熱期。

表2 納什系數(shù)及變形形式計算公式[20]Tab.2 Mathematical formulation of Nash and Sutcliffe efficiency and its modified forms

衡量模擬效果通常有兩種途徑，一是比較直觀的圖形，二是采用數(shù)值指標[5,6]。常用的圖形包括水文過程圖和流量頻率分布曲線，從圖形中可以直觀地查看模擬結(jié)果。其中，水文過程圖強調(diào)流量隨時間的變化，而流量頻率分布曲線側(cè)重于對流量的高估與低估問題；按照有無量綱，數(shù)值指標可以分為相對指標和絕對指標兩大類，其中，相對指標中應用最為廣泛的是納什系數(shù)，絕對指標包括均方根誤差、平均絕對誤差等。Willmott[26]指出判斷模擬效果時納什系數(shù)這一相對指標不宜單獨使用。Legates等[7]通過研究表明，平均絕對誤差略優(yōu)于均方根誤差，因此，根據(jù)Legates等[7]和Anderton等[27]的建議，在評價模型效果時選取納什系數(shù)(E)、相關(guān)系數(shù)(r)、平均絕對誤差(MAE)和平均值之差(DM)4個數(shù)值指標來衡量模擬效果。4個數(shù)值指標的計算公式見表3。

3 結(jié)果與討論

3.1 模擬結(jié)果

根據(jù)模擬結(jié)果繪制水文過程圖和流量頻率分布曲線。不同目標函數(shù)驗證期的模型水文過程圖和流量頻率分布曲線見圖4和圖5。圖4和圖5中，Qo為觀測流量，Qx為以x為目標函數(shù)的模擬流量。

由圖4可以看出，以納什系數(shù)及變形形式為目標函數(shù)的模擬結(jié)果中，模擬流量基本可以反映觀測流量隨時間的變化，但吻合程度不同：較好的是E和Esqrt，Eln和Eabs次之，Erel最差。流量頻率分布曲線反映模型對漲水段、洪峰、退水段和基流等的模擬效果。由圖5可以看出，對于洪峰流量，E出現(xiàn)高估，而Eln、Esqrt和Eabs較E略好，Erel則低估了洪峰流量；對于基流，各個目標函數(shù)的模擬效果相似，均低估了基流流量；漲水段和退水段除Erel最差外，E、Eln、Esqrt和Eabs表現(xiàn)差別不大，均小于觀測值。

表3 3個數(shù)值指標的計算公式Tab.3 Mathematical formulation of the three numerical criterion

圖4 驗證期鶯落峽水文過程圖Fig.4 Hydrograph of validation period of Yingluoxia

圖5 驗證期鶯落峽流量頻率分布曲線Fig.5 Flow duration curve (FDC) of validation period of Yingluoxia

雖然水文過程圖和流量頻率分布曲線可以直觀地看出模擬結(jié)果，但整體模擬效果難以客觀地直接衡量。因此，本文計算了5套模擬結(jié)果的4個數(shù)值指標，計算結(jié)果見表4。

從表4中可以看出：①4種變形形式的模擬結(jié)果中，相關(guān)系數(shù)一列的值相差不大；其余三列除Erel最差外，Eln、Esqrt和Eabs相差很小，Esqrt在各個數(shù)值指標中均為最優(yōu)；②E與Esqrt相比，相關(guān)系數(shù)有所提高，這表明Esqrt較E更能反映流量的動態(tài)變化，而其余指標差別很小；③5個模擬結(jié)果中，Erel的各數(shù)值指標值均為最差；④5個模擬結(jié)果中，平均值之差一列均為正數(shù)，表明模擬結(jié)果均低估了總體水量。

表4 不同目標函數(shù)的數(shù)值指標Tab.4 Numerical criterionof different objective functions

3.2 結(jié)果分析與討論

Legates等[7]指出，納什系數(shù)的計算公式采用平方形式，導致洪峰流量的影響過大，即納什系數(shù)對洪峰流量過于敏感、忽視基流，這是欠妥的。為了解決這一問題，產(chǎn)生了對數(shù)形式、開方形式、絕對值形式、相對值形式。

絕對值形式的納什系數(shù)改變了計算公式中采用的平方形式(見表2)。Legates等[7]建議使用絕對值形式的納什系數(shù)，因為平方形式將模擬流量與觀測流量的差值放大，而絕對值形式將賦予兩者差值相同的權(quán)重，因此減小了洪峰流量的影響。對數(shù)形式、開方形式和相對值形式則采用另一種方式：改變流量的絕對大小(見表2)。因為計算公式中采用這種處理，就降低了洪峰流量的影響。Pushpalatha等[5]、Krause等[6]、Houghton-Carr[11]、Oudin等[12]、Kim等[13]研究證明這些變形降低了納什系數(shù)對洪峰流量的敏感性。董磊華等[28]比較了對數(shù)均方誤和平方根均方誤為目標函數(shù)的模擬結(jié)果，文章指出：由對數(shù)函數(shù)曲線(y=lnx)和開方函數(shù)曲線(y=x1/2)可知，x值越小，y值遞增的梯度越大，同理，流量值越小模擬流量與實測流量的誤差放大越大，因此對數(shù)形式和開方形式削弱了洪峰流量的影響、增強基流流量作用。Krause等[6]提出了基于相對偏差的納什系數(shù)計算方式，即相對值形式。通過這種變形，觀測流量和模擬流量的差值用相對偏差來衡量，從而減小洪峰流量的影響、增強基流流量的影響。但是，從結(jié)果來看，這種變形形式的模擬效果在所有變形形式中最差。

上述研究表明變形形式彌補了納什系數(shù)對洪峰流量過度敏感這一不足。然而，本研究發(fā)現(xiàn)，變形形式對基流流量的模擬并不理想。造成這一結(jié)果的原因可能有兩個方面：一是觀測誤差的存在，包括降雨、溫度、徑流的觀測誤差；二是納什系數(shù)計算公式中參考值為平均流量。對于季節(jié)流量波動較小的流域，選擇平均值為參考值是合適的，但是對于流量波動大的流域則有待商榷。

4 結(jié) 語

本研究通過設(shè)定不同的目標函數(shù)，從納什系數(shù)、相關(guān)系數(shù)、平均絕對誤差和平均值之差4個指標評估模擬效果。對比各個目標函數(shù)驗證期的水文過程圖、流量頻率分布曲線和數(shù)值指標，可以看出以相對值形式為目標函數(shù)的模擬結(jié)果最差，而以納什系數(shù)和開方形式為目標函數(shù)的模擬結(jié)果最好；綜合比較以納什系數(shù)和開方形式為目標函數(shù)的模擬結(jié)果，兩者的數(shù)值指標差別很小，但前者會明顯高估洪峰流量，后者的相關(guān)系數(shù)更高，更能反映流量的動態(tài)變化。整體來講，當選擇納什系數(shù)及其變形形式為目標函數(shù)時，開方形式的模擬效果最好。

盡管如此，5個目標函數(shù)的模擬結(jié)果均低估了基流流量，同時，也低估了整體水量。造成這一結(jié)果的原因可能有兩個方面：一是觀測誤差；二是計算公式中參考值選擇為平均流量。鑒于觀測技術(shù)的限制以及觀測誤差存在的必然性，建議對納什系數(shù)公式中參考值的選擇開展進一步研究以提高模擬效果。

□

[1] 徐宗學,包為民,謝正輝,等. 水文模型[M].北京: 科學出版社, 2009.

[2] 劉金濤, 馮 杰, 張佳寶. 分布式水文模型在流域水資源開發(fā)利用中的應用研究進展[J]. 中國農(nóng)村水利水電，2007,(2):142-144.

[3] Muleta M K. Model performance sensitivity to objective function during automated calibrations [J].Journal of Hydrologic Engineering, 2012,17(6):756-767.

[4] Vrugt J A, Gupta H V, Bouten W, et al. A shuffled complex evolution metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters [J]. Water Resources Research, 2003,39(8):1 201.

[5] Pushpalatha R, Perrin C, Moine N L, et al. A review of efficiency criteria suitable for evaluating low flow simulations [J]. Journal of Hydrology, 2012,420-421:171-182.

[6] Krause P, Boyle D P, B?se F. Comparison of different efficiency criteria for hydrological model assessment [J]. Advances in Geosciences, 2005,5(5):89-97.

[7] Legates D R, Gregory J, Mccabe JR. Evaluating the use of ‘goodness-of-fit’ measures in hydrologic and hydrologic model validation [J].Water Resources Research, 1999,35(1):233-241.

[8] Gan T Y, Dlamini E M, Biftu G F. Effects of model complexity and structure, data quality, and objective functions on hydrologic modelling [J]. Journal of hydrology, 1997,192:81-103.

[9] W?hling T, Samaniego L, Kumar R. Evaluating multiple performance criteria to calibrate the distributed hydrological model of the upper Neckar catchment [J]. Environmental Earth Sciences, 2013,69(2):453-468.

[10] Schaefli B, Gupta H V. Do Nash value have value? [J]. Hydrological Processes, 2007,21:2 075-2 080.

[11] Houghton-carr H A. Assessment criteria for simple conceptual daily rainfall-runoff model [J]. Hydrological Sciences, 1999,44(2):237-261.

[12] Oudin L, Aneréassian V, Mathevet T, et al. Dynamic averaging of rainfall-runoff model simulations from complementary model parameterizations [J]. Water Resources Research, 2006,42(7):887-896.

[13] Kim H S, Lee S. Assessment of a seasonal calibration technique using multiple objectives in rainfall-runoff analysis [J].Hydrological Processes, 2014,28:2 159-2 173.

[14] Savenije H H G. HESS Opinions “Topography driven conceptual modelling (FLEX-Topo)”[J]. Hydrology and Earth System Sciences, 2010,14:2 681-2 692.

[15] Gao H, Hrachowitz M, Fenicia F, et al. Testing the realism of a topography-driven model (FLEX-Topo) in the nested catchments of the Upper Heihe, China[J].Hydrology and Earth System Sciences, 2014,18:1 895-1 915.

[16] Renn? C D, Nobre A D, Cuartas L A, et al. HAND, a new terrain descriptor using SRTM-DEM: Mapping terra-firme rainforest environment in Amazonia [J]. Remote Sensing of Environment, 2008,112:3 469-3 481.

[17] Nobre A D, Cuartas L A, Hodnett M, et al. Height Above the Nearest Drainage-a hydrologically relevant new terrain model [J]. Journal of Hydrology, 2011,404:13-29.

[18] Gharari S, Hrachowitz M, Fenicia F, et al. Hydrological landscape classification: investigating the performance of HAND based landscape classifications in a central European meso-scale catchment [J]. Hydrology and Earth System Sciences, 2011,15:3 275-3 291.

[19] Hamon W R. Estimating potential evapotranspiration [J]. Journal of the Hydraulics Division-Asce,1961,87:107-120.

[20] Vrugt J A, Gupta H V, Bastidas L A, et al. Effective and efficient algorithm for multiobjective optimization of hydrologic model [J]. Water Resources Research, 2003,39(8):1-19.

[21] Schoups G, Addams C L, Gorelick S M. Multi-objective calibration of a surface water-groundwater flow model in an irrigated agricultural region: Yaqui Valley, Sonora, Mexico [J]. Hydrology and Earth System Sciences, 2005,9:549-568.

[22] Tang Y, Reed P, Kollat J. Parallelization strategies for rapid and robust evolutionary multiobjective optimization in water resources applications [J]. Advances in Water Resources,2007,30(3):335-353.

[23] Parajka J, Merz R, Bl?schl G. Uncertainty and multiple objective calibration in regional water balance modelling: case study in 320 Austrian catchments [J]. Hydrological Processes, 2007,21(4):435-446.

[24] Fenicia F, Solomatine D P, Savenije H H G, et al. Soft combination of local models in a multi-objective framework [J]. Hydrology and Earth System Sciences, 2007,11:1 797-1 809.

[25] 王宇暉,雷曉輝,蔣云鐘.基于MOSCEM-UA算法的水文模型參數(shù)優(yōu)化[J].人民長江，2010,41(19):95:99.

[26] Willmott C J. On the validation of models [J]. Physical Geography, 1981,2(2):184-194.

[27] Anderton S, Latron J, Gallart F. Sensitivity analysis and multi-response, multi-criteria evaluation of a physically based distributed model [J]. Hydrological Processes, 2002,16:333-353.

[28] 董磊華, 熊立華. 水文模型中不同目標函數(shù)的影響比較分析[J].水文，2009, 29(3):24-27.