基于貝葉斯模型平均的水文模型不確定性及集合模擬

江善虎，任立良，劉淑雅，袁 飛 ，楊肖麗

(1 河海大學水文水資源與水利工程科學國家重點實驗室，南京 210098；2 河海大學水文水資源學院，南京 210098)

水文模型是對自然界復雜水循環過程的近似描述，是水文科學研究的重要手段和方法[1]。自1932年Sherman提出著名的水文單位線概念和1933年Horton提出著名的下滲定理至今，各國水文工作者在水文循環規律和水文過程模擬方面做了大量研究工作，獲得了豐碩的成果；特別是20世紀70年代以來，國內外水文學者研制了大量不同的水文模型，并在現代洪水預報、氣候變化與人類活動的水文響應、流域生態水文過程模擬和水資源規劃與管理等領域廣泛應用[2,3]。然而水文過程是一種受氣候、氣象、地形、地貌和下墊面等較多因素影響的復雜過程，蘊含著確定性的動態規律和不確定性的統計規律[4]。由于當前人類認知水平的有限，使得對水文過程的認識存在著一定的不完全性，而據此建立的水文模型也存在著諸多不確定性[5]?，F階段，根據水文模型進行流域水文過程模擬及預報時不可避免地存在著模型輸入、模型參數和模型結構的不確定性[6-11]。針對模型輸入(主要針對降水)的不確定性， Kavetski等[7]、Ajami等[8]引入一個服從正態分布的隨機乘子來定量考慮。水文模型參數的估計方法可分為兩類：第一類估計方法認為在給定流域和模型結構時，水文模型存在唯一的一組最優參數解；第二類估計方法認為在給定流域和模型結構時，水文模型的參數是一組隨機變量，它們服從一個確定的聯合概率分布[6, 9]。目前，已有大量的水文模型參數優化和不確定性分析方法，其中具有代表性的有第一類的遺傳算法(GA)、SCE-UA算法等和第二類的基于貝葉斯統計理論的GLUE方法、SCEM-UA算法等。對于模型結構的不確定性，由于不同模型具有不同的優缺點，單一模型的模擬效果總是有限的，需充分發揮每個模型的優勢，對每個模型分別考慮參數和輸入的不確定性，進而對不同模型的模擬結果進行集合來消減模型結構的不確定性[10-13]。因此，綜合考慮模型輸入、模型參數和模型結構的不確定性，進行基于多模型的流域水文過程集合模擬，在給定確定模擬值的同時，又提供概率預報區間，是當前減少水文模擬不確定性的有效方法[14-16]。

本文選用基于子流域的新安江模型 (XAJ)[17]、混合產流模型 (HYB)[18]和改進的HYMOD模型 (HYM)[19]，從模型參數估計的兩類方法出發，分別采用單純多邊形進化算法(SCE-UA)[20]和SCEM-UA算法[21]來進行模型參數優化和不確定性分析，探討兩類方法對水文模型模擬結果的影響；同時引入輸入誤差乘子來進行模型輸入的不確定性分析；最終采用貝葉斯模型平均(BMA)將3個模型的模擬結果集合進行概率預報。

1 研究方法

1.1 水文模型

選用基于子流域的新安江模型、混合產流模型和改進的HYMOD模型進行水文過程模擬。產流機制方面，新安江模型和HYMOD模型為蓄滿產流，混合產流模型考慮了蓄滿和超滲兩種產流機制的垂向組合。蒸散發計算，新安江模型和混合產流模型均考慮三層蒸散發，HYMOD模型僅考慮一層蒸散發。徑流劃分方面混合產流模型和HYMOD模型為二水源模型，新安江模型為三水源模型。匯流方面，把全流域劃分為15個單元子流域，坡面匯流采用線性水庫進行計算，河道匯流采用馬斯京根分段連續演算方法。新安江模型有16個參數，混合產流模型為14個參數，改進的HYMOD模型有8個參數。3個模型中，HYMOD模型最為簡單，混合產流模型為在新安江模型基礎上改進而來，產流機制方面較新安江模型更為復雜；模型參數方面，新安江模型和混合產流模型參數相當，HYMOD模型參數最少。3個模型的驅動數據都較為簡單，僅為降水和蒸發。表1列出了這3個模型的主要特征。

表1 3個水文模型主要特征比較Tab.1 The main characteristics of the three hydrological models

