梁冀雨，劉曙光，周正正，鐘桂輝，方 琦，甄億位

（同濟(jì)大學(xué)土木工程學(xué)院，上海 200092）

山洪災(zāi)害是山丘區(qū)由降雨、融雪等引起的洪水、泥石流河滑坡災(zāi)害的總稱(chēng)[1］，其導(dǎo)致泥沙超量補(bǔ)給現(xiàn)象壅高河道水位產(chǎn)生流量高估現(xiàn)象[2］，使流域出口斷面實(shí)測(cè)徑流大于響應(yīng)流域降雨的清水流量，破壞了流域水量平衡關(guān)系，增加了山丘區(qū)流域徑流模擬結(jié)果的輸入不確定性?，F(xiàn)有水文模型受限于其簡(jiǎn)化的物理過(guò)程、參數(shù)估算方法及資料觀測(cè)條件，在模擬復(fù)雜的山丘區(qū)小流域水文過(guò)程中也會(huì)產(chǎn)生來(lái)自模型輸入、模型參數(shù)及模型結(jié)構(gòu)的不確定性[3-4］。在校準(zhǔn)水文模型時(shí)，上述所有來(lái)源的不確定性將集中傳遞給模型參數(shù)[5-6］，降低了徑流模擬結(jié)果的準(zhǔn)確性并造成“異參同效”現(xiàn)象[4］。

貝葉斯方法通過(guò)模擬徑流的后驗(yàn)分布定量計(jì)算其不確定性，提供了較確定性模擬方法更加合理的系統(tǒng)描述，增加了徑流模擬結(jié)果的可靠性[7-8］。貝葉斯方法通常假定觀測(cè)徑流可分解為由水文模型輸出表示的確定性分量及代表殘差的隨機(jī)分量[9］。其中，隨機(jī)殘差的統(tǒng)計(jì)模型確定了似然函數(shù)及徑流預(yù)測(cè)分布的形式，直接影響著模型參數(shù)及貝葉斯徑流預(yù)測(cè)分布估計(jì)的魯棒性和準(zhǔn)確性[10］。模擬徑流殘差具有2個(gè)重要統(tǒng)計(jì)特征，即時(shí)間上的自相關(guān)性與異方差性（殘差方差隨流量增加）[11］。大量的研究證明[6，11-14］，濕潤(rùn)流域徑流殘差序列的自相關(guān)性可使用一階自回歸模型描述，而對(duì)于具有明顯豐枯變化的流域，則可利用Copula函數(shù)建立更加靈活的殘差自相關(guān)形式[6，11-14］。另一方面，對(duì)徑流殘差的異方差性的處理可分為2種思路：直接處理法和轉(zhuǎn)換方法。直接處理法通過(guò)建立殘差標(biāo)準(zhǔn)差與協(xié)變量（如流量）的函數(shù)關(guān)系，以此對(duì)殘差進(jìn)行標(biāo)準(zhǔn)化處理后，假定殘差為零均值、同分布[6，14-17］。轉(zhuǎn)換方法則是在對(duì)觀測(cè)和模擬徑流進(jìn)行函數(shù)變換后，直接假定轉(zhuǎn)換后的殘差服從均值為零的常方差正態(tài)分布[18］。

雖然貝葉斯徑流模擬方法能夠提升水文模型參數(shù)估計(jì)的魯棒性并給出徑流預(yù)測(cè)分布[19］，但是，目前的貝葉斯徑流模擬方法皆沒(méi)有考慮水文過(guò)程之外的不確定性來(lái)源[20］，特別是忽略了山丘區(qū)流域山洪災(zāi)害引起的流量放大作用，使其在山丘區(qū)流域應(yīng)用時(shí)仍存在以下不足：①假定的殘差分布無(wú)法描述山洪災(zāi)害導(dǎo)致的流量殘差極大值，影響了似然函數(shù)的合理性；②參數(shù)后驗(yàn)估計(jì)為補(bǔ)償山洪災(zāi)害導(dǎo)致的觀測(cè)徑流極大值而產(chǎn)生偏差[21］，降低了后驗(yàn)徑流預(yù)測(cè)分布的整體可靠性；③頻繁出現(xiàn)的山洪流量極大值增大了模擬殘差的方差，令徑流預(yù)測(cè)分布跨度過(guò)大、預(yù)測(cè)精度差。因此，在山洪頻發(fā)的山丘區(qū)流域，貝葉斯徑流模擬方法需要進(jìn)一步優(yōu)化和改進(jìn)。

本文以汶川地震后山洪災(zāi)害頻發(fā)的都江堰白沙河流域?yàn)檠芯繀^(qū)域，探究引發(fā)山洪災(zāi)害的臨界雨量，引入線性放大系數(shù)修正山洪災(zāi)害對(duì)徑流量的放大作用，耦合一階自回歸殘差模型與不同異方差轉(zhuǎn)換函數(shù)構(gòu)建了6種貝葉斯徑流模擬方案。選用適用于山區(qū)小流域的概念性水文模型GR4J[22］，模擬不同方案下流域出口日徑流過(guò)程。通過(guò)比較不同模擬方案的表現(xiàn)，分析在山區(qū)貝葉斯徑流模擬中考慮山洪災(zāi)害流量放大作用的必要性。

1 研究方法

1.1 模擬殘差模型

令t（t=1，2，…，n）時(shí)刻山洪災(zāi)害流量放大作用導(dǎo)致的附加洪水流量為QF，t，則該時(shí)刻的觀測(cè)流量Q?t可表示為

式中：QS，t為t時(shí)刻模擬流量；εt為t時(shí)刻的原始?xì)埐?。若假定QF，t是QS，t的線性函數(shù)，式（1）可改寫(xiě)為

式中：α為模擬流量線性放大系數(shù)，α≥1。傳統(tǒng)的自相關(guān)異方差殘差模型不考慮山洪對(duì)流量的放大作用，可視為α=1。進(jìn)一步假定山洪災(zāi)害觸發(fā)條件為前期累積雨量超過(guò)臨界雨量，有

式中：PC，t為t時(shí)刻的前期累積雨量；PF為相應(yīng)的山洪臨界雨量。為消除εt的異方差性，在計(jì)算殘差前，選定轉(zhuǎn)化函數(shù)Z將流量正態(tài)化，得到正態(tài)化殘差ηt，即

利用一階自回歸模型（AR1）描述ηt的自相關(guān)性，有

式中：?1為一階自回歸系數(shù)；yt為正態(tài)化流量殘差新息（innovation）。通常假定yt服從均值為零的同方差正態(tài)分布[14，17，23］，于是

其中σy為新息yt的標(biāo)準(zhǔn)差。式（1）至式（6）構(gòu)建了考慮山洪影響的徑流量殘差模型，其基本假定為流量殘差在去除異方差性、自相關(guān)性及山洪放大作用后服從均值為零的同方差正態(tài)分布。于是，殘差模型參數(shù)θε={α，?1，θT}，其中θT為轉(zhuǎn)換函數(shù)的參數(shù)。

1.2 似然函數(shù)

令水文模型參數(shù)為θH，根據(jù)貝葉斯原理，在給定輸入資料序列X?及觀測(cè)流量序列Q?時(shí)，θH、θε的后驗(yàn)分布可表示為

式中：p(θH，θε)為水文、殘差模型參數(shù)的聯(lián)合先驗(yàn)分布；p(Q?|θH，θε，X?)為似然函數(shù)，可看作殘差的聯(lián)合分布，即p(Q?|θH，θε，X?)=p(ε(θH，Q?，X?)|θε)。若給定異方差轉(zhuǎn)換函數(shù)Z，在式（1）至式（6）的假定下，可以得到似然函數(shù)

式中：Z′(Q?t)=dZ(Q?t)/dQ?t；fN為正態(tài)分布的概率密度函數(shù)。本文采用的轉(zhuǎn)換函數(shù)及其相應(yīng)的似然函數(shù)表達(dá)式見(jiàn)表1，其中A1、A2、A3、λ1、λ2為轉(zhuǎn)換函數(shù)參數(shù)。若λ2=0，則BC（Box-Cox）轉(zhuǎn)換變?yōu)長(zhǎng)G（對(duì)數(shù)）轉(zhuǎn)換。確定轉(zhuǎn)換函數(shù)后，似然函數(shù)p(Q?|θH，θε，X?)通過(guò)殘差新息構(gòu)建觀測(cè)數(shù)據(jù)與模型、殘差參數(shù)之間的聯(lián)系，并以此更新參數(shù)的先驗(yàn)分布。

表1 異方差轉(zhuǎn)換函數(shù)與似然函數(shù)Tab.1 Transformation functions and corresponding likelihood functions

1.3 徑流的預(yù)測(cè)分布

將式（8）代入式（7），利用馬爾科夫鏈蒙特卡羅（MCMC）抽樣方法，可以得到水文、殘差模型參數(shù)的后驗(yàn)分布p(θH，θε|X?，Q?)。本文選擇Vrugt[24］提出的 DREAM（Differential Evolution Adaptive Metropolis algorithm）算法作為抽樣算法，該方法采用自適應(yīng)隨機(jī)子空間采樣技術(shù)（即自適應(yīng)調(diào)整的交叉概率函數(shù)），并能在抽樣過(guò)程中去除無(wú)用的馬爾可夫鏈，提高了傳統(tǒng)MCMC方法的抽樣效率和解的精度，在水文模型的校準(zhǔn)方面得到了廣泛的應(yīng)用[24-26］。基于概念水文模型及殘差模型參數(shù)的后驗(yàn)分布，可以 生 成 模 擬 流 量 的 預(yù) 測(cè) 分 布Q?={Q?r，t；r=1，2，...，Nr，t=1，2，...，n}，其中Nr為采樣次數(shù)，n為樣本容量，其具體步驟如下：

（1）對(duì)t=1，2，...，n，從參數(shù)后驗(yàn)分布中抽取θH，r、θε，r，對(duì)給定的水文模型f，。計(jì)算

（2）根據(jù)式（1）-（6），計(jì)算新息標(biāo)準(zhǔn)差σy，r。

（3）生成t時(shí)刻的隨機(jī)殘差新息yr，t～N(0，σy，r)。

（4）根據(jù)式（5）計(jì)算ηr，t。

（5）根據(jù)式（4）對(duì)轉(zhuǎn)換流量進(jìn)行逆轉(zhuǎn)換，得到預(yù)測(cè)流量

完成Nr次采樣后，集合Q?就是徑流的預(yù)測(cè)分布，其每一個(gè)現(xiàn)實(shí)Q?r={Q?r，t；t=1，2，...，n}可以看作徑流過(guò)程的一個(gè)點(diǎn)估計(jì)。于是t時(shí)刻徑流的后驗(yàn)中位數(shù)估計(jì)Q?m，t即為流量預(yù)測(cè)分布Q?的橫截面Q?t={Q?r，t；r=1，2，...，Nr}的中位數(shù)，而Q?t的其余分位數(shù)可作為相應(yīng)置信水平下t時(shí)刻預(yù)測(cè)流量置信區(qū)間上下限。

2 研究區(qū)域與模擬方案設(shè)置

2.1 研究區(qū)域

以汶川地震后山洪災(zāi)害頻發(fā)的都江堰市白沙河流域?yàn)檠芯繀^(qū)域，其集水面積為363 km2。該流域位于亞熱帶濕潤(rùn)季風(fēng)氣候區(qū)，降雨充沛，多年平均降雨量約1 134 mm，降雨量年內(nèi)分配不均勻性顯著，5至9月降雨量占全年的80%。同時(shí)，降雨是該流域山洪災(zāi)害的主要誘發(fā)因素[27］。白沙河流域內(nèi)共有3個(gè)雨量站（楊柳坪站、大火地站和虹口站），其中楊柳坪站也是位于流域出口的水文站（見(jiàn)圖1）。選用白沙河流域2009—2018年的逐日降雨、蒸發(fā)及徑流序列進(jìn)行模擬。

圖1 白沙河流域概況Fig.1 Water system,gauging sites,and DEM in the Baisha River Basin

2.2 水文模型

選用Perrin等[22］建立的GR4J概念水文模型，該模型適用于山區(qū)小流域的日尺度徑流模擬[28］。GR4J模型含有1個(gè)產(chǎn)流單元以及1個(gè)匯流存儲(chǔ)單元，并使用2個(gè)單位線方程來(lái)描述降雨與徑流產(chǎn)生之間的時(shí)間間隔，能夠模擬流域從降水、蒸散發(fā)到產(chǎn)流、匯流的整體過(guò)程。模型僅含4個(gè)參數(shù)，分別為產(chǎn)流單元最大容量S1（mm）、地表水交換系數(shù)EX（mm·d-1）、匯 流 單 元 最 大 容 量S2（mm）及 單 位 線 系數(shù)U（d）。

2.3 山洪臨界雨量

根據(jù)白沙河流域內(nèi)31處山洪溝的監(jiān)測(cè)數(shù)據(jù)及流域山洪災(zāi)害調(diào)查報(bào)告，白沙河流域內(nèi)各山洪溝2009—2018年間發(fā)生山洪事件天數(shù)共347d，包含361次山洪事件。白沙河流域內(nèi)山洪溝主要分布于上游無(wú)人區(qū)，少量山洪溝發(fā)生洪水時(shí)僅會(huì)影響河道徑流過(guò)程，而不會(huì)構(gòu)成嚴(yán)重災(zāi)害。

楊柳坪站、大火地站和虹口站3站中，大火地站雨量序列存在缺測(cè)情況，為減少缺測(cè)資料影響，采用通過(guò)等權(quán)重的算術(shù)平均法估算流域平均逐日面雨量P?t，t=1，2，...，n。以山洪當(dāng)日的前DC日（含當(dāng)日）累積雨量為判別條件，計(jì)算日尺度山洪災(zāi)害判別準(zhǔn)確率

式中：AC為山洪災(zāi)害判別準(zhǔn)確率；F?為判別山洪發(fā)生的日期集合；F?為實(shí)測(cè)山洪發(fā)生的日期集合。選取DC=1，2，...，10d，計(jì)算流域前期累積面雨量

以累積雨量集合PC={PC，t；t=1，2，...，n}的百分位數(shù)作為臨界雨量PF，令Pe為相應(yīng)的百分比。根據(jù)田述軍等[27］的研究，汶川地震后，白沙河流域泥石流臨界雨量較震前有所減小，其中DC取1時(shí)，臨界雨量為51.2 mm，相應(yīng)Pe值為98.6。針對(duì)這一現(xiàn)象，以0.1為步長(zhǎng)，取Pe∈[90.0，99.9]，根據(jù)不同前期累積日數(shù)DC分組，計(jì)算山洪災(zāi)害判別準(zhǔn)確率AC，繪制每組AC樣本的箱型圖于圖2。

圖2 不同累積日數(shù)下山洪判別準(zhǔn)確率箱型圖Fig.2 Box plot of accuracy of flash flood identification

由圖2可知，山洪判別準(zhǔn)確率AC的均值隨前期累積日數(shù)DC的增加先增后減，其中當(dāng)DC取6d時(shí)，AC樣本具有最高的中位數(shù)0.63。選取AC均值最大的DC=5，6，7，8d共4組，繪制AC隨Pe的分布于圖3。如圖所示，當(dāng)DC=6d，Pe=96.4時(shí)，即山洪臨界雨量取前6日流域平均雨量18.9mm、累積雨量113.9mm時(shí)，AC得到最大值0.81。

圖3 山洪判別準(zhǔn)確率隨臨界雨量百分比分布圖Fig.3 Distribution of accuracy of flash flood identification with rainfall percentile

為進(jìn)一步確定山洪災(zāi)害對(duì)徑流的影響，根據(jù)圖3中的先驗(yàn)信息，本文選擇平均判別準(zhǔn)確率最高的DC=6d作為計(jì)算前期累積雨量及臨界雨量的日數(shù)，并將Pe作為殘差模型參數(shù)，在校準(zhǔn)模型時(shí)與其余參數(shù)一起估計(jì)。

2.4 徑流模擬方案

選取山洪災(zāi)害發(fā)生頻繁的2013年1月1日至2018年12月31日為模型校準(zhǔn)期，以2009年1月1日至2012年12月31日為模型驗(yàn)證期，模型校準(zhǔn)、驗(yàn)證期均以其第1年的資料對(duì)模型進(jìn)行預(yù)熱，調(diào)整模型的狀態(tài)參數(shù)S1、S2。為探究本文提出的殘差模型對(duì)貝葉斯徑流模擬方法表現(xiàn)的影響，采用6種模擬方案對(duì)白沙河流域楊柳坪站的流量序列進(jìn)行模擬，模擬方案的設(shè)置見(jiàn)表2、表3。

表2 方案設(shè)置及殘差模型參數(shù)范圍Tab.2 Simulation schemes and corresponding parameter ranges

表3 GR4J模型參數(shù)Tab.3 Parameter ranges of the GR4J model

表2給出了模擬方案的名稱(chēng)，使用的異方差轉(zhuǎn)換函數(shù)及殘差模型參數(shù)范圍，以后綴F表示方案采用考慮山洪影響的殘差模型。表3給出了GR4J的模型參數(shù)取值范圍。本文假定所有參數(shù)在其取值范圍內(nèi)服從均勻先驗(yàn)分布。同時(shí)，設(shè)定DREAM算法馬爾可夫鏈數(shù)為各方案參數(shù)總個(gè)數(shù)的2倍，參數(shù)后驗(yàn)分布樣本均取MCMC收斂后的2 000次抽樣。

2.5 徑流模擬結(jié)果評(píng)價(jià)指標(biāo)

根據(jù)選定的徑流模擬方案，為系統(tǒng)性比較不同異方差轉(zhuǎn)換方法及殘差模型的模擬表現(xiàn)，選取了不同的統(tǒng)計(jì)指標(biāo)來(lái)描述各徑流模擬方案的準(zhǔn)確性、預(yù)測(cè)分布可靠性、預(yù)測(cè)分布精度及殘差模型合理性。

2.5.1 模擬準(zhǔn)確性

通過(guò)徑流的后驗(yàn)中位數(shù)估計(jì)的納什效率系數(shù)（Nash-Sutcliffe model efficiency coefficient）描述徑流模擬的準(zhǔn)確性，其表達(dá)式

式中：-Q為觀測(cè)徑流序列的平均值。Nm的取值范圍為(-∞，1]，越接近1則模擬準(zhǔn)確度越高，反之模擬準(zhǔn)確度越差。

2.5.2 徑流預(yù)測(cè)分布可靠性

若徑流預(yù)測(cè)分布能捕捉觀測(cè)資料的分布特性，即徑流觀測(cè)值能看作從徑流預(yù)測(cè)分布中抽取的樣本，則稱(chēng)該預(yù)測(cè)分布是可靠的。采用預(yù)測(cè)QQ（predictive quantile-quantile plot）圖[18，29-30］可 視 化 表征徑流預(yù)測(cè)分布的可靠性，并采用Evin等[17］提出的可靠性指標(biāo)π對(duì)其定量計(jì)算，即

2.5.3 徑流預(yù)測(cè)分布精確性

徑流預(yù)測(cè)分布精確性指1.3節(jié)中由徑流預(yù)測(cè)分布分位數(shù)得到的置信區(qū)間寬度，可用精確性指標(biāo)Pc計(jì)算。

式中：std(Q?t)為徑流預(yù)測(cè)分布在t時(shí)刻的橫截面標(biāo)準(zhǔn)差。Pc≥0，Pc值越小，流量預(yù)測(cè)分布的精度越高、預(yù)測(cè)置信區(qū)間越窄。

3 徑流模擬結(jié)果分析

3.1 不同方案的徑流預(yù)測(cè)分布

根據(jù)表2、表3，利用貝葉斯方法對(duì)白沙河流域楊柳坪站的日徑流資料進(jìn)行模擬。不同方案去除模型預(yù)熱期后的驗(yàn)證期（2010年1月1日至2012年12月31日）徑流模擬結(jié)果如圖4所示。

圖4分別展示了考慮山洪災(zāi)害影響（圖4a、4c、4e）與未考慮山洪影響（圖4b、4d、4f）情況下模擬徑流分布的95%置信區(qū)間與其后驗(yàn)中位數(shù)估計(jì)，同時(shí)，圖4a與4b、圖4c與4d、圖4e與4f分別代表同一種異方差轉(zhuǎn)換函數(shù)。

比較圖4b、4d、4f可知，在未考慮山洪影響時(shí)，3種異方差轉(zhuǎn)換函數(shù)中，BC方案低流量模擬的95%置信區(qū)間最窄，然而在模擬高流量時(shí)出現(xiàn)了較大的不確定性。圖5給出了BC、BC_F方案下殘差模型參數(shù)λ2的后驗(yàn)核密度分布，如圖所示，BC方案λ2的后驗(yàn)眾數(shù)為-1.24，且其后驗(yàn)樣本均小于零。圖6給出了殘差模型后驗(yàn)眾數(shù)下（λ2=-1.24，A2=0.009，A1=0.017）BC與LG方案的轉(zhuǎn)換流量分布圖。如圖所示，相比于LG轉(zhuǎn)換，BC轉(zhuǎn)換大幅減小了高低流量差，這彌補(bǔ)了山洪災(zāi)害的流量放大作用，減小了殘差的整體方差，使BC方案的整體精度在未考慮山洪影響的徑流模擬方案中最優(yōu)。然而，BC方案按式（10）利用逆轉(zhuǎn)換生成流量預(yù)測(cè)分布時(shí)，有預(yù)測(cè)流量由于A2的后驗(yàn)估計(jì)趨于零，當(dāng)λ2＜0且QS，t較大時(shí)，λ2ηt+(QS，t+A2)λ2可能接近零，此時(shí)Q?t→01/λ2，導(dǎo)致流量預(yù)測(cè)分布頻繁采樣到不合理的大值，因此出現(xiàn)圖4d中的高流量區(qū)的高尖峰，降低了高流量區(qū)徑流預(yù)測(cè)分布的可靠性和精確性。LSH及LG方案的逆轉(zhuǎn)換函數(shù)較為穩(wěn)定[31］，不易產(chǎn)生不合理的高流量預(yù)測(cè)，但也無(wú)法彌補(bǔ)山洪災(zāi)害導(dǎo)致的殘差方差增大，導(dǎo)致其低流量區(qū)域的置信區(qū)間較寬。

圖4 楊柳坪站驗(yàn)證期不同方案徑流模擬結(jié)果Fig.4 Simulation performance of different schemes in validation period

圖5 λ2后驗(yàn)核密度分布Fig.5 Posterior kernel density ofλ2

圖6 LG、BC方案轉(zhuǎn)換流量對(duì)比Fig.6 Transformed runoff of LG and BC scheme

對(duì)比圖4a、4c、4e與圖4b、4d、4f，易見(jiàn)考慮山洪災(zāi)害對(duì)流量的放大作用后，模擬徑流的95%置信區(qū)間寬度皆明顯減小。其中，BC_F方案參數(shù)λ2的后驗(yàn)眾數(shù)為0.12（圖5），僅19%的λ2后驗(yàn)樣本小于零，顯著減小了流量預(yù)測(cè)分布采樣到異常值的概率，增加了BC逆轉(zhuǎn)換下高流量徑流預(yù)測(cè)分布的可靠性和精確性。

除轉(zhuǎn)換函數(shù)本身的性質(zhì)外，對(duì)應(yīng)于殘差標(biāo)準(zhǔn)差的殘差新息標(biāo)準(zhǔn)差σy對(duì)流量預(yù)測(cè)分布的精確性也有著直接影響，較大的σy值對(duì)應(yīng)較寬的徑流預(yù)測(cè)分布及較大的預(yù)測(cè)不確定性。表4給出了6種方案下σy的后驗(yàn)樣本的5%分位數(shù)、中位數(shù)及95%分位數(shù)。根據(jù)表4，采用考慮山洪影響的殘差模型后，LSH、BC、LG這3種方案的σy中位數(shù)分別減小84%、25%、23%，σy的90%置信區(qū)間寬度分別減小15%、90%、17%。

表4 不同方案后驗(yàn)σy分位數(shù)表Tab.4 Quantiles of posterior distribution ofσy

表5列舉了6種模擬方案的準(zhǔn)確性指標(biāo)Nm、可靠性指標(biāo)π及精確性指標(biāo)Pc。由表5可知，與前述一致，未考慮山洪影響時(shí)，BC方案在徑流預(yù)測(cè)分布精確性指標(biāo)Pc及準(zhǔn)確性指標(biāo)Nm上表現(xiàn)最優(yōu)。綜合考慮全部6種方案后，模擬準(zhǔn)確性最優(yōu)方案為L(zhǎng)SH_F，而預(yù)測(cè)分布可靠性及精確性的最優(yōu)方案皆為L(zhǎng)G_F。

表5 不同方案模擬結(jié)果評(píng)價(jià)指標(biāo)Tab.5 Performance metrics of different schemes

值得注意的是，除最優(yōu)方案外，亦可發(fā)現(xiàn)在采用考慮山洪影響的殘差模型后，徑流模擬結(jié)果評(píng)價(jià)指標(biāo)都有所改進(jìn)，其中，Nm平均提升19%，π平均降低32%，Pc平均降低62%，而對(duì)于不同的異方差轉(zhuǎn)換函數(shù)，LG方案3項(xiàng)評(píng)價(jià)指標(biāo)平均提升45%，LSH及BC方案分別平均提升43%和26%。

3.2 殘差新息的后驗(yàn)分布

以LG方案為例，分別繪制LG、LG_F方案標(biāo)準(zhǔn)化殘差新息v的直方圖于圖7、圖8。如圖7所示，與標(biāo)準(zhǔn)正態(tài)分布相比，LG方案的標(biāo)準(zhǔn)化殘差新息分布呈較強(qiáng)的正偏態(tài)，右側(cè)有長(zhǎng)尾并有較大的峰度。相較之下，圖8中LG_F方案的標(biāo)準(zhǔn)化殘差新息與標(biāo)準(zhǔn)正態(tài)分布更為一致。為進(jìn)一步比較，將LG和LG_F方案下v的均值、偏態(tài)系數(shù)及峰度系數(shù)列于表6。

圖7 LG方案標(biāo)準(zhǔn)化殘差新息直方圖Fig.7 Histogram of standardized innovation of LG scheme

圖8 LG_F方案標(biāo)準(zhǔn)化殘差新息直方圖Fig.8 Histogram of standardized innovation of LG_F scheme

表6 考慮山洪災(zāi)害影響前后標(biāo)準(zhǔn)化殘差新息v統(tǒng)計(jì)參數(shù)Tab.6 Statistics of standardized innovation of LG and LG_F schemes

由表6可知，相較于LG方案，LG_F方案下v的統(tǒng)計(jì)參數(shù)更接近殘差模型所假定的標(biāo)準(zhǔn)正態(tài)分布的統(tǒng)計(jì)參數(shù)。由于2種方案中v的標(biāo)準(zhǔn)差均為1，因此具有更大峰度系數(shù)的LG方案中v包含更多的特大值，根據(jù)式（3）-（6）可知，某些時(shí)間點(diǎn)的模擬徑流遠(yuǎn)小于觀測(cè)流量。分別記LG方案v中較大的前5%、3%為v5及v3。v5與山洪災(zāi)害發(fā)生時(shí)間重合的天數(shù)占v5總數(shù)的65%，v3相應(yīng)重合比例為83%，表明v出現(xiàn)特大值與山洪災(zāi)害發(fā)生具有較強(qiáng)的相關(guān)性。而這些特大值也是LG方案下v的分布出現(xiàn)較強(qiáng)正偏性的主要原因。同時(shí)，MCMC參數(shù)估計(jì)方法過(guò)多關(guān)注這些特大值[32］，導(dǎo)致參數(shù)的后驗(yàn)估計(jì)使模型模擬徑流QS整體偏大，造成低流量區(qū)域QS頻繁大于觀測(cè)流量Q?，使LG方案下v的均值小于零。

另一方面，LG_F方案通過(guò)流量線性放大系數(shù)α放大山洪災(zāi)害發(fā)生時(shí)的模擬徑流QS，減小了這些時(shí)間點(diǎn)上QS與觀測(cè)流量Q?的偏差；LG_F方案下v5、v3的均值較LG方案分別降低39%、57%。v的特大值的下降進(jìn)一步提升了MCMC參數(shù)估計(jì)的準(zhǔn)確性，降低的低流量區(qū)域負(fù)殘差的出現(xiàn)概率使v的均值由-0.342增至-0.017。綜上，考慮山洪影響的殘差模型通過(guò)識(shí)別并修正山洪災(zāi)害導(dǎo)致的轉(zhuǎn)換殘差特大值，使其更接近假定的同方差正態(tài)分布，LSH與BC方案下得到的結(jié)論與LG方案相似。

3.3 山洪修正參數(shù)的后驗(yàn)分析

相較于未考慮山洪災(zāi)害影響的自相關(guān)異方差殘差模型，本文通過(guò)引入?yún)?shù)Pe、α識(shí)別并修正山洪災(zāi)害對(duì)流量的放大作用。圖9、圖10給出了3種考慮山洪影響的殘差模型參數(shù)Pe、α的后驗(yàn)核密度分布。

圖9 各方案P e后驗(yàn)核密度分布Fig.9 Posterior kernel density of P e

圖10 LG_F方案α后驗(yàn)核密度分布Fig.10 Posterior kernel density ofα

根據(jù)圖9，LG_F方案Pe的后驗(yàn)分布集中于95至98之間，其最大后驗(yàn)估計(jì)為96.5，對(duì)應(yīng)的前6日臨界累積雨量為115.1mm。值得注意的是，該雨量并非傳統(tǒng)意義上的山洪臨界雨量，其含義是，當(dāng)前期雨量超過(guò)該值時(shí)，流域響應(yīng)徑流發(fā)生明顯增大，需要進(jìn)行修正。因此，該雨量應(yīng)稱(chēng)為“流量修正臨界雨量”。同時(shí)，根據(jù)圖10，LG_F方案α的后驗(yàn)眾數(shù)為6.0，在式（3）的假定下，即判斷山洪發(fā)生時(shí)，徑流量將被放大6倍，這一結(jié)果與胡凱衡等[33］得到的非線性流量修正系數(shù)結(jié)果相近。

LSH_F方案Pe、α的后驗(yàn)分布與LG_F方案相似，而B(niǎo)C_F方案的Pe后驗(yàn)眾數(shù)為97.7，判別驗(yàn)證山洪災(zāi)害發(fā)生的日數(shù)較LG_F、LSH_F方案減少約14d，山洪災(zāi)害識(shí)別準(zhǔn)確率AC為0.59。同時(shí)，BC_F方案α的后驗(yàn)眾數(shù)為2.8明顯小于LG_F、LSH_F方案。造成這一差異的一個(gè)可能原因是BC轉(zhuǎn)換函數(shù)Z2（表1）在λ2較小的情況下，減小了高流量觀測(cè)值與模擬值，因此平緩了山洪災(zāi)害導(dǎo)致的較大模擬殘差與普通殘差的區(qū)別，使殘差模型無(wú)法完全識(shí)別、修正某些場(chǎng)次的模擬殘差。同時(shí)，這種殘差平緩效應(yīng)也使得參數(shù)α的后驗(yàn)估計(jì)偏小。

4 結(jié)論

山洪災(zāi)害壅高河道水位帶來(lái)的流量放大作用影響了貝葉斯徑流模擬方法在山丘區(qū)流域的表現(xiàn)。以汶川地震后山洪災(zāi)害頻發(fā)的都江堰白沙河流域?yàn)檠芯繀^(qū)域，基于引發(fā)流域山洪災(zāi)害的臨界雨量，引入線性放大系數(shù)修正山洪災(zāi)害對(duì)流量的放大作用，耦合一階自回歸殘差模型與不同異方差轉(zhuǎn)換函數(shù)構(gòu)建了6種貝葉斯徑流模擬方案，利用GR4J模型模擬不同方案下流域出口日徑流過(guò)程，得到如下主要結(jié)論：識(shí)別并修正山洪災(zāi)害的流量放大作用后，①山丘區(qū)流域徑流殘差新息更接近殘差模型假定的同方差正態(tài)分布；②貝葉斯徑流模擬表現(xiàn)明顯提升，其中準(zhǔn)確性、可靠性及精確性指標(biāo)平均提升19%、32%、62%；③不同異方差轉(zhuǎn)換函數(shù)中，LG轉(zhuǎn)換得到了最優(yōu)的徑流預(yù)測(cè)分布可靠性及精確性，log-sinh轉(zhuǎn)換得到了徑流過(guò)程的最優(yōu)后驗(yàn)中位數(shù)估計(jì)，而B(niǎo)C轉(zhuǎn)換受限于其轉(zhuǎn)換函數(shù)特性，對(duì)山洪災(zāi)害的識(shí)別及修正效果較差，徑流模擬表現(xiàn)的提升有限。

提出的考慮山洪影響的貝葉斯徑流模擬方法理論基礎(chǔ)扎實(shí)，修正了山洪災(zāi)害對(duì)山丘區(qū)流域徑流量的放大作用，優(yōu)化了徑流模擬結(jié)果，可在汶川地震后山洪頻發(fā)的山丘區(qū)流域推廣應(yīng)用，為徑流預(yù)報(bào)提供更為準(zhǔn)確的科學(xué)依據(jù)。同時(shí)，為進(jìn)一步提升山丘區(qū)貝葉斯徑流模擬的表現(xiàn)，除臨界雨量、線性放大系數(shù)外，更加高效、準(zhǔn)確的山洪災(zāi)害識(shí)別及流量修正方法值得繼續(xù)深入研究。

作者貢獻(xiàn)聲明：

梁冀雨：數(shù)據(jù)收集整理、模型腳本編寫(xiě)、模擬結(jié)果分析及論文撰寫(xiě)。

劉曙光：提供研究思路及設(shè)計(jì)模擬實(shí)驗(yàn)方案。

周正正：數(shù)據(jù)校核及論文修訂

鐘桂輝：數(shù)據(jù)校核及論文修訂

方 琦：數(shù)據(jù)校核及論文修訂

甄億位：數(shù)據(jù)校核及論文修訂。