基于SCEM-UA算法和全局敏感性分析的水文模型參數優選不確定性研究*

曹飛鳳，張世強，許月萍，樓章華

(1.浙江省水利水電工程局，浙江 杭州 310020；2.浙江大學 建筑工程學院水文與水資源工程研究所，浙江 杭州 310058；3.中國科學院寒區旱區環境與工程研究所，甘肅 蘭州 730000)

水文模型參數的非線性和相關性導致參數分布存在多個局部最優解和不連續可導點，采用傳統的貝葉斯算法推求其顯性聯合概率分布十分困難，而Markov鏈蒙特卡羅法(Markov Chain Monte Carlo，MCMC)較傳統的統計方法(如矩估計法，極大似然估計法等)能更好的處理非線性復雜問題的求解，因此，近年來被廣泛應用于流域水文模型參數后驗分布的推算[1-6]。

Duan等[7]提出的SCE-UA算法致力于高效且具魯棒性的搜索參數全局最優解，許多研究都證明了SCE-UA算法的高效性[8-9]，但是該算法容易將搜索范圍聚焦到單個最佳參數周圍的小區域而陷入局部最優[10]。SCE-UA算法與其他傳統參數優選方法一樣更多關注參數優化的效率問題，而對模型參數的不確定性分析較為缺乏。由于水文數據和模型本身都有一定的置信區間，只關注參數優化效率而忽視對模型輸入、參數及結構等的不確定性分析，將導致模型輸出缺乏可靠性。基于上述考慮，Vrugt等[10]提出用MCMC理論的自適應Metropolis算法[11]取代SCE-UA算法中的坡降單純形法，形成了SCEM-UA (Shuffled Complex Evolution Metropolis-UA)參數全局優化算法。該算法采用自適應Metropolis采樣器，使新樣本的推薦分布不僅僅依靠先驗分布，而是在計算過程不斷更新，從而使該算法的計算效率大大提高，并能夠有效獲取參數的后驗概率分布。RSA(Regional Sensitivity Analysis)是全局敏感性分析的一種重要方法，由Hornberger和Spear于1981年提出的[12]，該方法側重分析參數空間分布的復雜性與相關性對模型輸出的響應，已被廣泛應用于水文模型參數的敏感性分析。如何將優化算法與參數的敏感性分析相結合，進一步提高模型參數優化的效率和可靠性，是目前水文模型研究方面的一個熱點。本文以岷江流域為案例，基于SCEM-UA算法和RSA方法研究了水文模型參數的后驗分布和不確定性，并在最后給出了置信區間為90%的模擬結果。

1 SCEM-UA算法

1.1 參數后驗概率密度估計

降雨徑流模型的不確定性主要包括輸入的不確定、參數的不確定和模型結構的不確定性[13]。考慮到模型誤差來源很難被定量化，一般常用如式(1)所示的簡單集總變量來表示，式中y代表模型輸出的N×1維向量；ε代表模型輸入的N×r維矩陣；θ代表n個模型參數向量。因此，對公式(1)問題的優化就轉化為求解E(θ)={e(θ)1,e(θ)2,…,e(θ)N}的最小解問題。假設殘差e(θ)相互獨立并服從高斯分布，可推求參數θ(t)的似然度，如式(2)所示，式中，λ表示殘差的分布類型，推導過程中取λ=0，即表示殘差服從正態分布；σ為殘差序列方差，∝表示正比于，假設p(θ)∝σ-1，可進一步求得如式(3)所示的θ(t)后驗概率密度[14]。

y=f(ε|θ)+e

(1)

(2)

(3)

1.2 算法原理

1.3 SCEM-UA算法的參數設置

(4)

2 案例研究

2.1 數據