1.2 輸入不確定性

水文模型的輸入數據是水文氣象要素，其中最重要的是降雨，其不確定性對水文過程模擬結果具有至關重要的影響。本研究采用一個誤差乘子來定量表述降雨輸入的不確定性。假定降水輸入誤差服從正態分布，用隨機乘數來修正實測降水量，每個實測降水量都對應一個乘數：

(1)

(2)

1.3 參數不確定性

采用SCE-UA算法和SCEM-UA算法來進行模型的參數優化和不確定性分析。SCE-UA是一種全局優化算法，該算法結合了現有算法(包括基因算法等)的一些優點，將基于確定性的復合型搜索技術和自然界中的生物競爭進化原理相結合，可以解決高維參數的全局優化問題，且不需要顯示目標函數或目標函數的偏導數[19]。SCE-UA算法可以有效、快速地搜索到水文模型參數的全局最優解，在水文模型參數優化中得到了廣泛的應用。SCEM-UA算法是Vrugt等于2003年對單純多邊形進化算法(SCE-UA)加以改進而提出的一種全局搜索算法，該算法不僅能搜索到一套具有最大似然性的參數組，同時也能估計參數不確定性[21]。SCEM-UA算法結合了馬爾科夫蒙特卡洛(MCMC)采樣技術，以Metropolis-Hastings算法(簡稱M-H算法)取代SCE-UA算法中的坡降算法，并估計出最優參數組及其后驗概率分布[22]。

1.4 貝葉斯模型平均

BMA是一種基于貝葉斯理論的多模型集合分析方法，它基于似然函數推斷單個模型的權重，保證模擬效果更好的模型具有更高的權重系數；通過多個水文模型競爭，BMA可產生更可靠的集合預報值。BMA不僅可以用于模型集合，而且可以用于計算單個模型和模型集合的不確定性。其基本原理簡介如下。

假設y為BMA集合預報流量，D=[X,Y]為實測數據(其中X表示輸入資料，Y表示實測的流量資料)，f=[f1,f2,…,fk]為k個模型預報的集合。BMA的概率預報表示如下：

(3)

(4)

(5)

1.5 不確定性區間估算

(3) 重復步驟(1)和(2)M次。M是在任意時刻t的樣本容量，本文中令M=1 000。將它們從小到大排序，BMA的95%預報區間就是2.5%和97.5%分位數之間的部分。

針對SCEM-UA算法優化模型的模擬結果，將模型參數收斂前的模擬結果和收斂后的模擬結果綜合起來考慮，因為樣本容量足夠大(本文為15 000組)，可以近似服從正態分布，根據正態總體均值的區間估計方法可計算出95%置信水平下樣本的置信區間[8]。

1.6 評估指標

采用Nash-Sutcliffe效率系數(NS)、偏差(bias,B)和均方根誤差(Root mean square error,RM)來評價模型模擬徑流的精度。其中NS系數表征模擬徑流和觀測徑流的吻合程度；B反映模擬徑流的系統偏差程度；RM表征模擬徑流和觀測徑流的誤差大小。采用覆蓋率(Containing ratio,CR)、平均帶寬(Average Bandwidth,BW)和平均偏移度(Average Deviation Amplitude,DA)來評價模型預報區間的準確度[23]。CR指預報區間覆蓋實測流量數據的比例，CR越大，表示預報區間的覆蓋率越高。BW表示預報區間的平均帶寬，對于一定的置信水平，在保證有較高的覆蓋率前提下，BW越窄越好。DA反映預報區間的中心線偏離實測流量過程線的程度，理論上DA越小，表示預報區間的對稱性越好。

2 研究數據

2.1 研究區域概況

洣水流域地處湖南省的東南部(東經112°52′～113°56′ ，北緯26°01′～27°10′)，為湘江的一級支流，其中甘溪水文站以上控制面積為9 977 km2(圖1)。流域地勢東南高、西北低，海拔高程自東南向西北由2 093 m降低到49 m，平均坡降為1.01‰。流域氣候類型屬亞熱帶季風型濕潤性，四季分明，降水充沛，多年平均年降水量達1 561.0 mm。流域多年平均氣溫18.0 ℃，多年平均風速1.8 m/s。

圖1 洣水流域地理位置及測站分布圖Fig.1 Location of Mishui basin and distribution of measurement stations

