基于SODA方法的HyMOD模型不確定性分析

，, 2，

(1.中國長江三峽集團(tuán)公司 三峽樞紐建設(shè)運(yùn)行管理局，湖北 宜昌 443133；2.華中科技大學(xué) 水電與數(shù)字化工程學(xué)院，武漢 430074；3.長江水利委員會(huì) 水文局，武漢 430010)

基于SODA方法的HyMOD模型不確定性分析

李帥1，文小浩1, 2，杜濤3

(1.中國長江三峽集團(tuán)公司 三峽樞紐建設(shè)運(yùn)行管理局，湖北 宜昌 443133；2.華中科技大學(xué) 水電與數(shù)字化工程學(xué)院，武漢 430074；3.長江水利委員會(huì) 水文局，武漢 430010)

為了改進(jìn)水文建模過程中的不確定性處理，采用一種融合全局優(yōu)化和數(shù)據(jù)同化(Simultaneous Optimization and Data Assimilation, SODA)的混合框架，對HyMOD模型進(jìn)行了不確定性分析，并與經(jīng)典SCEM-UA方法進(jìn)行了比較。SODA方法具有如下特點(diǎn)：①具備較高的參數(shù)搜索效率和尋優(yōu)能力；②明確考慮包括輸入、輸出、參數(shù)以及模型結(jié)構(gòu)在內(nèi)的重要不確定性來源。SODA方法在渭河流域的實(shí)例應(yīng)用結(jié)果表明:與SCEM-UA方法相比，SODA方法不僅顯著提高了預(yù)報(bào)精度，而且推求出了性質(zhì)更為優(yōu)良的預(yù)報(bào)區(qū)間。SODA方法的成功應(yīng)用，有助于模型概念的改進(jìn)及對水文系統(tǒng)功能的理解。

SODA; HyMOD模型; 水文模型不確定性; SCEM-UA; 集合卡爾曼濾波(EnKF)

1 研究背景

水文模型的不確定性歷來都是水文科學(xué)研究中極具挑戰(zhàn)性的課題，尤其是在全球水文及氣候變化的研究中，水文模型不確定性研究更是成為水文科學(xué)發(fā)展中的前沿和熱點(diǎn)[1]。一般而言，水文模型中的參數(shù)無法直接觀測，往往需要通過歷史資料率定得到。早期提出的模型率定智能算法[2-4]，精力主要集中在尋求如何利用計(jì)算機(jī)的計(jì)算速度和尋優(yōu)能力來解決人工調(diào)參耗時(shí)且過于主觀的缺點(diǎn)。然而，這些方法往往試圖依靠一組模型參數(shù)最優(yōu)值來獲取模型輸出，這顯然是不切實(shí)際的。越來越多的研究發(fā)現(xiàn)，不同的模型參數(shù)組合可以使建模者獲得相同或者幾乎接近的模型輸出，即所謂的“異參同效”現(xiàn)象[5]。鑒于此，基于經(jīng)典Bayesian理論[6-7]、pseudo-Bayesian理論[8]、集理論[9]以及多重準(zhǔn)則理論[10]的方法被廣泛地應(yīng)用于模型參數(shù)及預(yù)報(bào)精度的不確定性推求。特別是基于經(jīng)典Bayesian理論的SCEM-UA (Shuffled Complex Evolution Metropolis-University of Arizona) 方法[6]，通過參數(shù)空間的有效探索，使Markov鏈能夠朝著高概率密度區(qū)進(jìn)化，從而推導(dǎo)出具有顯著統(tǒng)計(jì)特征的水文模型參數(shù)的后驗(yàn)分布[11]。然而，這種方法的主要缺點(diǎn)就在于在人為忽略模型輸入不確定性的同時(shí)，將模型輸出和結(jié)構(gòu)誤差簡單地集總成一個(gè)單一的零均值、同方差的白噪音項(xiàng)。

20世紀(jì)60年代發(fā)展起來、90年代開始應(yīng)用于陸面過程的順序數(shù)據(jù)同化(Sequential Data Assimilation, SDA)方法，明確考慮包括輸入、輸出、參數(shù)以及模型結(jié)構(gòu)在內(nèi)的重要不確定性來源，具有克服普遍存在于傳統(tǒng)模型率定策略中這種缺陷的潛力[12]。與傳統(tǒng)的模型率定策略相比，SDA方法在新的觀測數(shù)據(jù)獲得的前提下，通過連續(xù)更新模型狀態(tài)來實(shí)時(shí)評價(jià)并改進(jìn)模型預(yù)報(bào)精度。而在眾多SDA方法中，基于全局卡爾曼理論的EnKF(Ensemble Kalman Filter)方法因其概念簡單、計(jì)算效率高以及應(yīng)用相對容易等特點(diǎn)，被廣泛地應(yīng)用于土壤水分反演[13]、模型參數(shù)估計(jì)[14-15]、洪水預(yù)報(bào)[16]以及雪水當(dāng)量估計(jì)[17]等方面。