本文選取岷江高場站以上流域作為研究例子。流域面積為135 400 km2，多年平均降雨量為1 265 mm，屬于濕潤地區。選取2000-2004年日平均流量用于模型率定，2005-2007年的資料用于模型驗證。為了減輕初始條件對模型計算的影響，2000-2004年的前20%時間序列數據不用于模型率定過程(Warm up 設置20%)。模型的輸入為根據流域內的氣象站和降水遙測站根據泰森多邊形計算的日面雨量數據，流域內馬爾康和小金站的日氣溫觀測數據，用于與模擬流量數據對比的是高場站2000-2007年的日流量數據。

2.2 模型

模型包括產流和匯流部分。產流部分用于將總降雨量轉化為有效降雨量，匯流部分將有效降雨量轉化為可與觀測流量對比的流域出口斷面流量。本研究采用Ye等人[17]提出的改進流域濕度損失方程計算產流,該模型是對Jakeman[18]提出的流域濕度指標模型的改進，匯流部分采用兩個平行線性水庫模型進行計算，改進后的模型具有更嚴格的統計意義，能被廣泛應用于濕潤地區的降雨徑流計算，詳見文獻[17]。需要率定的模型參數為9個，包括6個產流模型參數(τ、pref、mf、l、p、Sin)及3個匯流模型參數(kq、ks、q)。表1給出了模型中9個參數的先驗取值范圍。以上參數范圍的設置基于模型參數的物理意義并結合其在其它與研究區具有相似環境地區的應用經驗，以及基于作者對岷江流域地區水文模型的多次計算分析所得出的較為合適的范圍取值。模型的實現基于英國帝國理工學院土木與環境系自行研發的RRMT(Rainfall-Runoff Modelling Toolbox)[19]降雨徑流模型工具包。

2.3 目標函數及似然值定義

采用如式(5)所示的Nash-Sutcliffe法作為本研究的目標函數，用于評價實測流量與計算流量的擬合程度。Nash-Sutcliffe法可以集合整個率定過程的殘余誤差，峰值流量對結果的貢獻較大。NSE*越靠近零，擬合度越高。本研究全局敏感性分析基于式(6)所示的似然函數，該似然函數的定義沒有嚴格意義上的統計依據，可被視為一種相對概率，用于表達參數組能描述現實水文系統的概率情況。式(6)中i表示第i次取樣，N表示NSE*小于給定臨界值的對應參數組個數。

(5)

(6)

表1 降雨徑流概念模型參數取值及其物理意義

3 結果與討論

3.1 基于SCEM-UA算法的參數優選分析

圖1 比例縮小系數演化過程

圖2 參數τ和mf以及q對應的Markov鏈演化過程

圖3 參數后驗分布直方圖

選取由SCEM-UA算法取樣獲得的對應NSE*小于臨界值1的模型參數組，可得如圖3所示的邊緣后驗概率分布。圖3中每個參數的取值范圍被等分成20個直方圖，每個直方圖所對應的概率根據式(2)或式(3)求得。從圖中可知，參數l、p、kq、ks及q的后驗分布成偏正態分布，其余4個參數的后驗分布基本為不規則分布，但都具有峰值。這說明SCEM-UA算法能確保Markov鏈收斂于高概率后驗分布區，從而使參數后驗分布存在顯著的峰值。

SCEM-UA算法取樣獲取30 000個樣本，選取比例縮小系數小于1.2的樣本，進行參數后驗分布統計，結果如表2所示。表2 共統計了參數后驗分布均值、標準差、最大值、最小值。SCEM-UA算法獲得參數的后驗分布而不是一組唯一的最優參數解，一般認為SCEM-UA算法的對應最優解應該位于最高后驗概率分布區內，因此本文對比分析了SCEM-UA算法獲得的參數最大后驗概率對應參數值和SCE-UA算法獲取的最佳參數值。表中的均值和SCEM-UA算法獲得的參數最大后驗概率對應參數值較為接近，說明參數分布是朝著高概率區進行收斂。表2中最大值和最小值基本和參數取值范圍相一致，也說明了Markov鏈具有各態遍歷性。SCEM-UA列和SCE-UA列對應的參數值有些很接近，相對誤差較大的原因是由參數本身的辨識性不強決定的。

3.2 敏感性分析結果比較

