無資料地區(qū)水文模型參數(shù)估算方法比較研究

徐長江，楊無雙，汪青靜，陳 華

(1.長江水利委員會水文局，武漢 430010；2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

0 引 言

隨著科學(xué)技術(shù)的進(jìn)步和人們對于防洪安全意識的提高，人們對水文預(yù)報(bào)的準(zhǔn)確性以及時(shí)空尺度的要求越來越高，這促使了無資料地區(qū)水文模擬的研究越來越受重視。2003年針對發(fā)展中國家的無水文觀測資料或者資料不全的全球重點(diǎn)研究計(jì)劃(PUB計(jì)劃)開展以來，關(guān)于無資料地區(qū)的水文模擬研究無論在理論還是應(yīng)用上都取得了很多成果[1-3]。但由于流域的多樣性、水文機(jī)理的復(fù)雜性、模型的局限性和資料的限制等原因，還有很多難題未解決，所以在未來相當(dāng)長的一段時(shí)間內(nèi)，無資料流域的水文模擬仍然是水文研究的重難點(diǎn)之一[4]。

Parajka[5]等人總結(jié)了PUB計(jì)劃實(shí)施十年間(2003-2013年)的部分研究成果，發(fā)現(xiàn)大多數(shù)研究區(qū)域都集中在北美洲、歐洲和澳大利亞。作為世界上人口最多的發(fā)展中國家，我國人均水資源十分有限，所以水資源規(guī)劃和管理需要更準(zhǔn)確的水文預(yù)報(bào)，另外，地域遼闊、氣候條件多變的特性也使得在許多無資料流域以及PUB科學(xué)問題方面的研究需要探索水文預(yù)報(bào)新的理論與方法[4, 6]，因此在中國開展無資料地區(qū)水文預(yù)報(bào)的研究是十分必要的。本文在PUB計(jì)劃研究成果的基礎(chǔ)上，將4種無資料地區(qū)水文模型參數(shù)估算方法應(yīng)用到湘江流域，驗(yàn)證各方法的可行性，為其在中國無資料地區(qū)的應(yīng)用提供依據(jù)與參考。

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域

湘江流域?yàn)殚L江流域洞庭湖水系，地處北緯24°31′～29°00′，東經(jīng)110°30′～114°00′之間，位于亞熱帶季風(fēng)氣候區(qū)，降水充足，但年內(nèi)分布不均，主要集中在3-7月間。流域水系發(fā)達(dá)，河網(wǎng)密布，干流總長度約為856 km，流域面積為9.48 萬km2。本文根據(jù)實(shí)測水文資料以及流域地形特征，利用湘江流域DEM圖(數(shù)字高程模型)劃分子流域，從中選取了11個(gè)子流域?yàn)檠芯繀^(qū)域。

1.2 數(shù)據(jù)與處理

1.2.1 氣象數(shù)據(jù)

本文研究所用氣象數(shù)據(jù)是洞庭湖區(qū)水文年鑒中摘錄的2006-2014年共9年的日降水、日蒸發(fā)和日徑流實(shí)測數(shù)據(jù)。其中，流域內(nèi)雨量站有274個(gè)站點(diǎn)，蒸發(fā)站有14個(gè)站點(diǎn)，考慮到有些雨量站和蒸發(fā)站分布在流域邊界上或流域外且分布不均勻，因此利用泰森多邊形法對各子流域的面雨量、面蒸發(fā)量進(jìn)行了整理計(jì)算[7]，結(jié)果如表1所示。

表1 流域氣象特征值Tab.1 Meteorological Characteristics in Xiangjiang River Basin

1.2.2 流域下墊面特征值

流域下墊面特征值是指一個(gè)流域的坡度、植被、土壤、地形地貌、流域形狀等物理性質(zhì)，決定著流域的蓄水量、透水性等特性的不同，因而對流域產(chǎn)匯流存在著很大的影響。本文根據(jù)現(xiàn)有資料共選取了8種下墊面特征值，如表2所示。

2 模型與方法

2.1 新安江模型

新安江模型是在1973年由河海大學(xué)趙人俊教授領(lǐng)導(dǎo)的研究組設(shè)計(jì)的國內(nèi)第一個(gè)完整的流域水文模型。該模型在國內(nèi)外濕潤半濕潤地區(qū)的模擬效果普遍較好，應(yīng)用十分廣泛[9-11]。本文研究使用三水源新安江模型，共15個(gè)參數(shù)，每個(gè)參數(shù)都具有明確的物理意義[12,13]。