因此，Vrugt等[18]提出將2種方法融合起來，構(gòu)造一個(gè)可以同時(shí)進(jìn)行參數(shù)優(yōu)化和數(shù)據(jù)同化(Simultaneous Optimization and Data Assimilation, SODA)的集成框架，即外部循環(huán)采用SCEM-UA方法進(jìn)行參數(shù)尋優(yōu)，同時(shí)內(nèi)部循環(huán)采用EnKF方法進(jìn)行狀態(tài)變量遞歸估計(jì)。本文在系統(tǒng)研究SODA方法的基礎(chǔ)上，結(jié)合董磊華等[6]、衛(wèi)曉婧等[11]、Xiong等[19]提出的覆蓋率、平均帶寬以及平均偏移幅度等指標(biāo),對SODA方法與傳統(tǒng)SCEM-UA方法進(jìn)行綜合評價(jià)和全面比較，試圖從算法結(jié)構(gòu)和模擬結(jié)果比較上解釋SODA方法的獨(dú)特之處，以期為后續(xù)水文模型不確定性分析方法的改進(jìn)提供一定的思路。

2 研究方法

2.1 SCEM-UA方法

SCEM-UA方法是Vrugt等人于2003年提出的一種融合SCE-UA (Shuffled Complex Evolution-UA)[2]算法和MCMC(Markov Chain Monte Carlo)方法的全局優(yōu)化算法。針對SCE-UA算法只關(guān)注參數(shù)優(yōu)化效率而忽視對模型輸入、結(jié)構(gòu)和參數(shù)不確定性分析的不足，SCEM-UA算法采用基于MCMC理論的Metropolis-Hasting搜索策略[20-21]取代SCE-UA算法中的坡降單純形法，進(jìn)而在進(jìn)化產(chǎn)生子代的過程中，通過構(gòu)造Markov鏈，使參數(shù)推薦分布有效地朝目標(biāo)后驗(yàn)概率分布演化，從而大大提高了參數(shù)后驗(yàn)概率分布的搜索效率[22]。SCEM-UA算法的基本原理和步驟為[23]：首先，在選定模型參數(shù)的可行空間內(nèi)，隨機(jī)產(chǎn)生參數(shù)集種群;其次，利用貝葉斯統(tǒng)計(jì)推斷對每個(gè)參數(shù)集計(jì)算其后驗(yàn)概率密度，并按降序排列；隨后，將整個(gè)種群劃分為q個(gè)復(fù)合形Ck(k=1,2,…,q)，每個(gè)Ck對應(yīng)一條Markov鏈Sk(k=1,2,…,q)，鏈的起始點(diǎn)為后驗(yàn)概率密度第k高的個(gè)體；然后，利用多元正態(tài)分布形成新的候選點(diǎn)，并用Metropolis-annealing準(zhǔn)則檢驗(yàn)候選點(diǎn)是否可加入當(dāng)前復(fù)合形，如果可行，則用該點(diǎn)隨機(jī)取代復(fù)合體中已經(jīng)存在的點(diǎn);最后，將更新后的復(fù)合形摻混，按照后驗(yàn)概率密度由高到低的順序重新進(jìn)行排列，如此反復(fù)迭代，直到滿足收斂條件，從而獲取有效的參數(shù)后驗(yàn)分布信息。

2.2 集合卡爾曼濾波(EnKF)方法

集合卡爾曼濾波(Ensemble Kalman Filter，EnKF)是一種基于MCMC的順序數(shù)據(jù)同化方法。該方法可以根據(jù)順序時(shí)刻的觀測資料，對狀態(tài)向量中的模型變量和模型參數(shù)進(jìn)行實(shí)時(shí)修正，并實(shí)時(shí)跟蹤流域水文行為的動(dòng)態(tài)變化[24]。由于順序數(shù)據(jù)同化框架綜合考慮了水文模型中包括參數(shù)、輸入、輸出以及模型結(jié)構(gòu)在內(nèi)的重要不確定性來源，因此，EnKF在減少水文預(yù)報(bào)不確定性方面表現(xiàn)出了一定的優(yōu)勢[25]。

圖1 用于狀態(tài)變量遞歸估計(jì)的集合卡爾曼濾波(EnKF)方法流程Fig.1 Flowchart of the ensemble Kalman filter(EnKF) used to recursively estimate state variables