表2 參數后驗分布統計情況1)

1)SCEM-UA列代表最大后驗概率對應的參數取值，SCE-UA列代表由SCE-UA算法獲取的最佳參數取值。

圖4 全局敏感性分析圖

圖5 模型率定過程下90%置信區間的計算值與觀測值的擬合情況 (20000101-20041231，Warm up=20%)

圖6 模型驗證過程下90%置信區間的計算值與觀測值的擬合情況 (20050101-20071231)

選取由SCEM-UA算法取樣獲得的NSE*小于臨界值1的模型參數組，基于式(6)所定義的似然值可得到如圖4所示的敏感性分析圖。圖4中參數組根據對應目標函數值從小到大等分成10組，每組根據式(6)計算對應的似然值，將似然值歸一化后獲得累積頻率分布。圖4紫色曲線對應似然度最高的邊緣累積分布曲線，藍色曲線代表似然度最低的邊緣累積分布曲線，兩曲線的水平距離相隔越遠代表此參數敏感性越強，反之，敏感性越弱。由圖5可知，流域蓄水量與有效降雨量之間的非線性相關系數p的敏感性最強，有效降雨量的臨界值系數l的敏感性次之，其次為參數ks及q。對比分析圖2-4可知，參數敏感性和后驗分布情況密切相關，敏感性強的參數其Markov鏈在演化過程中收斂于范圍較窄的高概率后驗區，其邊緣后驗概率密度分布存在明顯的峰值，相反，敏感性弱的參數對應的Markov鏈分布在較寬的區間內，其后驗概率分布較為平坦且無規律可循，從而導致模型參數的不確定性大大增強。

3.3 模型輸出不確定性分析

由SCEM-UA算法獲取30 000個參數組，設置目標函數值臨界值為1(NSE*)，小于臨界值的對應參數組為8 764個，對8 764個參數組分別進行計算，每個步長(1 440 min)對應8 764個計算流量值，根據這些流量數據獲取每個步長的流量分布情況，并設定置信區間為90%(置信度上限為95%，下限為5%)的模型計算不確定性區間。由圖5和圖6可知，置信區間為90%的計算結果基本較好地包含了觀測值，根據誤差情況來看，模擬精度較為良好，但整個率定、驗證過程觀測流量的峰值有躍出90%置信區間的情況，率定期第三個峰處、驗證期第一及第三個峰處的擬合情況有待進一步提高，分析原因，模型結構本身和觀測流量值(系統性偏差)存在的誤差以及單目標函數的選取都會影響實測流量與計算流量的擬合結果。此外，流域上的許多人類活動影響也是造成上述現象的重要原因之一，例如流域中的水利工程設施(中小型水庫、水土保持措施工程、跨流域調水工程等)，這些水利工程的規模影響流域產流參數和產流的整個過程，最終導致水文資料的不一致性。因此，系統科學的降雨-徑流模型研究，還應該對這些水利工程的控制流域面積、蓄水能力、流域分布位置、建造時期、管理調度方式等情況進行詳細研究。

4 結 論

1)本文以岷江流域作為研究實例，應用SCEM-UA算法和RSA方法分析了集總式降雨徑流模型的參數后驗分布及參數敏感性情況，對由于模型參數誤差引起的模型輸出不確定性進行了估計。結果發現，參數的強敏感性往往對應著參數邊緣后驗概率密度分布存在明顯的峰值，相反，參數的弱敏感性使得后驗概率分布較為平坦，也即參數的辨識性不強，從而導致參數選取的不確定性大大增強；

2)SCEM-UA算法的取樣無需給定先驗分布，更適用于參數先驗信息較少的復雜水文模型參數后驗分布的推算；

3)由參數后驗分布統計可知，Markov鏈朝著高概率后驗分布區進行收斂，且具有各態遍歷性；

4)相比SCE-UA算法而言，SCEM-UA算法更側重模型的不確定性分析，不僅可以搜索模型參數后驗分布情況，而且可推求一定置信區間下的模型計算值與觀測值的擬合情況，為模型參數及輸出的不確定性分析奠定基礎，并可指導模型結構的進一步優化。