表2 流域下墊面特征值 Tab.2 Underlying surface characteristics in Xiangjiang River Basin

2.1.1 SCE-UA優(yōu)化算法

SCE-UA(Shuffled Complex Evolution)算法是一種非常有效的解決概念性水文模型參數(shù)率定問題的全局優(yōu)化算法[14-16]。該算法主要基于信息共享和自然生物進(jìn)化的理念[17]，整合了全局抽樣以及復(fù)合形進(jìn)化等概念的優(yōu)勢[18]，這些特性能夠確保充分利用樣本信息，大大提高了算法的收斂效率[16, 19]，因此被廣泛地應(yīng)用于概念性水文模型參數(shù)的率定問題[20-22]。SCE-UA優(yōu)化算法的參數(shù)比較多，但多數(shù)可以根據(jù)已有研究結(jié)果確定，僅有復(fù)合型個(gè)數(shù)p需要根據(jù)實(shí)際問題確定[17]。p值越大，收斂于全局平均的概率越大，但計(jì)算量也越大，反之亦然[23]。本文取經(jīng)驗(yàn)值p=2。

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

本文采用納西效率系數(shù)(NSE)和水量相對誤差(RE)作為模型模擬精度的評價(jià)指標(biāo)。

(1)

(2)

式中：Qsim,i為模擬流量序列；Qobs,i為實(shí)測流量序列。

2.2 Sobol’敏感性分析方法

由于模型參數(shù)較多，一些不敏感參數(shù)對水文模擬的結(jié)果影響較小，本文區(qū)分出較敏感參數(shù)和不敏感參數(shù)，只對敏感參數(shù)用各移植方法進(jìn)行移植，可以減少參數(shù)間的非獨(dú)立性，從而降低模型的不確定性。Sobol’是一種非常有效的基于方差原理的全局敏感性定量分析方法[24]。該方法不僅能定量描述單個(gè)獨(dú)立參數(shù)的敏感性，而且考慮到不同參數(shù)間相互作用對模型輸出產(chǎn)生影響，能有效地分析非線性模型中參數(shù)之間相互作用的敏感性[25]，較局部敏感性分析更接近實(shí)際，因此近年來受到了廣泛的應(yīng)用[25-27]。

Sobol’方法首先將參數(shù)抽樣樣本轉(zhuǎn)化到一個(gè)多維單元體中，將模型目標(biāo)函數(shù)分解為單個(gè)模型參數(shù)與參數(shù)間相互作用的子函數(shù)之和，根據(jù)方差分解原理，模型總方差也可以分解為單個(gè)模型參數(shù)與參數(shù)間相互作用的方差之和，通過蒙特卡洛擬合求解方法對目標(biāo)函數(shù)的總方差以及各子項(xiàng)的偏方差進(jìn)行估算，最后根據(jù)各子項(xiàng)偏方差對總方差的貢獻(xiàn)比例來表示模型參數(shù)及其相互作用的敏感程度[28]。

抽樣方法上，本研究采用的是拉丁超立方抽樣方法，與傳統(tǒng)的蒙特卡洛抽樣方法相比，該方法能夠使產(chǎn)生的隨機(jī)點(diǎn)更均勻地分布在采樣區(qū)域內(nèi)，能夠有效地用采樣點(diǎn)反映隨機(jī)變量的整體分布[29, 30]。

2.3 參數(shù)移植方法

目前對無資料地區(qū)徑流預(yù)報(bào)常用的方法為區(qū)域化方法。該方法是在一定范圍內(nèi)，將有資料流域的水文信息比如水文模型參數(shù)通過一定的方法移植到無資料流域，從而實(shí)現(xiàn)對無資料流域的預(yù)報(bào)，由于區(qū)域化方法能夠利用流域下墊面特征值等信息，所以能夠有效地降低不確定性影響，從而提高預(yù)報(bào)精度。常用的區(qū)域化方法有多元回歸法、空間鄰近法、流域物理特性相似法以及全局平均法[31]。

(1)多元回歸法：在有資料流域，以新安江模型參數(shù)為因變量，以流域下墊面物理特征值為自變量建立多元回歸方程，根據(jù)多元回歸方程和無資料地區(qū)的下墊面特征值反推出模型參數(shù)值，進(jìn)而達(dá)到參數(shù)移植的目的。本文通過逐步回歸法建立敏感參數(shù)與流域下墊面特征值之間的關(guān)系。