EnKF方法流程如圖1所示。其具體步驟如下：

(1) 生成狀態(tài)初始集合。從試驗(yàn)?zāi)Ｐ蚼個(gè)狀態(tài)變量的先驗(yàn)分布中，隨機(jī)抽取N組樣本ψj，j=1,…,N，并將其存儲(chǔ)在矩陣A=[1∶m,1∶N]中，即

A=(ψ1,…,ψN)?Rm×N。

(1)

(2) 傳遞狀態(tài)。將A中的N組集合元素通過非線性模型結(jié)構(gòu)f(·)隨時(shí)間t向前傳遞，即

(2)

(3)

(4) 計(jì)算平均集合預(yù)報(bào)誤差。當(dāng)t時(shí)刻的觀測值已知時(shí)，計(jì)算平均集合預(yù)報(bào)誤差zt為

(4)

(5)

D=(y1,…,yN)?R1×N。

(6)

(6) 狀態(tài)更新。在每一觀測步，用如下方程更新A中的每個(gè)集合元素，即

A=A+Kt[D-H(A)] 。

(7)

式中Kt為卡爾曼增益矩陣，并由式(8)近似逼近，即

(8)

(7) 檢查停止準(zhǔn)則。如果t小于同化總時(shí)長T，則返回步驟(2)，否則終止。

2.3 SODA方法

基于以上2種方法，一種同時(shí)進(jìn)行全局優(yōu)化和數(shù)據(jù)同化的雜交方法應(yīng)運(yùn)而生。圖2給出了SODA方法的流程。

圖2 SODA方法流程Fig.2 Flowchart of the simultaneous optimization and data assimilation (SODA) method

從圖2中可以看出，為了成功執(zhí)行SODA算法，使SCEM-UA算法和EnKF方法能夠有機(jī)地結(jié)合起來，需要用合適的判斷準(zhǔn)則來描述由EnKF方法同化得到的模型平均集合預(yù)報(bào)流量與觀測流量之間的接近程度。為了充分提取必要數(shù)據(jù)信息并將其有效地轉(zhuǎn)化到模型參數(shù)估計(jì)上，SODA方法不僅融合了SCEM-UA算法參數(shù)全局搜索以及EnKF方法計(jì)算效率上的優(yōu)勢，而且對水文模型輸入、輸出、參數(shù)以及模型結(jié)構(gòu)的不確定性整體上提供了一個(gè)較好的處理。本文采用水文預(yù)報(bào)中常用的Nash-Sutcliffe 模型效率系數(shù)NSE作為后驗(yàn)概率密度函數(shù)[26]，即

(9)

3 應(yīng)用實(shí)例

3.1 試驗(yàn)流域

試驗(yàn)流域選用黃河最大支流渭河流域。渭河位于33°50′N—37°18′N，104°00′E—110°20′E之間，流域總面積13.48萬 km2，干流全長818 km[27]。渭河流域地處干旱地區(qū)和濕潤地區(qū)的過渡地帶，多年平均氣溫10～13 ℃之間，屬大陸性季風(fēng)氣候。多年平均降水量572 mm，多年平均水面蒸發(fā)量660～1 600 mm，陸面蒸發(fā)量500 mm左右。本文選取該流域1999年1月1日—2004年12月31日共6 a的日降雨徑流資料用于模型率定，2005年1月1日—2009年12月31日的資料用于模型驗(yàn)證。與此同時(shí)，為了減弱初始條件對模型的影響，率定期前3個(gè)月的降雨徑流資料用來預(yù)熱模型。

3.2 HyMOD模型

HyMOD模型是一個(gè)基于超滲蓄滿產(chǎn)流機(jī)制的集總式概念性水文模型[28]。該模型因其結(jié)構(gòu)簡單，參數(shù)個(gè)數(shù)較少，便于分析模型不確定性，故在水文模型不確定性評價(jià)中得到了較為廣泛的應(yīng)用[14-15, 29-30]。其基本思想是基于流域蓄水能力曲線,即流域被認(rèn)為是無數(shù)個(gè)點(diǎn)的集合，且每點(diǎn)都有確定的土壤蓄水容量。由于土壤類型及深度在整個(gè)流域上存在一定的空間變異，因此，土壤蓄水容量在整個(gè)流域上也是變化的。圖3給出了HyMOD模型的結(jié)構(gòu)。從圖3中可以清晰地歸納出模型的產(chǎn)匯流機(jī)制：流域在經(jīng)歷過某場降雨過程(P(t))后，超過流域最大土壤蓄水能力Cmax的部分(ER1(t))將無法滲流，其作為超滲部分直接排走，進(jìn)入快速徑流響應(yīng)水箱；而滲入土壤中的凈雨(扣除雨期蒸發(fā)ET(t))，在滿足不同位置土壤蓄水容量c(t)以后，超出的部分(ER2(t))按照快慢速徑流劃分比例系數(shù)α，分別分配給快速響應(yīng)和慢速響應(yīng)水箱。通過快速流和慢速流2組平行水箱調(diào)蓄過后得到的徑流之和形成流域出口總流量Q(t)。模型各參數(shù)物理意義及取值范圍如表1所示。

