賀新月,曾獻(xiàn)奎,王 棟

(南京大學(xué)地球科學(xué)與工程學(xué)院,江蘇 南京 210023)

不確定性表示由于缺乏足夠信息而未能準(zhǔn)確描述或預(yù)測某事件的狀態(tài),其廣泛存在于工程、金融、地球科學(xué)等多個領(lǐng)域[1]。近幾十年來,流域水文過程的模擬預(yù)測及其不確定性分析已經(jīng)成為水文領(lǐng)域的研究熱點之一[2-3]。

水文系統(tǒng)的復(fù)雜性、資料的缺乏及人類認(rèn)識的局限性等原因使水文模型參數(shù)存在較大的不確定性[4]。隨著抽樣算法的不斷提升,馬爾可夫鏈蒙特卡洛(MCMC)模擬方法被廣泛用于分析復(fù)雜、高維的水文模型參數(shù)不確定性問題。Panday等[5]采用MCMC方法對喜馬拉雅地區(qū)Tamor河流域的融雪徑流模型進(jìn)行參數(shù)識別,結(jié)果表明MCMC方法可減少模型參數(shù)的不確定性,提高模型預(yù)測性能。魯帆等[6]以漢江流域丹江口水庫為例,將水文頻率分布線型的未知參數(shù)看作隨機(jī)變量,采用MCMC方法估計未知參數(shù)和設(shè)計洪水的后驗分布,結(jié)果表明MCMC模擬可有效推算出水文頻率分布線型參數(shù)及設(shè)計洪量的后驗分布。

MCMC方法通過似然函數(shù)描述模擬值與觀測值的擬合程度,確定模型參數(shù)的接受概率。在傳統(tǒng)的MCMC模擬中,通常假設(shè)模型模擬值與觀測值的差值(模擬殘差)無自相關(guān)性,將似然函數(shù)中的殘差協(xié)方差矩陣簡化為對角型結(jié)構(gòu)。Lu等[7]在采用最大似然估計法識別鈾反應(yīng)運(yùn)移柱體實驗?zāi)Ｐ蛥?shù)的研究中,發(fā)現(xiàn)模擬殘差時間序列(殘差序列)存在顯著的自相關(guān)性,采用簡化的殘差協(xié)方差矩陣計算似然函數(shù),將對模型表現(xiàn)的合理評價產(chǎn)生影響。Pang[8]為分析二元時間序列隨時間的動態(tài)變化情況,提出一種基于AR模型的線性多級模型來模擬該時間序列,通過修正傳統(tǒng)MCMC方法中的殘差協(xié)方差矩陣,有效識別模型參數(shù)。

本研究將融合AR模型的改進(jìn)MCMC方法(AR-MCMC)應(yīng)用于水文模擬的不確定性分析中,并以新疆提孜那甫河流域的融雪徑流模擬為例,對比評價傳統(tǒng)MCMC和AR-MCMC方法進(jìn)行參數(shù)不確定性分析的表現(xiàn)。

1 方 法

1.1 MCMC模擬

馬爾可夫鏈蒙特卡洛(MCMC)模擬方法通過構(gòu)建平穩(wěn)分布的馬爾可夫鏈搜索模型參數(shù)的概率分布空間[9]。在搜索過程中,不斷融合觀測信息D計算似然函數(shù),定量評價模型在參數(shù)θ下對觀測數(shù)據(jù)D的擬合表現(xiàn),使得馬爾可夫鏈從模型參數(shù)θ的先驗分布逐漸收斂至其后驗分布[10]。在MCMC模擬中一般采用式(1)計算似然函數(shù)。

(1)

式中:L(θ|D)——參數(shù)θ的高斯似然函數(shù);N——觀測數(shù)據(jù)的數(shù)量;f(θ)——參數(shù)取為θ時的模型輸出值;Σ——模擬殘差協(xié)方差矩陣。

MCMC模擬的一般步驟如下[11]:(a) 根據(jù)經(jīng)驗設(shè)置模型參數(shù)θ[θ1,θ2, …,θd]的先驗分布,d表示模型參數(shù)維數(shù)。(b) 基于某種分布函數(shù)在參數(shù)θ的先驗分布范圍內(nèi)抽取樣本θi(i=1,…,M),其中M為馬爾可夫鏈長度;運(yùn)行模型,計算參數(shù)樣本θi的L(θi|D)。(c) 對比L(θi|D)與上一步馬爾可夫鏈樣本θi-1的L(θi-1|D),計算θi的接受概率,更新馬爾可夫鏈。(d) 根據(jù)當(dāng)前所有樣本信息更新步驟(b)中的樣本抽取算法。(e) 重復(fù)步驟(b)～(d)至馬爾可夫鏈達(dá)到收斂標(biāo)準(zhǔn),將獲得的馬爾可夫鏈樣本θi和f(θi) 統(tǒng)計得到模型參數(shù)θ與模擬結(jié)果f(θ)的后驗分布。當(dāng)前,有多種抽樣算法可實現(xiàn)步驟(b)～(d),本文采用DREAMZS算法[12],該算法適用于復(fù)雜條件下的高維參數(shù)識別與不確定性分析。