(2)空間鄰近法：直接以地理上相近的有資料流域的模型參數(shù)作為本流域的模型參數(shù)。本文中兩個(gè)流域間的距離利用流域質(zhì)心間的經(jīng)緯度進(jìn)行計(jì)算。

(3)流域物理特性相似法：通過比較各流域的相似性，以與無資料流域最相似流域的模型參數(shù)作為無資料流域的模型參數(shù)。本文采用下式計(jì)算流域相似度，Z的范圍為0～1，越小表示兩個(gè)流域越相似[32]：

(3)

式中：XG和XU分別為有資料流域和無資料流域的下墊面特征值；ΔX為某屬性在各流域中最大值和最小值的差值；k為下墊面特征值個(gè)數(shù)。

(4)全局平均法：將研究區(qū)域內(nèi)，各有資料流域模型參數(shù)的平均值移植到無資料流域進(jìn)而進(jìn)行水文模擬，本文中該方法僅針對兩個(gè)較大流域使用。本文對選取的11個(gè)子流域，依次假定其中一個(gè)流域?yàn)椤盁o資料”流域，其余為有資料流域，通過以上區(qū)域化方法將有資料流域的模型敏感參數(shù)移植到“無資料”流域，其余不敏感參數(shù)則用各有資料流域模型參數(shù)的平均值代替，檢驗(yàn)各方法的可行性及其優(yōu)劣。

3 結(jié)果與討論

3.1 參數(shù)敏感性分析結(jié)果

不同目標(biāo)函數(shù)側(cè)重不同的徑流特性，因此不同目標(biāo)函數(shù)下的參數(shù)敏感性也會有所差異[27]，本文以納西效率系數(shù)(NSE)和水量相對誤差(RE)作為目標(biāo)函數(shù)對新安江模型參數(shù)敏感性進(jìn)行分析(圖1)。

圖1 不同流域各參數(shù)的一階敏感度和總敏感度Fig.1 The first-order sensitivity and total sensitivity of each parameter in different watersheds

圖1以大西灘和湘潭流域?yàn)槔嫵隽瞬煌繕?biāo)函數(shù)下各參數(shù)的一階敏感度和總敏感度，文中指定當(dāng)一階敏感度超過閾值(10%)時(shí)，該參數(shù)為敏感參數(shù)[27]。結(jié)果顯示，當(dāng)目標(biāo)函數(shù)為NSE時(shí)，主要的敏感參數(shù)有KG、CG、N、NK，其中CG控制著地下水的徑流過程，對洪水過程線影響較大，因此較為敏感；而KG控制著地下水占比；N和NK決定著地面匯流單位線的形狀和大小，對洪峰流量的模擬都有很大影響，因此較敏感。當(dāng)目標(biāo)函數(shù)為RE時(shí)，主要的敏感參數(shù)有CG和KE，其中KE控制著水量平衡，是影響產(chǎn)流量最為重要的參數(shù)，對水量計(jì)算非常重要，因此具有較強(qiáng)的敏感性；而CG控制著地下退水過程的快慢，對最終的流量過程有很大影響，因此也較敏感。表3匯總了各流域敏感性參數(shù)。

表3 各流域敏感性參數(shù)Tab.3 Sensitivity parameters of each Sub-basin

3.2 參數(shù)移植結(jié)果

3.2.1 模型率定

建立基于SCE-UA全局優(yōu)化算法的三水源新安江模型，參數(shù)率定結(jié)果見表4。從表4中可以看出， 11個(gè)子流域的納西效率系數(shù)均能達(dá)到0.75以上，說明模擬結(jié)果較好。另以大西灘站和衡陽站為例畫出了2013年兩個(gè)站的日流量過程示意圖，圖示結(jié)果也表明新安江模型能夠較好地模擬大西灘站和衡陽站的流量過程，如圖2和圖3所示。以上結(jié)果表明新安江模型能夠較準(zhǔn)確地模擬湘江流域的降雨徑流關(guān)系，可用于無資料流域水文模型參數(shù)估算方法的研究。

表4 各子流域模擬結(jié)果Tab.4 Simulation results of each sub-basin

圖2 大西灘站新安江模型2013年日流量過程線Fig.2 The 2013 daily flow hydrograph by Xin'anjiang Model at Daxitan Station