2.2 研究數據

收集洣水流域35個雨量站、1個蒸發站(五里牌，E601蒸發皿)和1個流量站(甘溪水文站)2000-2008年逐日降雨、潛在蒸發和流量觀測數據。采用臨近三站點反距離加權法(IDW)將點雨量插值到流域面上分布。流域地形資料采用美國地質調查局(USGS)的全球30″數字高程模型數據(DEM)。應用ArcGIS提取構建分布式流域水文模型所需的基本空間信息，包括流域邊界、柵格水流流向、各柵格至流域出口斷面的河長等，并將流域劃分為15個子流域。土地覆被資料采用美國馬里蘭大學全球1 km土地覆被資料，統計得洣水流域主要植被類型為林地占61.8%，草地占24.6%，耕地占13.3%，城鄉和水體占0.3%。

3 結果對比

利用洣水流域2000-2008年連續9年的日降雨-徑流資料進行水文過程模擬和不確定性分析，其中2000-2005年為模型的率定期，2006-2008年為模型的驗證期。SCE-UA算法中，采用RM指標作為目標函數來衡量優選參數的優劣，反映模型模擬值和實測值之間的誤差。SCEM-UA算法中，初始均勻分布樣本個數為1000，復合型個數為5，判斷后驗分布推斷完畢的樣本個數為5000，最大運行次數為15000。本研究采用多種情景進行水文過程模擬和不確定性分析：① 僅考慮模型參數的不確定性，采用SCE-UA算法和SCEM-UA算法分別進行單個模型的參數優選和不確定性分析；② 引入一個誤差乘子來定量表述降雨輸入的不確定性，對每一模型同時考慮模型參數和模型輸入的不確定性；③ 采用BMA將3個模型的模擬結果集合，綜合考慮模型參數、模型輸入和模型結構的不確定性。具體結果分析如下。

3.1 BMA集合和單個模型的模擬值對比

表2為單個模型只考慮模型參數不確定性，單個模型考慮模型參數和模型輸入不確定性，以及利用BMA法綜合考慮模型參數、模型輸入和模型結構不確定性的徑流模擬精度比較(NS、B、RM)。圖2和3分別為BMA法集合SCE-UA優化和SCEM-UA優化的3個模型的模擬結果。從表2和圖2，圖3中可以看出：基于SCE-UA和SCEM-UA算法優選參數，3個模型在洣水流域都能取得較好的模擬結果，率定期模型的NS都在0.85以上，驗證期模型的NS都在0.80以上；模型在率定期的模擬精度要優于驗證期的模擬精度，說明模型模擬結果在驗證期存在更大的不確定性。3個模型考慮輸入的不確定性后，模擬精度和不考慮輸入不確定時相當，沒有大的提高(僅驗證期新安江模型和混合產流模型的NS系數有所提高)，這可能是由于：① 洣水流域雨量站點較為密集，降水輸入較為準確；② 3個模型中均存在降水蒸散發折減系數參數，在一定程度上描述了模型輸入的不確定性；③ 考慮模型輸入的不確定性，引入了兩個新的模型參數，一定程度上增加了模型參數的不確定性。不管是SCE-UA算法，還是SCEM-UA算法，采用BMA法綜合考慮模型參數、模型輸入和模型結構不確定性的模擬結果都要優于單個模型，在率定期和驗證期，BMA都能取得最高的NS(SCE-UA優化的模型BMA集合后，率定期NS為0.91，驗證期為0.88；SCEM-UA優化的模型BMA集合后，率定期NS為0.92，驗證期為0.87)、最小的B和最小的RM，尤其在驗證期BMA集合的NS要較單個模型有很大提高，說明BMA集合多模型提高預報精度的有效性。

表2 不同模擬情景下SCE-UA和SCEM-UA優化的水文模型模擬精度比較Tab.2 Precision performance of the streamflow simulation series at different simulation cases

注：參數表示僅考慮參數不確定性；參數+輸入表示同時考慮參數和輸入不確定性；參數+輸入+結構表示綜合考慮模型參數、輸入和結構不確定性；黑體字表示最優值。

3.2 BMA集合和單個模型的預報區間對比