1.2 AR模型

AR模型是一種被廣泛用于刻畫時間序列X(X1,X2,X3,…,Xn)自相關(guān)性的線性回歸模型。利用AR模型分析時間序列的相關(guān)性,需估算該序列的自相關(guān)函數(shù)值(ρ)與偏相關(guān)函數(shù)值(k)[14]。其中ρ表示序列中2個變量之間的總相關(guān)性,k表示2個變量排除了其他中間變量影響后的相關(guān)性。

根據(jù)k的衰減情況可確定AR模型的階數(shù)P[15]。若k在i個時間間隔后開始穩(wěn)定在(-2/n0.5, 2/n0.5)范圍內(nèi),則序列的自回歸階數(shù)P為i。確定模型階數(shù)P后,可依據(jù)式(2)計算第i個與第j個變量的自相關(guān)系數(shù)cij,其取值僅與i、j變量的時間間隔|i-j|有關(guān),可定義為時間間隔l的函數(shù)g(l)。

(2)

1.3 AR-MCMC方法

融合AR模型與MCMC方法的不確定性分析方法(AR-MCMC),采用AR模型描述殘差序列rN的相關(guān)性,修正MCMC似然函數(shù)中的殘差協(xié)方差矩陣。其具體步驟如下:(a) 構(gòu)建水文模型,運(yùn)行傳統(tǒng)MCMC進(jìn)行參數(shù)不確定性分析;(b) 計算似然函數(shù)最大時模擬殘差序列的ρ與k,確定自回歸階數(shù)P;(c) 重新運(yùn)行MCMC,針對每次迭代中的模擬殘差序列,建立其AR模型,修正殘差協(xié)方差矩陣結(jié)構(gòu)為Cek并計算似然函數(shù)L(θ|D),更新馬爾可夫鏈至總鏈長為M;(d) 對比傳統(tǒng)MCMC及AR-MCMC方法的不確定分析結(jié)果。

步驟(c)中,Cek的計算公式如下:

Cov(i,j)=E(ri-E(ri))E(rj-E(rj))=cijσiσj

(3)

式中Cov(i,j)——Cek中i行j列位置處的元素,表示殘差序列中第i個與第j個殘差之間的協(xié)方差;ri、rj——i時刻和j時刻的殘差(i,j=1,2,…,N);σi、σj——i時刻和j時刻觀測誤差的標(biāo)準(zhǔn)差。

1.4 模型評價方法

選擇模型邊緣似然值、預(yù)測區(qū)間的性質(zhì)以及對觀測序列的擬合精度來評價水文模型在傳統(tǒng)MCMC和AR-MCMC方法進(jìn)行不確定性分析的表現(xiàn)。邊緣似然值(L0)表示似然函數(shù)在先驗概率分布下的期望,可用于評價模型表現(xiàn)或估計模型權(quán)重[16]。預(yù)測區(qū)間性質(zhì)的評價指標(biāo)包括區(qū)間平均寬度W、對觀測數(shù)據(jù)的覆蓋率R、區(qū)間對稱性系數(shù)S。似然函數(shù)最大值對應(yīng)的納什系數(shù)E可用于評價模型對觀測序列的擬合精度。以上各項評價指標(biāo)計算公式如下:

(4)

(5)

(6)

(7)

(8)

本文選擇95%置信水平下的分布區(qū)間(2.5%～97.5%)作為預(yù)測區(qū)間。L0越大說明模型表現(xiàn)越好,S越接近1表明區(qū)間對稱性越好,E越大說明模型模擬序列與實測序列擬合效果越好。

2 案 例 研 究

2.1 研究區(qū)概況

研究區(qū)位于新疆提孜那甫河流域,流域地形復(fù)雜,海拔跨度1 472～6 352 m,面積約為5 518 km2。區(qū)域內(nèi)氣溫年內(nèi)變化大,蒸發(fā)強(qiáng)烈,降水稀少,但冰雪資源豐富,冰雪融水是河流的主要補(bǔ)給來源。研究對象為區(qū)域內(nèi)融雪及降水形成的地表徑流。