圖3 HyMOD模型結(jié)構(gòu)Fig.3 Structure of HyMOD model

參數(shù)物理意義取值范圍Cmax/mm流域土壤最大蓄水能力[0．00，500．00]bexp土壤蓄水容量空間變化率[0．10，2．00]α快慢速徑流劃分比例因子[0．10，0．99]Rs/d慢速響應(yīng)水箱殘留時(shí)間[0．01，0．10]Rq/d快速響應(yīng)水箱殘留時(shí)間[0．10，0．99]

HyMOD模型共有5個(gè)狀態(tài)變量，分別為流域土壤含水量S(t)、3個(gè)快速徑流響應(yīng)水箱蓄水量S1(t)、S2(t)和S3(t)及1個(gè)慢速徑流響應(yīng)水箱蓄水量Ss(t)。其中，S(t)存在一定的物理約束，即不能超過最大土壤含水量Smax(如圖3中斜陰影部分面積所示)。

3.3 模型結(jié)構(gòu)、輸入和輸出誤差項(xiàng)的確定方法

SODA方法與SCEM-UA方法的核心區(qū)別就在于SODA方法在融合SCEM-UA方法參數(shù)全局搜索的同時(shí)，還引入了EnKF數(shù)據(jù)同化操作，而EnKF方法成功實(shí)施的關(guān)鍵取決于模型結(jié)構(gòu)、輸入以及輸出誤差的合理選取[18]。基于Evensen經(jīng)典論述[31]推薦，本文中模型結(jié)構(gòu)誤差項(xiàng)(即式(2)中的ωt-1)按照式(10)確定,即

(10)

其中：

(11)

(12)

式中：st為模型的隨機(jī)驅(qū)動(dòng)項(xiàng)；ηt-1為服從均值為0、方差為1的正態(tài)分布抽樣向量，即ηt-1～N(0,1)；系數(shù)φ∈[0,1)，表示一階自相關(guān)；Δt表示時(shí)間步長；l表示在每個(gè)時(shí)間單元Δt中的步數(shù)，lΔt=1；σ為模型誤差的標(biāo)準(zhǔn)差，本文中σ=0.5。

本文輸入變量(面平均降水和潛在蒸散發(fā))的誤差項(xiàng)由式(13)和式(14)確定[16]，即:

Ptrue=Pinput+δP，

δP～N(0,(0.15Pinput+0.2)2) ;

(13)

Etrue=EinputδE，δE～U(0,2) 。

(14)

式中:Pinput和Einput分別為觀測面平均降水量(mm)和潛在蒸散發(fā)量(mm)；Ptrue和Etrue為考慮了輸入誤差之后的面平均降水量(mm)和潛在蒸散發(fā)量(mm)。

本文中流量觀測誤差的估計(jì)遵照Georgakakos等[32]的假設(shè)，認(rèn)為流量觀測誤差的標(biāo)準(zhǔn)差為觀測值Qobs的0.1倍，即

Qtrue=Qobs+δQ，δQ～N(0,(0.1Qobs)2) 。

(15)

式中:Qobs為觀測流量(m3/s)；Qtrue為考慮了觀測誤差之后的流量(m3/s)。

4 精度指標(biāo)

4.1 預(yù)報(bào)精度評定指標(biāo)

為了全面評價(jià)SODA策略的整體表現(xiàn)并與傳統(tǒng)SCEM-UA方法進(jìn)行比較，除水文預(yù)報(bào)中常被用來評定預(yù)報(bào)精度的Nash-Sutcliffe效率系數(shù)NSE(見式(9))外，均方根誤差(RMSE)、Pearson相關(guān)系數(shù)R和總徑流量相對誤差BIAS被用來分析模擬結(jié)果,即：

(16)

(17)

(18)

4.2 確定性區(qū)間優(yōu)良性評定指標(biāo)

本文中主要采用覆蓋率(CR)、平均帶寬(B)和平均偏移幅度(D)3個(gè)常用指標(biāo)來分析比較SODA和SCEM-UA算法的預(yù)報(bào)區(qū)間，即

(19)