表3給出了不同方案下單個模型和BMA集合模擬徑流的不確定性區間優良性評價指標(CR、BW、DA)。圖2和3分別展示了BMA法集合SCE-UA優化和SCEM-UA優化的3個模型模擬徑流的95%置信區間。從表2和圖2，圖3中可以看出：不管是SCE-UA算法，還是SCEM-UA算法，采用BMA方法綜合考慮模型參數、模型輸入和模型結構不確定性的預報不確定性區間要優于單個模型，BMA的模擬徑流不確定性區間具有更高的覆蓋率和更優的平均偏移幅度，但它較單個模型具有更大的平均帶寬；這與Xiong等[23]和董磊華等[10]的研究結論一致。橫向比較表3中SCE-UA算法和SCEM-UA算法的模擬結果，可以發現，SCEM-UA算法的概率預報性區間具有更高的CR、更小的BW和更優的DA；其中采用BMA方法集合SCEM-UA算法優化的3個模型模擬結果，其率定期95%置信區間的CR達到95.62%，驗證期95%置信區間的CR達到95.17%，幾乎包含了整個實測徑流序列，取得了最優的結果?？傮w上，采用BMA方法綜合考慮模型參數、模型輸入和模型結構不確定性，較單個模型能夠提供更可靠的預報不確定性區間；與SCE-UA算法相比，SCEM-UA算法能更好進行參數的不確定性分析，基于SCEM-UA算法和BMA集合得到的概率預報區間更為準確，具有更大優勢。

圖2 貝葉斯模型平均集合SCE-UA優化的3個模型的模擬結果Fig.2 The BMA mean from the SCE-UA based simulations of the three hydrological models, and the 95% confidence interval

圖3 貝葉斯模型平均集合SCEM-UA優化的3個模型的模擬結果Fig.3 The BMA mean from the SCEM-UA based simulations of the three hydrological models, and the 95% confidence interval

不同情景SCEUA(抽樣1000次)CR/%RM/(m3·s)DA/(m3·s)SCEMUACR/%RM/(m3·s)DA/(m3·s)率定期XAJ(參數)59.31152.8758.2078.65169.1752.31XAJ(參數+輸入)74.86200.1560.7879.06169.7851.79HYB(參數)75.05258.6674.0880.34222.2564.41HYB(參數+輸入)81.07273.0165.3978.97225.2367.60HYM(參數)71.40225.7063.0885.26237.9762.91HYM(參數+輸入)68.57212.4963.2987.68254.0664.01BMA(參數+輸入+結構)90.19315.6056.7095.62271.1555.03驗證期XAJ(參數)62.32183.6471.2180.47188.6864.41XAJ(參數+輸入)73.81220.5074.6381.48190.3163.07HYB(參數)71.99289.9582.8480.66244.4074.13HYB(參數+輸入)82.66285.4471.7180.38249.2477.14HYM(參數)68.61270.2677.8886.77261.2376.62HYM(參數+輸入)69.16252.2376.8488.96278.2477.31BMA(參數+輸入+結構)90.97348.5669.7495.17303.0466.06

注：參數表示僅考慮參數不確定性；參數+輸入表示同時考慮參數和輸入不確定性；參數+輸入+結構表示綜合考慮模型參數、輸入和結構不確定性；黑體字表示最優值。

4 結 論

本文選用3個基于子流域的半分布式水文模型(新安江模型、混合產流模型和HYMOD模型)，采用貝葉斯模型平均法(BMA)開展水文模型不確定性及集合模擬實例研究。模擬結果表明： SCE-UA算法和SCEM-UA算法優選的參數使得3個模型都能取得較好的模擬結果，但SCEM-UA算法能給出模型參數的后驗概率分布，進行參數的不確定性分析，得到模型的概率預報區間，具有更大的優勢；模型在率定期的模擬結果要優于驗證期的模擬結果，說明模型在驗證期具有更大的不確定性；考慮模型輸入的不確定性，模擬徑流的精度沒有較大提高，但計算的概率預報區間精度有一定的提高；采用BMA方法，綜合考慮模型參數、模型輸入和模型結構的不確定性，其模擬結果優于單個模型，尤其在驗證期BMA集合模擬結果精度提高更大；此外，BMA集合模擬的概率預報性區間也優于單個模型，說明了BMA集合模擬對提高水文預報精度和進行不確定性分析的有效性；比較SCE-UA算法和SCEM-UA算法的BMA集合模擬結果，可以發現，基于SCEM-UA算法和BMA集合得到的概率預報區間更為準。綜合而言，采用SCEM-UA算法考慮模型參數的不確定性，再利用BMA集合各模型結果，能夠得到最為精確的預報值和最為可靠的概率預報區間。

