鄭言峰

(青島市水利勘測設(shè)計研究院有限公司，山東 濟(jì)南 250013)

水文模型的不確定性評價已經(jīng)成為水文學(xué)者研究的熱點(diǎn)問題之一[1-4]。水文模型不確定性的來源是多方面的，總體上可以分為以下3類：實(shí)測資料(降雨、蒸發(fā)、徑流等)的不確定性、模型結(jié)構(gòu)的不確定性以及模型參數(shù)的不確定性。為了綜合分析以上因素對模型不確定性的影響，Beven和 Binley[5]于1992年提出了通用似然不確定性估計方法，即GLUE方法。該方法將Monte Carlo隨機(jī)取樣技術(shù)與Bayesian理論結(jié)合起來，原理簡單，易于操作。然而，GLUE方法在很多實(shí)際應(yīng)用中的結(jié)果[6]表明預(yù)測區(qū)間的覆蓋率遠(yuǎn)低于給定的置信水平。XIONG等[7]通過引入“模擬偏差曲線”對GLUE方法進(jìn)行了改進(jìn)，在一定程度上提高了區(qū)間覆蓋率，但同時也增加了區(qū)間帶寬。

數(shù)據(jù)同化方法將模型和觀測這兩種研究手段進(jìn)行了結(jié)合[8]，從而實(shí)現(xiàn)了在觀測資料更新的同時可以更新模型的計算結(jié)果。它的主要思想是將最新的觀測數(shù)據(jù)融合到所使用的模型之中，通過不斷更新的觀測數(shù)據(jù)來減小狀態(tài)變量模擬誤差，繼而使用更新之后的狀態(tài)變量重新預(yù)報從而來達(dá)到提高預(yù)報精度的目的。作為數(shù)據(jù)同化方法的典型代表，集合卡爾曼濾波(Ensemble Kalman Filter，EnKF)已被廣泛應(yīng)用到海洋、陸面和大氣數(shù)據(jù)同化領(lǐng)域之中[9]。EnKF方法融合了蒙特卡洛(Monte Carlo，MC)采樣技術(shù)，可以分別考慮模型觀測資料的不確定性以及模型結(jié)構(gòu)的不確定性。2005年，Moradkhani等[10]在EnKF的基礎(chǔ)上又加入了模型參數(shù)的更新(即dual EnKF)，用以描述模型參數(shù)的不確定性。

為了探索解決在提高區(qū)間覆蓋率的同時降低預(yù)測區(qū)間平均帶寬問題，本文將dual EnKF方法與GLUE方法2種水文模型不確定性研究方法進(jìn)行了耦合，得到了dual E-G方法，并以新安江模型為例，研究了dual E-G方法在漢江上游3個子流域(即甲河流域、天河流域和灞河流域)中的適用性。該成果可為水文模型預(yù)報、洪水決策風(fēng)險分析等方面的深入研究提供有益參考。

1 研究區(qū)概況

漢江又稱漢水，是長江中下游最大的一條支流，發(fā)源于秦嶺南麓，干流流經(jīng)陜西和湖北兩省，于武漢市匯入長江。漢江干流全長1 570多km，流域面積大約15.9萬km2。干流丹江口以上為上游，流域面積大約9.5萬km2。漢江流域河流密布，溝壑交織，分布著大小河流及山溝數(shù)以千計。漢江流域幅員廣闊，光、熱和水資源空間差異比較大，氣候特征多變，其降水主要來源于東南、西南2股暖濕氣流，因此降水變率比較大，是中國旱、澇災(zāi)害多發(fā)的地區(qū)之一。作為南水北調(diào)中線工程水源地的丹江口水庫，不僅承擔(dān)著湖北武漢、襄陽等23個縣市1億多人口及124萬多hm2耕地的防洪任務(wù)，還承擔(dān)著河南、河北、北京、天津等四省市的20多座城市的供水任務(wù)，意義重大。但是21世紀(jì)以來，漢江流域出現(xiàn)的干旱問題，嚴(yán)重影響了流域內(nèi)的生態(tài)環(huán)境與供水工程。研究漢江流域上降雨徑流模擬的不確定性，不僅可以為決策者提供充足的風(fēng)險信息，還可以實(shí)現(xiàn)水資源的合理配置，從而保證南水北調(diào)中線工程的順利實(shí)施。