(20)

(21)

其中，

(22)

式中：Qlow,t為預(yù)測區(qū)間下界(m3/s)；Qup,t為預(yù)測區(qū)間上界(m3/s)。由上述公式不難看出，覆蓋率CR越接近于100%，平均帶寬B越窄，且平均偏移度D越小，不確定性區(qū)間性質(zhì)越優(yōu)良。

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

HyMOD模型參數(shù)的不確定性分析結(jié)果歸納于圖4—圖6，以及表2中，并作如下分析討論。

5.1 似然函數(shù)及模型參數(shù)的統(tǒng)計(jì)學(xué)性質(zhì)

為了高效運(yùn)行SODA算法來有效估計(jì)HyMOD模型參數(shù)，本文中SODA算法的初始參數(shù)設(shè)置為：初始樣本取樣數(shù)為250組，復(fù)合形個(gè)數(shù)為5個(gè)，設(shè)置Markov鏈長為50。用于狀態(tài)遞歸估計(jì)的EnKF集合個(gè)數(shù)為50，共計(jì)3 000次迭代計(jì)算被用來估計(jì)模型參數(shù)的最優(yōu)值，使其收斂到一個(gè)較為合理的不確定性區(qū)間并推求其隱含的后驗(yàn)概率分布。

圖4 SODA方法得到的HyMOD模型后驗(yàn)概率密度函數(shù)值及參數(shù)變化情況Fig.4 Variations of values of the posterior probability density function and the HyMOD model parameters obtained by the SODA method

圖4給出了由SODA方法推導(dǎo)的模型后驗(yàn)概率密度函數(shù)及各參數(shù)隨迭代次數(shù)增加的變化情況。從圖4(a)中可以發(fā)現(xiàn)，經(jīng)過近250次迭代以后，后驗(yàn)概率密度函數(shù)NSE基本穩(wěn)定在80%以上，達(dá)到水文預(yù)報(bào)規(guī)范的乙等標(biāo)準(zhǔn)[33]。對應(yīng)于似然函數(shù)值的變化，模型各參數(shù)隨著迭代次數(shù)增加同樣表現(xiàn)出不同的波動(dòng)特征。從圖4中的(b)—(f)可看出，250次迭代以后，后驗(yàn)概率密度函數(shù)維持在較高水平并向更高水平進(jìn)一步改善，而各參數(shù)也在一定的不確定性區(qū)間波動(dòng)并不斷窄化這一區(qū)間，這從側(cè)面體現(xiàn)了“異參同效”現(xiàn)象。在迭代繼續(xù)進(jìn)行的過程中，快速流水箱參數(shù)Rq相比于其他參數(shù)，以最快速度向最小不確定性區(qū)間收斂，因此表現(xiàn)出最高的可識(shí)別度，這與Moradkhani得到的結(jié)論一致[14]。與此相對，慢速流水箱參數(shù)Rs則表現(xiàn)出最低的可識(shí)別度。

由圖3可知，土壤蓄水空間變化參數(shù)bexp與土壤最大蓄水能力參數(shù)Cmax存在較強(qiáng)的相互作用，表現(xiàn)出一定的互補(bǔ)關(guān)系，這可以從圖4(b)和圖4(c)中參數(shù)不確定性區(qū)間一定程度的互補(bǔ)演變過程中得到驗(yàn)證。

5.2 參數(shù)后驗(yàn)分布

采用SCEM-UA和SODA算法推求的HyMOD模型參數(shù)的后驗(yàn)分布如圖5所示。從圖5中可以得出：第一，對于多數(shù)參數(shù)而言，SCEM-UA和SODA算法推求得到的參數(shù)空間的最高概率密度區(qū)間位置并不相同。假設(shè)模型結(jié)構(gòu)正確，且模型輸入和輸出數(shù)據(jù)在沒有任何誤差的情況下被觀測到，那么在模型率定的過程中就不需要進(jìn)行模型狀態(tài)校正，進(jìn)而由SCEM-UA和SODA算法識(shí)別出的模型參數(shù)空間的后驗(yàn)概率分布模式應(yīng)該完全一致。然而，在模型率定過程中倘若沒有適當(dāng)考慮模型結(jié)構(gòu)和輸入誤差，將會(huì)引起參數(shù)估計(jì)的墮落來補(bǔ)償這些誤差?；谝陨险摂?，本文認(rèn)為，由SODA算法識(shí)別出的模型參數(shù)更加不易被模型結(jié)構(gòu)和輸入誤差所污染，因而能夠更好地表征物理系統(tǒng)的潛在特性。