2.2 SRM模型

SRM模型是一種經(jīng)驗性融雪徑流模型,將徑流量分為3個部分:融雪徑流、降雨徑流以及原始徑流衰減后的徑流,模型基本結(jié)構(gòu)如下[17]:

(9)

式中:Q——日徑流量,m3/s;Cs——融雪徑流系數(shù);Cr——降水徑流系數(shù);a——度日因子,表示單位溫度、單位時間的融雪深度,cm/(℃·d);T——度日數(shù),℃·d;ΔT——利用溫度直減率對度日數(shù)進(jìn)行的校正值,℃;S——積雪覆蓋率;P——降雨形成的徑流深,cm;A——流域或流域分區(qū)的面積,km2;K——退水系數(shù);10 000/86 400——單位換算系數(shù)。其余符號含義見文獻(xiàn)[18]。

本次研究中,需要識別的模型參數(shù)包括a、Cs、Cr、K的擬合參數(shù)x和y。

2.3 數(shù)據(jù)

研究區(qū)高程數(shù)據(jù)采用NASA(national aeronautics and space administration)的ASTER數(shù)據(jù)庫的DEM數(shù)據(jù)。氣溫數(shù)據(jù)通過測站溫度和全球平均溫度直減率 (0.65℃/100 m)推求獲得。降水?dāng)?shù)據(jù)采用TRMM(tropical rainfall measuring mission)衛(wèi)星產(chǎn)品。另外利用MODIS(moderate-resolution imaging spectrometer) 8 d積雪合成產(chǎn)品(MOD10A2, h24v25) 進(jìn)行線性插值得到每日積雪覆蓋率。徑流數(shù)據(jù)來自研究區(qū)內(nèi)的水文觀測站——玉孜門勒克站。此次研究將2006年1月1日至12月31日的徑流觀測數(shù)據(jù)用于模型識別,2007年1月1日至12月31日的徑流觀測數(shù)據(jù)用于模型驗證。

3 結(jié)果與討論

3.1 殘差序列的P值

根據(jù)步驟中殘差序列ρ與k的計算結(jié)果,繪制兩者的變化曲線(圖1)。由圖1可見:ρ值在時間間隔為1 d時接近0.75,之后隨著時間間隔的增加呈現(xiàn)單邊遞減趨勢,最終接近0,說明模擬殘差序列具有顯著的自相關(guān)性,序列中2個殘差變量之間為正相關(guān)關(guān)系,并且這種正相關(guān)關(guān)系隨著時間間隔的增加逐漸減弱;根據(jù)AR模型計算得到k的截尾置信區(qū)間為(-0.105,0.105),結(jié)合k變化曲線可知,時間間隔大于1 d時k值穩(wěn)定在其截尾置信區(qū)間內(nèi),故此次研究實例中殘差序列的自回歸階數(shù)P為1。

圖1 殘差ρ與k變化曲線Fig.1 Variation curves of residuals’ ρ and k

3.2 模型參數(shù)識別

本次研究中需要識別的SRM模型未知參數(shù)包含x、y、Cs、Cr、a,各參數(shù)的先驗分布范圍依次為(0.85,1.2)、(-0.15, 0.3)、(0.1, 0.6)、(0.1, 0.6)、(0.05, 0.5)。根據(jù)模擬結(jié)果統(tǒng)計各參數(shù)經(jīng)MCMC方法與AR-MCMC方法識別后的后驗概率密度分布如圖2所示,似然函數(shù)最大值對應(yīng)的SRM參數(shù)取值如表1所示。采用MCMC與AR-MCMC進(jìn)行SRM模擬時,獲得的參數(shù)y和Cr后驗分布(如分布范圍、形狀)比較相似,而參數(shù)x、Cs和a的后驗分布差別較大。因此,MCMC與AR-MCMC方法采用不同結(jié)構(gòu)的殘差協(xié)方差矩陣計算似然函數(shù),導(dǎo)致2種方法獲得的參數(shù)識別結(jié)果不一致。

圖2 基于傳統(tǒng)MCMC和AR-MCMC的SRM參數(shù)后驗分布Fig.2 Posterior distributions of SRM parameters based on MCMC and AR-MCMC

表1 似然函數(shù)最大值對應(yīng)的模型參數(shù)值

表2 基于MCMC與AR-MCMC的識別與驗證結(jié)果