因此，選擇了漢江上游的3個子流域(即甲河流域、天河流域和灞河流域)作為研究對象。3個子流域均為濕潤地區(qū)，且資料長度均為8年。前5年(1980—1984年，共1 827天)的流域?qū)崪y平均日降雨、日蒸發(fā)和日徑流作為率定期，后3年(1985—1987年，共1 095天)的流域?qū)崪y平均日降雨、日蒸發(fā)和日徑流作為檢驗(yàn)期。流域水文資料的具體情況見表1，3個子流域在漢江中的位置及流域見圖1。

表1 研究流域資料

圖1 甲河、天河及灞河在漢江中的位置及其流域

2 模型簡介及評價指標(biāo)

2.1 新安江模型簡介

對于研究水文模型的不確定性而言，概念性水文模型優(yōu)于系統(tǒng)模型和物理模型。概念性模型主要利用一些簡單的物理概念和經(jīng)驗(yàn)公式，來近似描述水流在流域的運(yùn)動狀態(tài)。系統(tǒng)模型主要是將所研究的流域或區(qū)間看作一種動力系統(tǒng)，利用輸入與輸出資料，建立某種數(shù)學(xué)關(guān)系，這種模型只關(guān)心模擬結(jié)果的精度，并不考慮輸入與輸出之間的物理關(guān)系。物理模型主要根據(jù)水流連續(xù)方程和動量方程來求解水流在流域的時間和空間變化規(guī)律，雖然模型物理意義更加明確，但其對資料要求過于嚴(yán)格，一般流域很難滿足。

新安江(三水源)模型是河海大學(xué)(原華東水利學(xué)院)趙人俊等在1973年提出來的降雨徑流模型，簡稱新安江模型。由于新安江模型采用蓄滿產(chǎn)流進(jìn)行產(chǎn)流計算，因此主要適用于濕潤、半濕潤地區(qū)。近幾十年，新安江模型不斷改進(jìn)，已成為中國應(yīng)用較為廣泛的一個流域水文模型。

新安江模型作為國內(nèi)常用的一種概念性水文模型，大部分模型參數(shù)都具有明顯的物理意義，不僅模擬精度高，而且對水文資料的要求比較低，對于日模型來說，只需要流域日蒸發(fā)資料、日降雨資料以及日徑流資料即可。新安江模型共有15個模型參數(shù)。①K：流域蒸散發(fā)能力Ep與蒸發(fā)皿實(shí)測值EM之比，簡稱蒸散發(fā)折算系數(shù)；②IMP：一般將整個流域面積分為透水面積與不透水面積，IMP是指全流域面積中不透水面積所占的比例；③B：蓄水容量曲線的指數(shù)；④WUM：流域上層蓄水容量(mm)；⑤WLM：流域下層蓄水容量(mm)；⑥WDM：流域深層蓄水容量(mm)；⑦C：深層蒸散發(fā)系數(shù)；⑧SM：自由水蓄水容量(mm)，反映表層土蓄水能力的大小，直接決定著地面徑流和地下徑流所占的比重；⑨EX：自由水蓄水容量曲線的指數(shù)；⑩KG：地下水出流系數(shù)；KI：壤中流出流系數(shù)，主要取決于表層土的滲透性；CG：地下徑流的消退系數(shù)，反映地下徑流的消退規(guī)律；CI：壤中流的消退系數(shù)，反映壤中流的消退規(guī)律；nn：反映流域調(diào)蓄能力的參數(shù)，相當(dāng)于線性水庫的個數(shù)或水庫的調(diào)節(jié)次數(shù)；KK：線性水庫的蓄泄系數(shù)，相當(dāng)于流域匯流時間的參數(shù)，具有時間因次。

2.2 似然函數(shù)的選擇

在對水文模型進(jìn)行率定時，需要定義一個評價模擬實(shí)測流量好壞的標(biāo)準(zhǔn)，即似然函數(shù)。似然函數(shù)必須滿足一定的特點(diǎn)：當(dāng)模擬系列與實(shí)測系列沒有相似性時它的值應(yīng)為0，并且相似性越高，它的值應(yīng)越大。不同的似然函數(shù)率定得到的參數(shù)也會有所不同，有的似然函數(shù)對高水的模擬效果較好，有的似然函數(shù)對低水的模擬效果較好，而納西效率系數(shù)因?qū)φw模擬效果較好而被廣泛應(yīng)用于水文模型參數(shù)估計研究中。因此，采用納西效率系數(shù)R2作為似然函數(shù)，表達(dá)見式(1)：

(1)