圖5 SCEM-UA和SODA方法推求的HyMOD模型參數(shù)后驗(yàn)分布Fig.5 Marginal posterior probability distributions of the HyMOD model parameters derived for the Weihe River Basin using the SCEM-UA algorithm and the SODA framework

圖6 檢驗(yàn)期2種方法推求的徑流過程及95%預(yù)報(bào)區(qū)間Fig.6 Streamflow processes with 95% forecast bounds derived by the SCEM-UA and SODA methods in the validation period

第二，當(dāng)允許在模型率定過程中加入模型狀態(tài)變量遞歸校正以后，HyMOD模型參數(shù)反而變得不易識(shí)別。由圖5可以看出，相比SCEM-UA算法，SODA算法推求得到的參數(shù)圍繞后驗(yàn)概率分布模式的離散程度明顯增加。從系統(tǒng)理論的角度來講，這種情況的出現(xiàn)，可以用SODA算法中部分徑流觀測“信息”被用來進(jìn)行狀態(tài)變量的遞歸估計(jì)來解釋。

5.3 結(jié)果比較

SODA算法的有效性還可以用許多其他不同的方式來證明，其中，一個(gè)最直接的方式就是與基于經(jīng)典Bayesian理論的SCEM-UA方法比較水文預(yù)報(bào)中洪水流量過程線的不確定性區(qū)間。檢驗(yàn)期典型峰值和非典型峰值徑流預(yù)測區(qū)間如圖6所示。

由圖6可看出，SODA方法推導(dǎo)的預(yù)測區(qū)間較窄，并能夠更精確地覆蓋實(shí)測數(shù)據(jù)，對于HyMOD模型不確定性的估計(jì)效果要明顯優(yōu)于傳統(tǒng)的SCEM-UA方法。基于SCEM-UA和SODA方法的HyMOD模型率定期及檢驗(yàn)期徑流預(yù)報(bào)表現(xiàn)如表2所示。從表2中可以看出，不論是率定期還是檢驗(yàn)期，SODA方法比SCEM-UA方法表現(xiàn)明顯更優(yōu)。

表2 基于SCEM-UA和SODA方法的HyMOD模型徑流預(yù)報(bào)結(jié)果統(tǒng)計(jì)Table 2 Summary statistics of streamflow forecast resultsusing SCEM-UA and SODA methods

6 結(jié) 論

基于MCMC理論的SCEM-UA算法原理簡單，易于操作。其通過融合Metropolis-Hastings算法、控制隨機(jī)搜索、生物競爭進(jìn)化以及復(fù)合摻混等方法的優(yōu)點(diǎn)，使參數(shù)推薦分布有效地朝目標(biāo)后驗(yàn)概率分布演化，且在高概率密度區(qū)識(shí)別出最有可能的參數(shù)值。與此同時(shí)，SCEM-UA算法交換并行進(jìn)化序列間信息的能力，增加了子鏈通過可行參數(shù)空間的穿越速度，大大提高了參數(shù)后驗(yàn)分布的搜索效率。然而SCEM-UA算法過分依賴模型結(jié)構(gòu)，且將模型輸入-輸出形式上的不確定性主要分配給參數(shù)估計(jì)不確定性與一個(gè)殘差“模型”誤差項(xiàng)的組合。而基于全局優(yōu)化和數(shù)據(jù)同化的SODA方法，在SCEM-UA算法的基礎(chǔ)上，針對水文模型強(qiáng)非線性、閾值型的結(jié)構(gòu)特點(diǎn)，引入順序數(shù)據(jù)同化方法中概念簡單、計(jì)算效率相對較高的集合卡爾曼濾波(EnKF)方法進(jìn)行狀態(tài)變量的遞歸估計(jì)，保證了在不失參數(shù)搜索效率和尋優(yōu)能力的前提下，又較好地考慮了模型參數(shù)、輸入、輸出以及模型結(jié)構(gòu)的不確定性，顯著改善了模型參數(shù)估計(jì)和預(yù)測區(qū)間。

試驗(yàn)結(jié)果表明：就HyMOD模型而言，SODA方法能夠推導(dǎo)出性質(zhì)更優(yōu)良的不確定性區(qū)間，從而更好地分析模型參數(shù)的不確定性，有助于改進(jìn)水文建模理念并更加深入地理解水文系統(tǒng)功能。

[1] 林凱榮, 陳曉宏. 基于FCM-SCEMUA的水文模型參數(shù)不確定性估計(jì)方法[J]. 水利學(xué)報(bào),2010,41(10):1186-1192.[2]DUAN Q, SOROOSHIAN S, GUPTA V. Effective and Efficient Global Optimization for Conceptual Rainfall-runoff Models[J]. Water Resources Research, 1992, 28(4): 1015-1031.