□

[1] 徐宗學, 程 磊. 分布式水文模型研究與應用進展[J]. 水利學報, 2010,41(9):1 009-1 017.

[2] 徐宗學. 水文模型[M]. 北京:科學出版社, 2009.

[3] Beven K. Rainfall-Runoff modeling[M]. Chichester: Wiley, 2000.

[4] 嚴登華, 袁 喆, 王 浩, 等. 水文學確定性和不確定性方法及其集合研究進展[J]. 水利學報, 2013,44(1):73-82.

[5] 梁忠民, 戴 榮, 李彬權. 基于貝葉斯理論的水文不確定性分析研究進展[J]. 水科學進展, 2010,21(2):274-281.

[6] 王 浩, 李 揚, 任立良, 等. 水文模型不確定性及集合模擬總體框架[J]. 水利水電技術, 2015,46(6):21-26.

[7] Kavetski D, Kuczera G, Franks S W. Bayesian analysis of input uncertainty in hydrological modeling: 1. Theory. Water Resources Research, 2006,42(3).

[8] Ajami N K, Duan Q Y, Sorooshian S. An integrated hydrologic Bayesian multimodel combination framework: Confronting input, parameter, and model structural uncertainty in hydrologic prediction[J]. Water Resources Research, 2007,43(1).

[9] 王書功. 水文模型參數估計方法及參數估計不確定性研究[M]. 鄭州: 黃河水利出版社, 2010.

[10] 董磊華, 熊立華, 萬 民. 基于貝葉斯模型加權平均法的水文模型不確定性分析[J]. 水利學報, 2011,42(9):1 065-1 073.

[11] Liang Z M, Wang D, Guo Y, et al. Application of Bayesian Model Averaging Approach to Multi-model ensemble Hydrologic Forecasting[J]. Journal of Hydrologic Engineering, 2013,18(11):1 426-1 436.

[12] Deng L H, Xiong L H, Zheng Y F. Uncertainty analysis of coupling multiple hydrologic models and multiple objective functions in Han River, China[J]. Water Science and Technology, 2013,68(3):506-513.

[13] 杜新忠, 李敘勇, 王慧亮, 等. 基于貝葉斯模型平均的徑流模擬及不確定性分析[J]. 水文, 2014,34(3):6-10.

[14] 任立良, 江善虎, 袁 飛, 等. 水文學方法的研究與詮釋[J]. 水科學進展, 2011,22(4):586-592.

[15] 陸桂華, 吳 娟, 吳志勇. 水文集合預報試驗及其研究進展[J]. 水科學進展, 2012,23(5):728-734.

[16] 袁 喆, 嚴登華, 楊志勇, 等. 集合建模在徑流模擬和預測中的應用[J]. 水利學報, 2014,45(3):351-359.

[17] 趙人俊. 流域水文模擬-新安江模型與陜北模型[M]. 北京: 水利電力出版社, 1984.

[18] Ren L L, Zhang W, Li C H, et al. Comparision of Runoff Parameterization Schemes with Spatial Heterogenetity across Different Temporal Scales in Semihumid and Semiarid Regions[J]. Journal of Hydrologic Engineering, 2008,13(5):400-409.

[19] 宋文獻, 江善虎, 楊春生, 等. SCE-UA、SCEM-UA算法對改進的HYMOD模型模擬結果的影響[J]. 水電能源科學, 2013,31(11):17-20.

[20] Duan Q, S Sorooshian, V Gupta. Effective and efficient global optimization for conceptual rainfall-runoff models[J]. Water Resources Research, 1992,28(4):1 015-1 031.

[21] Vrugt J A, H V Gupta, W Bouten, 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.

[22] Roberts G O, Tweedie R L. Geometric convergence and centra limit theorems for multidim ensional Hastings and Metropoli algorithms[J]. Biometrika, 1996,83(1):95-110.

[23] Xiong L H, Wan M, Wei X J, et al. Indices for assessing the prediction bounds of hydrological models and application by generalized likelihood uncertainty estimation[J]. Hydrological Science Journal, 2009,54(5):852-871.