當(dāng)R2的值等于1時，模擬流量序列完全重現(xiàn)了相應(yīng)的實(shí)測流量序列，此時的模擬效果最好；當(dāng)R2的值等于0時，模擬流量序列對同期觀測的擬合程度與實(shí)測流量序列長時間尺度上的平均值對實(shí)測流量序列的擬合程度相當(dāng)；R2的值越小，模型模擬的序列對相應(yīng)觀測的擬合程度越差[11]。

2.3 預(yù)測區(qū)間評價指標(biāo)

為了反映水文模型模擬和預(yù)報實(shí)測流量的不確定性，每個實(shí)測流量的預(yù)測結(jié)果不應(yīng)該只是單純的點(diǎn)估計，而應(yīng)該是給定置信水平下的區(qū)間估計。XIONG等[12]從區(qū)間覆蓋率、區(qū)間帶寬和區(qū)間不對稱性3個方面詳細(xì)介紹了7個用于評價預(yù)測區(qū)間好壞的指標(biāo)。因此，從3個方面中各選擇一個比較常用的指標(biāo)(區(qū)間覆蓋率、平均帶寬和平均不對稱程度)來評價預(yù)測區(qū)間的性質(zhì)。它們的定義如下。

a)區(qū)間覆蓋率(CR)。區(qū)間覆蓋率是指實(shí)測流量落入預(yù)報區(qū)間的個數(shù)與實(shí)測流量數(shù)據(jù)總個數(shù)的比值。自從Beven提出GLUE方法以后，區(qū)間覆蓋率就被廣泛應(yīng)用于評價預(yù)測區(qū)間的好壞。CR越大，表示預(yù)測區(qū)間覆蓋率越高，CR最為理想的值為1，這意味著預(yù)測區(qū)間覆蓋了所有的實(shí)測流量數(shù)據(jù)。

b)平均帶寬(B)。平均帶寬是所有流量的預(yù)測區(qū)間上下界差值的平均值，單位：m3/s。

(2)

QU(j)、QL(j)分別表示第j時刻的流量預(yù)測區(qū)間上下界。對于指定的置信水平，在保證有較高的覆蓋率前提下，預(yù)報區(qū)間的平均帶寬越窄越好，此時的預(yù)測區(qū)間的不確定性也較小。

c)平均不對稱程度(S)。

(3)

理想的預(yù)測區(qū)間是不僅能覆蓋大部分的實(shí)測流量，而且預(yù)測區(qū)間的平均帶寬值較小，并且預(yù)測下界與實(shí)測流量之間的距離近似等于預(yù)測上界與實(shí)測流量之間的距離。

3 不確定性研究方法

3.1 GLUE方法

GLUE方法中的一個很重要的觀點(diǎn)是：導(dǎo)致模型模擬結(jié)果好壞的不是模型單個參數(shù)，而是模型的參數(shù)組合。目前GLUE方法已成為水文模型不確定研究的主流方法[1-3,13-14]。GLUE方法的主要思路如下：選定模型及似然函數(shù)，并設(shè)定似然閾值；在預(yù)先設(shè)定的參數(shù)分布(一般假定為均勻分布)取值空間內(nèi)，利用Monte Carlo隨機(jī)取樣獲得一定數(shù)量的有效參數(shù)組(即似然函數(shù)值高于似然閾值的參數(shù)組)；在給定的置信水平下，求出每個實(shí)測流量的預(yù)測區(qū)間。

GLUE方法共有3個參數(shù)，即似然函數(shù)閾值TH、有效參數(shù)組的個數(shù)M、置信水平α。每個參數(shù)的取值如下：TH=0.5，M=2000，α=0.8。

3.2 dual EnKF方法

在EnKF方法中假定水文模型的參數(shù)不隨時間變化，模型的參數(shù)可以通過模型率定獲得。但是由于模型驅(qū)動數(shù)據(jù)和觀測數(shù)據(jù)存在一定誤差，不能保證模型的模擬效果不會隨著時間的變化而有所差異。因此，有必要將模型參數(shù)的更新納入預(yù)報的過程之中[10,15-16]。模型參數(shù)的更新過程類似于狀態(tài)變量的更新過程。

dual EnKF方法的具體步驟如下。

步驟一選擇水文模型，準(zhǔn)備必要的驅(qū)動數(shù)據(jù)(流域平均日降雨資料P、流域平均日蒸發(fā)資料E)和流域?qū)崪y日流量資料Q，確定一組有效參數(shù)組θ0。