[3]KENNEDY J, EBERHART R. Particle Swarm Optimization[C]∥IEEE. Proceedings of IEEE International Conference on Neural Networks (IV). Perth, Australia. November 27-December 1, 1995: 1942-1948.

[4] HOLLAND J H. Adaptation in Natural and Artificial Systems[M]. Ann Arbor: University of Michigan Press, 1975.

[5]BEVEN K, BINLEY A. The Future of Distributed Models: Model Calibration and Uncertainty Prediction[J]. Hydrological Process, 1992, 6(3): 279-298.[6] 董磊華, 熊立華, 萬 民. 基于貝葉斯模型加權(quán)平均方法的水文模型不確定性分析[J]. 水利學(xué)報(bào), 2011, 42(9): 1065-1074.[7] VRUGT J A, GUPTA H V, BOUTEN W,etal. 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.

[8] FREER J, BEVEN K, AMBROISE B. Bayesian Estimation of Uncertainty in Runoff Prediction and the Value of Data: An Application of the GLUE Approach[J]. Water Resources Research, 1996, 32(7): 2161-2173.

[9] VRUGT J A, BOUTEN W, WEERTS A H. Information Content of Data for Identifying Soil Hydraulic Parameters from Outflow Experiment[J]. Soil Science Society of America Journal, 2001, 65(1): 19-27.[10]VRUGT J A, GUPTA H V, BASTIDAS L A,etal. Effective and Efficient Algorithm for Multi-objective Optimization of Hydrological Models[J]. Water Resources Research,2003,39(8),1214,doi:10.1029/2002WR001746.

[11]衛(wèi)曉婧, 熊立華, 萬 民, 等. 融合馬爾科夫鏈-蒙特卡洛算法的改進(jìn)通用似然不確定性估計(jì)方法在流域水文模型中的應(yīng)用[J]. 水利學(xué)報(bào), 2009, 40(4): 464-473.[12]王 文, 寇小華．水文數(shù)據(jù)同化方法及遙感數(shù)據(jù)在水文數(shù)據(jù)同化中的應(yīng)用進(jìn)展[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版), 2009, 37(5): 556-562.

[13]RIECHLE R H, MCLAUGHLIN D B, ENTEKHABI D. Hydrologic Data Assimilation with the Ensemble Kalman Filter[J]. Monthly Weather Review, 2002, 130(1):103-114.

[14]MORADKHANI H, SOROOSHIAN S, GUPTA H V,etal. Dual State-parameter Estimation of Hydrological Models Using Ensemble Kalman Filter[J]. Advances in Water Resources, 2005, 28(2): 135-147.[15]WANG D, CHEN Y, CAI X. State and Parameter Estimation of Hydrologic Models Using the Constrained Ensemble Kalman Filter[J]. Water Resources Research, 2009, 45(11), W11416, doi: 10.1029/2008WR007401.[16]WEERTS A H, EL SERAFY G Y H. Particle Filtering and Ensemble Kalman Filtering for State Updating with Hydrological Conceptual Rainfall-runoff Models[J]. Water Resources Research, 2006, 42(9), W09403, doi: 10.1029/2005WR004093.

[17]SU H, YANG Z, NIU G,etal. Enhancing the Estimation of Continental-scale Snow Water Equivalent by Assimilating MODIS Snow Cover with the Ensemble Kalman Filter[J]. Journal of Geophysical Research, 2008, 113, D08120, doi: 10.1029/2007JD009232.

[18]VRUGT J A, DIKS C G H, GUPTA H V,etal. Improved Treatment of Uncertainty in Hydrologic Modeling: Combining the Strengths of Global Optimization and Data Assimilation[J]. Water Resources Research, 2005, 41(1), W01017, doi: 10.1029/2004WR003059.

[19]XIONG L, WAN M, WEI X,etal. Indices for Assessing the Prediction Bounds of Hydrological Models and Application by Generalized Likelihood Uncertainty Estimation[J]. Hydrological Sciences Journal, 2009, 54(5): 852-871.

[20]HASTINGS W K. Monte-Carlo Sampling Methods Using Markov Chains and Their Applications[J]. Biometrika, 1970, 57(1): 97-109.

[21]METROPOLIS N, ROSENBLUTH A W, ROSENBLUTH M N,etal. Equations of State Calculations by Fast Computing Machines[J]. Journal of Chemical Physics, 1953,21(6):1087-1092.

[22]曹飛鳳, 嚴(yán)齊斌, 張世強(qiáng). 改進(jìn)SCEM-UA算法在概念性降雨-徑流模型參數(shù)優(yōu)選中的應(yīng)用[J]. 系統(tǒng)工程理論與實(shí)踐, 2012, 32(6): 1362-1368.