圖3 衡陽站新安江模型2013年日流量過程線Fig.3 The 2013 daily flow hydrograph by Xin'anjiang Model at Hengyang Station

3.2.2 區(qū)域化方法移植結(jié)果

將建立的區(qū)域化方法應(yīng)用到無資料流域，結(jié)果如表5以及圖4所示。可以看出，通過參數(shù)移植，各“無資料”流域的納西效率系數(shù)均在0.65以上，說明各區(qū)域化方法在湘江流域是適用的，且能得到較好的結(jié)果。另外對比各流域之間的結(jié)果可以看出，各區(qū)域化方法在不同無資料地區(qū)的表現(xiàn)存在差異，如在朗梨、射埠和下馬渡3個(gè)子流域內(nèi)，多元回歸法應(yīng)用結(jié)果最好；而在全州流域，結(jié)果正好相反；另外在湘鄉(xiāng)、神山頭、歐陽海、衡陽和湘潭5個(gè)子流域內(nèi)，多元回歸法和物理特性相似法的結(jié)果相似，且優(yōu)于空間鄰近法和全局平均法；而在大西灘流域，3種區(qū)域化方法的結(jié)果大致相同，且和實(shí)測資料模擬結(jié)果很接近。

表5 參數(shù)移植結(jié)果Tab.5 The Results of Parameter Transplant

圖4 衡陽站參數(shù)移植2013年日流量過程線Fig.4 The 2013 daily flow hydrograph by parameter transplant at hengyang station

4 結(jié) 語

本文在總結(jié)PUB計(jì)劃研究成果的基礎(chǔ)上，研究了空間鄰近法、物理特性相似法、多元回歸法和全局平均法4種無資料地區(qū)參數(shù)移植方法在湘江流域的可行性。結(jié)果表明，各區(qū)域化方法在湘江流域上是適用的，能夠較好地模擬出湘江流域內(nèi)不同子流域的降雨徑流關(guān)系，但在不同的流域上表現(xiàn)各異，優(yōu)劣不同。但就總體來看，多元回歸法的移植效果普遍較好，空間鄰近法和流域物理特性相似法次之，全局平均法較差。此外，雖然區(qū)域化方法結(jié)果較實(shí)測資料模擬結(jié)果的精度相對偏低，但對無資料地區(qū)，通過這些方法進(jìn)行水文預(yù)報(bào)仍不失為一種有效的解決方法，并可為其在中國無資料地區(qū)的研究和應(yīng)用提供參考。

□

參考文獻(xiàn)：

[1] Bao Z, Zhang J, Liu J, et al. Comparison of regionalization approaches based on regression and similarity for predictions in ungauged catchments under multiple hydro-climatic conditions [J]. Journal of Hydrology, 2012, 466:37-46.

[2] Samuel J, Coulibaly P, Metcalfe R A. Estimation of Continuous Streamflow in Ontario Ungauged Basins: Comparison of Regionalization Methods [J]. Journal of Hydrologic Engineering, 2011,16(5):447-459.

[3] Zhang Y, Chiew F H S. Relative merits of different methods for runoff predictions in ungauged catchments [J]. Water Resources Research, 2009,45.

[4] 楊大文, 劉志雨, 張建云, 等. 中國的PUB科學(xué)問題和研究進(jìn)展[C]∥ 全國第三屆水問題研究學(xué)術(shù)研討會. 西安, 2005:39-45.

[5] Parajka J, Viglione A, Rogger M, et al. Comparative assessment of predictions in ungauged basins-Part 1: Runoff-hydrograph studies [J]. Hydrology and Earth System Sciences, 2013,17(5):1 783-1 795.

[6] 楊大文, 夏 軍, 張建云. 中國PUB研究與發(fā)展[C]∥ 第二屆全國水問題研究學(xué)術(shù)研討會. 北京, 2004:47-54.

[7] 徐 晶, 林 建, 姚學(xué)祥, 等. 七大江河流域面雨量計(jì)算方法及應(yīng)用[J]. 氣象，2001,27(11):13-16,51.

[8] 施 征, 包為民, 瞿思敏. 基于相似性的無資料地區(qū)模型參數(shù)確定[J]. 水文，2015,35(2):33-38.

[9] Lu H, Hou T, Horton R, et al. The streamflow estimation using the Xinanjiang rainfall runoff model and dual state-parameter estimation method [J]. Journal of Hydrology, 2013,480:102-114.