步驟二確定集合成員個數(shù)N、日蒸發(fā)的測量誤差超參數(shù)γE、日降雨的測量誤差超參數(shù)γP和日徑流的測量誤差超參數(shù)γQ。

步驟三分別對日蒸發(fā)、日降雨和日徑流資料加上一系列高斯白噪聲擾動。

(4)

(5)

(6)

步驟四計算每個成員的狀態(tài)變量先驗(yàn)值。

(7)

t=1時，

步驟五計算每個成員的模擬流量。

(8)

步驟六計算每個模型參數(shù)的增益。

(9)

步驟七計算每個成員t+1時刻的模型參數(shù)更新值。

(10)

步驟八使用更新之后的模型參數(shù)重新計算狀態(tài)變量值。

(11)

步驟九使用更新之后的模型參數(shù)重新計算模擬流量值。

(12)

步驟十計算每個狀態(tài)變量的增益。

(13)

步驟十一計算每個成員的狀態(tài)變量后驗(yàn)值。

(14)

dual EnKF方法共有4個參數(shù)，即集合成員個數(shù)N和3個超參數(shù)γE、γP、γQ。隨著集合成員個數(shù)N的增加，集合成員的代表性逐漸增強(qiáng)、分析誤差逐漸減小，但卻加重了計算負(fù)擔(dān)。因此，需要選擇一個合適的集合數(shù)既能滿足集合成員代表性的要求，又要考慮計算負(fù)擔(dān)問題[17]。超參數(shù)的大小直接影響著驅(qū)動數(shù)據(jù)和觀測數(shù)據(jù)誤差的大小。Moradkhani等[10]通過研究發(fā)現(xiàn)當(dāng)N=50、γE=γP=γQ=0.1時，模擬效果較好。因此，本文中的4個參數(shù)的取值與Moradkhani的研究結(jié)果保持一致。

3.3 dual E-G方法

目前大多數(shù)基于GLUE方法的水文模型不確定性研究結(jié)果表明：隨著覆蓋率的增高，平均帶寬反而增大[7,12,14]，無法提供有效的決策依據(jù)。

為了探索解決在提高區(qū)間覆蓋率的同時降低預(yù)測區(qū)間平均帶寬問題，將dual EnKF方法與GLUE方法2種水文模型不確定性研究方法進(jìn)行了耦合，得到了dual E-G方法。具體耦合方法如下：①首先按照GLUE方法產(chǎn)生似然函數(shù)大于閾值TH的M組有效參數(shù)組；②將每組模型參數(shù)產(chǎn)生的服從均值為參數(shù)本身，方差為參數(shù)的0.1倍的N個參數(shù)組作為每個成員的參數(shù)初始值；③然后按照dual EnKF方法進(jìn)行模型參數(shù)更新；④最后將得到的N個模擬流量系列的均值作為該組參數(shù)的最優(yōu)模擬系列，則M組有效參數(shù)組便得到M個模擬流量系列；⑤將M個模擬流量系列同一時刻的所有模擬流量Q(i,j)分別從小到大排序；⑥對于給定的置信水平α(0<α<1)，分別計算第j時刻實(shí)測流量的下界預(yù)測值QL(j)和上界預(yù)測值QU(j)。

4 實(shí)例計算與分析

為了分析耦合效果，以新安江模型為例，模型參數(shù)取值范圍見表2。對于給定的置信水平α(0<α<1)，采用3個評價指標(biāo)(覆蓋率、平均帶寬和平均不對稱性)評價了dual E-G方法在3個流域(甲河、天河、灞河)的預(yù)測區(qū)間性質(zhì)，并與GLUE方法進(jìn)行了對比，見表3、4。由表4可知，相比于GLUE方法，dual E-G方法在甲河、天河和灞河流域中計算得到的率定期的區(qū)間覆蓋率分別提高了28%、41%和29%，平均帶寬分別降低了8%、18%和6%，平均不對稱程度分別降低了62%、75%和44%；檢驗(yàn)期的區(qū)間覆蓋率分別提高了31%、53%和32%，平均帶寬分別降低了4%、1%和2%，平均不對稱程度分別降低了67%、80%和59%。綜上所述，dual E-G方法可以同時改善3種不確定性指標(biāo)。

表2 新安江模型參數(shù)取值范圍

表3 GLUE方法和dual E-G方法在不同流域上的預(yù)測區(qū)間性質(zhì)