[23]陳南祥, 甘 甜, 杜青輝. SCEM-UA優(yōu)化算法在水文頻率分析中的應(yīng)用[J]. 西北農(nóng)林科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011, 39(8): 210-214.[24]南統(tǒng)超, 吳吉春．集合卡爾曼濾波估計(jì)水文地質(zhì)參數(shù)的局域化修正[J]. 水科學(xué)進(jìn)展, 2010, 21(5): 613-621.[25]LIU Y, GUPTA H V. Uncertainty in Hydrologic Modeling: Toward an Integrated Data Assimilation Framework[J]. Water Resources Research, 2007, 43, W07401, doi: 10.1029/2006WR005756.[26]NASH J E, SUTCLIFFE J V. River Flow Forecasting Through Conceptual Models[J]. Journal of Hydrology, 1970, 10(3): 282-290.[27]秦 濤, 汪 妮,解建倉, 等. 渭河流域?qū)氹u段水體納污能力及總量控制定額分析計(jì)算[J]. 長江科學(xué)院院報(bào), 2009, 26(8): 20-23.

[28]MOORE R J. The Probability-distributed Principle and Runoff Production at Point and Basin Scales[J]. Hydrological Sciences Journal, 1985, 30(2): 273-297.

[29]MORADKHANI H, HSU K L, GUPTA H,etal. Uncertainty Assessment of Hydrologic Model States and Parameters: Sequential Data Assimilation Using the Particle Filter[J]. Water Resources Research, 2005, 41(5), W05012, doi: 10.1029/2004WR003604.

[30]VRUGT J A, TER BRAAK C J F, CLARK M P,etal. Treatment of Input Uncertainty in Hydrologic Modeling: Doing Hydrology Backward with Markov Chain Monte Carlo Simulation[J]. Water Resources Research, 2008, 44(12), W00B09, doi: 10.1029/2007WR006720.

[31]EVENSEN G. Sequential Data Assimilation with a Nonlinear Quasi-geostrophic Model Using Monte Carlo Methods to Forecast Error Statistics[J]. Journal of Geophysical Research, 1994, 99(C5): 10143-10162.

[32]GEORGAKAKOS K P. A Generalized Stochastic Hydrometeorological Model for Flood and Flash Flood Forecasting: 2. Case Studies[J]. Water Resources Research, 1986, 22(13): 2096-2106.

[33]SL 250—2000,水文情報(bào)預(yù)報(bào)規(guī)范[S]. 北京: 中國水利水電出版社, 2000.

(編輯：陳 敏)

Uncertainty Assessment of HyMOD Model Using the Method ofSimultaneous Optimization and Data Assimilation

LI Shuai1, WEN Xiao-hao1, 2, DU Tao3
(1.Three Gorges Construction and Operation Management Department, China Three Gorges Corporation, Yichang 443133, China; 2.School of Hydropower and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, China; 3.Bureau of Hydrology, Changjiang Water Resources Commission, Wuhan 430010, China)

To improve the treatment of uncertainty in hydrological modeling, a hybrid framework of simultaneous optimization and data assimilation (SODA) was adopt to assess the uncertainty of HyMOD model in this paper, andthen was compared with the classical method of the Shuffled Complex Evolution Metropolis-UA (SCEM-UA). The strengths of the SODA can be described as follows: (1) high parameter search efficiency and explorative capabilities; (2) explicit treatment of the various important sources of uncertainty (i.e., input, output, parameter and model structure uncertainties) associated with the application of hydrological models. The results of the SODA applied in the Weihe River Basin demonstrate that in comparison to the performances of SCEM-UA, the SODA could notably improve the streamflow prediction efficiency, and also could derive more accurate prediction interval. The successful application of the SODA is helpful to improving model concepts and understanding of the functioning of hydrological systems.

simultaneous optimization and data assimilation (SODA);HyMOD;hydrological modeling uncertainty; Shuffled Complex Evolution Metropolis-University of Arizona(SCEM-UA); Ensemble Kalman Filter(EnKF)

2016-05-25；

：2016-09-06

國家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2016YFC0402306-01);國家自然科學(xué)基金項(xiàng)目(51539009)

李 帥(1987-)，男，湖北漢川人，工程師，博士,主要從事水文預(yù)報(bào)和水庫調(diào)度方面的研究與管理工作,(電話)0717-6763354(電子信箱)li_shuai@ctg.com.cn。

10.11988/ckyyb.20160519

2017,34(9):6-13,23

P333

：A

：1001-5485(2017)09-0006-08