[10] 鄧元倩, 李致家, 劉甲奇, 等. 基于SCE-UA算法新安江模型在灃河流域的應(yīng)用[J]. 水資源與水工程學(xué)報(bào)，2017,28(3):27-31.

[11] 朱求安, 張萬昌. 新安江模型在漢江江口流域的應(yīng)用及適應(yīng)性分析[J]. 水資源與水工程學(xué)報(bào)，2004,15(3):19-23.

[12] 王佩蘭, 趙人俊. 新安江模型(三水源)參數(shù)的客觀優(yōu)選方法[J]. 河海大學(xué)學(xué)報(bào)，1989,(4):65-69.

[13] 趙人俊, 王佩蘭. 新安江模型參數(shù)的分析[J]. 水文，1988,(6):2-9.

[14] Wang Y, Yu P, Yang T. Comparison of genetic algorithms and shuffled complex evolution approach for calibrating distributed rainfall-runoff model [J]. Hydrological Processes, 2010,24(8):1 015-1 026.

[15] 宋星原, 舒全英, 王海波, 等. SCE-UA、遺傳算法和單純形優(yōu)化算法的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)(工學(xué)版)，2009,42(1):6-9,15.

[16] Duan Q Y, Gupta V K, Sorooshian S. Shuffled complex evolution approach for effective and efficient global minimization[J]. Journal of Optimization Theory and Applications, 1993,76(3):501-521.

[17] Duan Q Y, Sorooshian S, Gupta V K. Optimal use of the sce-ua global optimization method for calibrating watershed models [J]. Journal of Hydrology, 1994,158(3-4):265-284.

[18] Nelder J A, Mead R. A simplex-method for function minimization [J]. Computer Journal, 1965,7(4):308-313.

[19] Jeon J, Park C, Engel B A. Comparison of Performance between Genetic Algorithm and SCE-UA for Calibration of SCS-CN Surface Runoff Simulation [J]. Water, 2014,6(11):3 433-3 456.

[20] 馬海波, 董增川, 張文明, 等. SCE-UA算法在TOPMODEL參數(shù)優(yōu)化中的應(yīng)用[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2006,(4):361-365.

[21] 唐運(yùn)憶, 欒承梅. SCE-UA算法在新安江模型及TOPMODEL參數(shù)優(yōu)化應(yīng)用中的研究[J]. 水文，2007,27(6):33-35, 74.

[22] 董潔平, 李致家, 戴健男. 基于SCE-UA算法的新安江模型參數(shù)優(yōu)化及應(yīng)用[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2012,40(5):485-490.

[23] Madsen H. Automatic calibration of a conceptual rainfall-runoff model using multiple objectives[J]. Journal of Hydrology, 2000,235(3-4):276-288.

[24] Sobol I M. Sensitivity Estimates for Nonlinear Mathematical Models [J]. Mathematical Modelling and Computational Experiment, 1993,1(4):407-414.

[25] 齊 偉, 張 弛, 初京剛, 等. Sobol方法分析TOPMODEL水文模型參數(shù)敏感性[J]. 水文, 2014,34(2):49-54.

[26] Yatheendradas S, Wagener T, Gupta H, et al. Understanding uncertainty in distributed flash flood forecasting for semiarid regions [J]. Water Resources Research, 2008,44.

[27] 張小麗, 彭 勇, 徐 煒, 等. 基于 Sobol 方法的新安江模型參數(shù)敏感性分析[J]. 南水北調(diào)與水利科技，2014,(2):20-24,33.

[28] Sobol I M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates [J]. Mathematics and Computers in Simulation, 2001,55(1-3):271-280.

[29] Stein M. Large Sample Properties of Simulations Using Latin Hypercube Sampling [J]. Technometrics, 1987,29(2):143.

[30] Saltelli A. Making best use of model evaluations to compute sensitivity indices [J]. Computer Physics Communications, 2002, 145(PII S0010-4655(02)00280-12):280-297.

[31] 許崇育, 陳 華, 郭生練. 變化環(huán)境下水文模擬的幾個(gè)關(guān)鍵問題和挑戰(zhàn)[J]. Journal of Water Resources Research, 2013,(2):85-95.

[32] Parajka J, Merz R, Bloschl G. A comparison of regionalisation methods for catchment model parameters [J]. Hydrology and Earth System Sciences, 2005,9(3):157-171.