表4 dual E-G方法相比于GLUE方法的改善程度 %

為了便于比較不同方法計算出來的預(yù)測區(qū)間的性質(zhì)，本文選用了甲河流域率定期中的1984年和檢驗(yàn)期中的1985年的預(yù)測區(qū)間圖來比較GLUE方法和dual E-G方法預(yù)測區(qū)間的優(yōu)良性，見圖2?？梢钥闯觯cGLUE方法相比，EnKF-G方法率定期模擬洪峰能力提高，檢驗(yàn)期的低水部分模擬效果更好一些，且預(yù)測區(qū)間改善比較明顯，不僅提高了區(qū)間覆蓋率、降低了區(qū)間平均不對稱程度，還在一定程度上降低了平均帶寬。

a)GLUE方法率定期

d)dual E-G方法檢驗(yàn)期

為了說明dual E-G方法降低不確定性的原理，以GLUE方法篩選出的新安江模型的第1組有效參數(shù)為例進(jìn)行分析。對于GLUE方法，直接采用第1組有效參數(shù)計算的模擬流量系列計算納西效率系數(shù)；對于dual E-G方法，將由第1組有效參數(shù)生成的服從均值為模型參數(shù)本身，方差為模型參數(shù)的0.1倍的50個參數(shù)組作為每個成員的初始參數(shù)組，然后再按照dual EnKF方法進(jìn)行參數(shù)向量更新，從而計算出50個成員的模擬流量系列，將集合均值作為各時刻流量的最優(yōu)預(yù)報值，組成最優(yōu)模擬系列，進(jìn)而計算納西效率系數(shù)，結(jié)果見表5。由表5可知，dual E-G方法模擬精度明顯高于GLUE方法，甲河流域率定期納西效率系數(shù)提高了0.26，達(dá)到了0.89；灞河流域檢驗(yàn)期納西效率系數(shù)提高了一倍以上，達(dá)到0.62。

表5 第1組有效參數(shù)模擬系列納西效率系數(shù)對比

GLUE方法能夠充分考慮模型參數(shù)的空間不確定性，篩選出大于似然閾值的有效參數(shù)組，但是由于未能對參數(shù)進(jìn)行實(shí)時更新，隨著時間的變化，模型參數(shù)的時間不確定性逐漸累積，導(dǎo)致模擬精度有所降低；dual EnKF方法能夠充分考慮模型參數(shù)的時間不確定性，對于給定位置的有效參數(shù)組能夠通過實(shí)時數(shù)據(jù)及誤差擾動進(jìn)行不斷更新，以達(dá)到最優(yōu)的模擬效果，但是無法體現(xiàn)下墊面變化導(dǎo)致的參數(shù)空間不確定性影響。dual E-G方法將GLUE方法考慮的模型參數(shù)空間不確定性與dual EnKF方法考慮的模型參數(shù)時間不確定性進(jìn)行了充分結(jié)合，使得模擬效果更加接近實(shí)際，從而表現(xiàn)出更高的模擬精度和更優(yōu)的預(yù)測區(qū)間。

5 結(jié)語

單純采用GLUE方法進(jìn)行水文模型不確定性分析時不可避免地會產(chǎn)生以增加平均帶寬為代價來提高區(qū)間覆蓋率或以降低區(qū)間覆蓋率為代價來減小平均帶寬的問題，主要原因是由于“異參同效現(xiàn)象”的存在導(dǎo)致最優(yōu)參數(shù)很難獲得。dual EnKF方法可以根據(jù)每個成員的狀態(tài)變量先驗(yàn)值，通過迭代更新使得每組參數(shù)都能夠盡可能達(dá)到最優(yōu)，從而提高每個成員的模擬效果。耦合后的dual E-G方法，能夠充分考慮模型參數(shù)的空間不確定性和時間不確定性，可以使得M組有效參數(shù)組都能夠盡可能達(dá)到最優(yōu)，不僅可以提高模型預(yù)報精度，還可以改善預(yù)測區(qū)間性質(zhì)，即可以同時改善區(qū)間覆蓋率、平均帶寬和平均不對稱程度。

以新安江模型為例，研究了dual E-G方法在漢江3個子流域上的適用性。該成果可為水文模型預(yù)報、洪水決策風(fēng)險分析等方面的深入研究提供有益參考。然而水文模型的不確定性是一個非常復(fù)雜的系統(tǒng)問題，dual E-G方法在不同模型、不同流域上的適用性仍有待進(jìn)一步研究。