參考文獻：

[1]VRUGT J A, BOUTEN W, GUPTA H V, et al. Toward improved identifiability of hydrologic model parameters: The information content of experimental data [J]. Water Resources Research, 2002,38(12), 1312, doi: 10.1029/ 2001WR001118.

[2]MARSHALL L, NOTT D, SHARMA A. A comparative study of Markov chain Monte Carlo methods for conceptual rainfall-runoff modeling [J]. Water Resources Research, 2004, 40, w02501, doi:10.1029/2003WR002378.

[3]王建平， 程聲通， 賈海峰. 基于MCMC法的水質模型參數不確定性研究[J].環境科學, 2006, 27(1):24-30.

[4]VRUGT J A, BRAAK J F, CLARK M P, et al. Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlo simulation [J]. Water Resources Research, 2008, 44, W00B09, doi: 10.1029/2007WR006720.

[5]BLASONE R S, MADSEN H, ROSBJERG D. Uncertainty assessment of integrated distributed hydrological models using GLUE with Markov chain Monte Carlo sampling [J]. Journal of Hydrology, 2008,353：18-32.

[6]BLASONE R S, VRUGT J A, MADSEN H, et al. Generalized likelihood uncertainty estimation (GLUE) using adaptive Markov Chain Monte Carlo sampling [J].Advances in Water Resources, 2008,31(4): 630-648.

[7]DUAN Q Y, GUPTA V K, SOROOSHIAN S. Shuffled complex evolution approach for effective and efficient global minimization [J]. Journal of Optimization Theory and Application, 1993, 76(3): 501-521.

[8]BOYLE D P, GUPTA H V, SOROOSHIAN S. Toward improved calibration of hydrological models: Combing the strengths of manual and automatic methods [J]. Water Resources Research, 2000, 36(12): 3663-3674.

[9]馬海波，董增川，張文明，等.SCE-UA算法在TOPMODEL參數優選中的應用[J]. 河海大學學報:自然科學版, 2006, 34(4): 361-365.

[10]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), 1201, doi:10.1029/2002WR001642.

[11]HAARIO H, SAKSMAN E, TAMMINEN J. An adaptive Metropolis algorithm [J]. Bernoulli, 2001, 7(2):223-242.

[12]HORNBERGER G, SPEAR R. An approach to the preliminary analysis of environmental systems [J]. Environmental Management, 1981,12(1):7-18.

[13]曹飛鳳，許月萍，玄英姬，等.流域決策的不確定性分析[J].浙江大學學報：工學版，2009，43(1)：188-192.

[14]THIEMANN M, TROSSET M, GUPTA H, et al. Bayesian recursive parameter estimation for hydrological models [J]. Water Resources Research, 2001, 37 (10)：2521-2535.

[15]衛曉婧， 熊立華， 萬民，等. 融合Markov鏈-蒙特卡洛算法的改進通用似然不確定性估計估計方法在流域水文模型中的應用[J]. 水利學報, 2009, 40(4):464-480.

[16]GELMAN A, RUBIN D B. Inference from iterative simulation using multiple sequences [J].Statistic Science, 1992(7):457-472.

[17]YE W, BATES B C, VINEY N R, et al. Performance of conceptual rainfall-runoff models in low-yielding ephemeral catchments [J]. Water Resources Research, 1997, 33(1):153-166.

[18]JAKEMAN A J, HORNBERGER G M. How much complexity is warranted in a rainfall-runoff model? [J]. Water Resources Research, 1993, 29(8): 2637-2649.

[19]WAGENER T, LEES M J. Rainfall-Runoff Modelling Toolbox user manual [R]. UK: Department of Civil and Environmental Engineering, Imperial College London, 2001.

[20]WAGENER T, KOLLAT J. Numerical and visual evaluation of hydrological and environmental models using the Monte Carlo analysis toolbox [J]. Environmental Modeling and Software, 2007,22:1021-1033.