水文模型不確定性研究方法耦合分析

鄭言峰

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

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

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

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

1 研究區(qū)概況

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

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

表1 研究流域資料

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

2 模型簡(jiǎn)介及評(píng)價(jià)指標(biāo)

2.1 新安江模型簡(jiǎn)介

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

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

新安江模型作為國(guó)內(nèi)常用的一種概念性水文模型，大部分模型參數(shù)都具有明顯的物理意義，不僅模擬精度高，而且對(duì)水文資料的要求比較低，對(duì)于日模型來(lái)說(shuō)，只需要流域日蒸發(fā)資料、日降雨資料以及日徑流資料即可。新安江模型共有15個(gè)模型參數(shù)。①K：流域蒸散發(fā)能力Ep與蒸發(fā)皿實(shí)測(cè)值EM之比，簡(jiǎn)稱蒸散發(fā)折算系數(shù)；②IMP：一般將整個(gè)流域面積分為透水面積與不透水面積，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)于線性水庫(kù)的個(gè)數(shù)或水庫(kù)的調(diào)節(jié)次數(shù)；KK：線性水庫(kù)的蓄泄系數(shù)，相當(dāng)于流域匯流時(shí)間的參數(shù)，具有時(shí)間因次。

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

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

(1)

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

2.3 預(yù)測(cè)區(qū)間評(píng)價(jià)指標(biāo)

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

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

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

(2)

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

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

(3)

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

3 不確定性研究方法

3.1 GLUE方法

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

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

3.2 dual EnKF方法

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

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

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

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

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

(4)

(5)

(6)

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

(7)

t=1時(shí)，

步驟五計(jì)算每個(gè)成員的模擬流量。

(8)

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

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

(14)

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

3.3 dual E-G方法

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

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

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

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

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

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

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

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

a)GLUE方法率定期

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

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

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

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

5 結(jié)語(yǔ)

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

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