Table 2 Calibration and verification results based on MCMC and AR-MCMC

3.3 模型評價結(jié)果

采用不同結(jié)構(gòu)的殘差協(xié)方差矩陣(Cε不考慮殘差序列自相關(guān)性、Cek考慮殘差序列自相關(guān)性)計算的邊緣似然值L0如圖3所示。由圖3知:2種結(jié)構(gòu)的殘差協(xié)方差矩陣對應(yīng)的邊緣似然值均隨著樣本數(shù)的增多逐漸收斂,但Cek對應(yīng)的邊緣似然值更大。因此,在AR-MCMC模擬過程中采用Cek計算似然函數(shù),能獲得更大的邊緣似然值,即模型表現(xiàn)更好[18]。

圖3 不同殘差協(xié)方差矩陣對應(yīng)的邊緣似然值Fig.3 Marginal likelihoods for different residual covariance matrixes

根據(jù)MCMC與AR-MCMC 方法進(jìn)行SRM模型不確定性分析的徑流量模擬結(jié)果,各項評價指標(biāo)的計算結(jié)果見表2,融雪徑流預(yù)測區(qū)間繪圖見圖4。結(jié)合圖表分析:在識別期,AR-MCMC方法比MCMC方法獲得的流量預(yù)測區(qū)間平均寬度更大,分別為18.67 m3/s和12.31 m3/s;預(yù)測區(qū)間對實測數(shù)據(jù)的覆蓋率更高,分別為74.79%和59.73%;預(yù)測區(qū)間對稱性表現(xiàn)更好,對稱性系數(shù)分別為0.87和0.66。驗證期MCMC方法和AR-MCMC方法的評價結(jié)果與識別期一致,因此AR-MCMC方法較MCMC方法能推求出性質(zhì)更優(yōu)的預(yù)測區(qū)間。另外,采用傳統(tǒng)MCMC方法進(jìn)行SRM模擬時,識別期與模擬期的納什效率系數(shù)分別為0.84、0.87。而采用AR-MCMC方法對應(yīng)識別期與模擬期的納什效率系數(shù)分別為0.86、0.89。由此可見,采用AR-MCMC方法進(jìn)行融雪徑流模擬具有更優(yōu)的擬合效果。

圖4 基于傳統(tǒng)MCMC方法和AR-MCMC方法的SRM在識別期與驗證期的徑流模擬表現(xiàn)Fig.4 Runoff simulation performances of SRM based on traditional MCMC and AR-MCMC during calibration and verification periods

綜上,采用AR-MCMC方法進(jìn)行SRM模擬時,通過考慮殘差序列之間的相關(guān)關(guān)系,改進(jìn)似然函數(shù)中的殘差協(xié)方差矩陣,更新計算似然函數(shù),能夠更好地進(jìn)行模型參數(shù)不確定性分析。因此,AR-MCMC方法能夠更好地對研究區(qū)融雪徑流過程進(jìn)行模擬預(yù)測。

4 結(jié) 語

在傳統(tǒng)的水文模型不確定性分析中通常忽略模擬殘差序列的自相關(guān)性,將似然函數(shù)中的殘差協(xié)方差矩陣簡化為對角型結(jié)構(gòu),因此可通過修正殘差協(xié)方差矩陣進(jìn)一步提高模型不確定性分析效果。本文將融合AR模型的MCMC方法(即AR-MCMC方法)應(yīng)用于水文模型不確定性分析中,利用AR模型刻畫水文模型殘差序列的自相關(guān)性,修正殘差協(xié)方差矩陣,更新計算MCMC方法中的似然函數(shù)。

通過新疆提孜那甫河流域融雪徑流模擬的案例研究,分別利用傳統(tǒng)MCMC方法與AR-MCMC方法進(jìn)行融雪徑流模擬及預(yù)測的不確定性分析,結(jié)果發(fā)現(xiàn),SRM模擬的殘差序列具有顯著的自相關(guān)性,其對應(yīng)的自回歸階數(shù)P為1。采用AR模型修正似然函數(shù)中殘差協(xié)方差矩陣后,得到了更大的模型邊緣似然值,即提高了模型表現(xiàn)。此外,與傳統(tǒng)MCMC方法相比,利用AR-MCMC方法進(jìn)行SRM模擬及不確定性分析所獲得的徑流量預(yù)測區(qū)間對觀測數(shù)據(jù)的包含率更高、區(qū)間對稱性更好,似然函數(shù)最大值對應(yīng)模型的融雪徑流擬合效果更好。