王 凱 蔣曉蕾

（淮河水利委員會(huì)水文局（信息中心） 蚌埠 233001 河海大學(xué) 南京 210098）

1 研究背景

水文預(yù)報(bào)是一種重要的防洪非工程措施，直接為防汛抗旱和水資源管理服務(wù)。目前廣泛使用的水文預(yù)報(bào)模型大多是確定性的，模型以確定預(yù)報(bào)值的形式輸出給用戶。實(shí)際上，由于水文過程的影響因素復(fù)雜多變，加之人類認(rèn)識(shí)水平的有限，使得水文預(yù)報(bào)過程中不可避免地存在著諸多的不確定性。目前，水文預(yù)報(bào)不確定性量化已逐漸成為了研究熱點(diǎn)，概率預(yù)報(bào)作為其重要的表現(xiàn)形式，已成為了水文預(yù)報(bào)的發(fā)展趨勢(shì)。

當(dāng)前國(guó)內(nèi)外關(guān)于洪水概率預(yù)報(bào)的研究主要包括兩類途徑：第一類是全要素耦合途徑，分別量化降雨—徑流過程各個(gè)環(huán)節(jié)的主要不確定性，如降雨輸入不確定性、模型結(jié)構(gòu)不確定性、模型參數(shù)不確定性等，并進(jìn)行耦合，實(shí)現(xiàn)概率預(yù)報(bào)。如Kavetski等采用“潛在變量”雨深乘子（stormdepthmultiplier）反映降雨輸入的不確定性，并將模型的敏感性參數(shù)隨機(jī)化，應(yīng)用馬爾可夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）采樣方法求解流量后驗(yàn)分布，提出貝葉斯總誤差分析（Bayesian total error analysis，BATEA）方法。在國(guó)內(nèi)，李明亮等基于層次貝葉斯模型（Bayesian hierarchical model），構(gòu)建聯(lián)合概率密度函數(shù)以考慮模型參數(shù)和降雨輸入不確定性，并采用MCMC方法進(jìn)行求解。梁忠民等借用抽站法原理推求降雨量的條件概率分布，進(jìn)而實(shí)現(xiàn)考慮輸入不確定性的洪水概率預(yù)報(bào)。全要素耦合途徑能夠溯源預(yù)報(bào)不確定性，但計(jì)算相對(duì)耗時(shí)，無法滿足實(shí)時(shí)預(yù)報(bào)的需要。

第二類是總誤差分析途徑，即從確定性預(yù)報(bào)結(jié)果入手，直接對(duì)預(yù)報(bào)不確定性進(jìn)行量化分析，推求預(yù)報(bào)量的分布函數(shù)，實(shí)現(xiàn)概率預(yù)報(bào)。Krzysztofowicz等提出的貝葉斯預(yù)報(bào)系統(tǒng)（Bayesian forecastingsystem，BFS）最具代表性，其中水文不確定性處理器（hydrologic uncertainty processor，HUP）在正態(tài)空間中對(duì)似然函數(shù)進(jìn)行了線性假定，為推求預(yù)報(bào)量后驗(yàn)分布的解析表達(dá)提供了可能。王善序詳細(xì)介紹了BFS的理論，并指出它可以綜合考慮預(yù)報(bào)過程的不確定性，不限定預(yù)報(bào)模型的內(nèi)部結(jié)構(gòu)。近年來，很多研究表明，不同流量量級(jí)的預(yù)報(bào)不確定性存在差異。王晶晶等在原始空間中開展預(yù)報(bào)誤差規(guī)律研究，發(fā)現(xiàn)不同流量量級(jí)的誤差服從不同的分布函數(shù)，為此采用極小熵確定各分布線型，采用極大熵進(jìn)行參數(shù)估計(jì)，進(jìn)而降低了福建池潭水庫管理的風(fēng)險(xiǎn)。Van Steenbergen等針對(duì)不同預(yù)見期、不同流量量級(jí)預(yù)報(bào)誤差統(tǒng)計(jì)規(guī)律的差異，采用頻率學(xué)方法，構(gòu)建了三維誤差矩陣，量化預(yù)報(bào)不確定性。

本文在第二類研究途徑框架下，結(jié)合研究區(qū)域（淮河王家壩斷面）預(yù)報(bào)誤差的統(tǒng)計(jì)特征，假定誤差均值隨流量量級(jí)呈分段線性變化，在貝葉斯預(yù)報(bào)系統(tǒng)（BFS）的基礎(chǔ)上，對(duì)其水文不確定性處理器（HUP）進(jìn)行改進(jìn)，提出改進(jìn)的貝葉斯概率預(yù)報(bào)系統(tǒng)（PCA-HUP），定量評(píng)估單個(gè)水文模型預(yù)報(bào)結(jié)果的可靠度，實(shí)現(xiàn)概率預(yù)報(bào)。

2 模型方法

考慮誤差異分布的洪水概率預(yù)報(bào)方法不涉及預(yù)報(bào)的中間環(huán)節(jié)，只對(duì)預(yù)報(bào)結(jié)果進(jìn)行分析，定量評(píng)估預(yù)報(bào)結(jié)果的不確定性，并在此基礎(chǔ)上實(shí)現(xiàn)洪水概率預(yù)報(bào)。在貝葉斯預(yù)報(bào)系統(tǒng)（BFS）的基礎(chǔ)上，對(duì)其水文不確定性處理器（HUP）進(jìn)行改進(jìn)，提出改進(jìn)的貝葉斯概率預(yù)報(bào)系統(tǒng)（PCA-HUP），定量評(píng)估單個(gè)水文模型預(yù)報(bào)結(jié)果的可靠度，實(shí)現(xiàn)概率預(yù)報(bào)。

2.1 HUP基本原理

水文不確定性處理器（Hydrologic UncertaintyProcessor，HUP）是貝葉斯概率預(yù)報(bào)系統(tǒng)（BFS）的主要組成部分，用以分析除降雨之外的其他所有不確定性。其特點(diǎn)是，不需要直接處理預(yù)報(bào)模型的結(jié)構(gòu)與參數(shù)，而是從預(yù)報(bào)結(jié)果入手，分析其與實(shí)測(cè)水文過程的誤差，再利用貝葉斯公式估計(jì)預(yù)報(bào)變量的概率分布，從而實(shí)現(xiàn)水文模型預(yù)報(bào)結(jié)果的不確定性分析及概率預(yù)報(bào)。其工作流程如圖1所示。

2.2 改進(jìn)的概率洪水預(yù)報(bào)PCA-HUP模型

為推求預(yù)報(bào)量后驗(yàn)分布的解析解，傳統(tǒng)HUP模型結(jié)合亞高斯模型，在正態(tài)空間中對(duì)先驗(yàn)分布式和似然函數(shù)式進(jìn)行線性假設(shè)，并采用最小二乘法對(duì)相關(guān)參數(shù)進(jìn)行估計(jì)。

然而，由于似然函數(shù)式的自變量之間存在明確的線性關(guān)系，必然導(dǎo)致回歸方程的多重共線性問題。若采用傳統(tǒng)最小二乘法進(jìn)行參數(shù)估計(jì)，會(huì)使得估計(jì)的回歸系數(shù)不唯一，也使得回歸方程不穩(wěn)定（原始數(shù)據(jù)的極小變化可造成參數(shù)估計(jì)值和標(biāo)準(zhǔn)差的明顯變化）。因此，項(xiàng)目研究結(jié)合主成分分析技術(shù)（Principal components analysis，PCA），對(duì)傳統(tǒng) HUP模型進(jìn)行改進(jìn)，提出PCA-HUP模型。

主成分回歸的基本思想：對(duì)原始回歸變量進(jìn)行主成分分析，將線性相關(guān)的自變量，轉(zhuǎn)化為線性無關(guān)的新的綜合變量，采用新的綜合變量建立模型回歸方程。

（1）主成分分析

設(shè) X=（X1，…Xp）T是 P 維隨機(jī)向量，均值 E（X）=μ，協(xié)方差陣 D（X）=Σ。

考慮它的線性變換：

用矩陣表示為：

圖1 水文不確定性處理器工作流程示意圖

由式（2）可以將 P 個(gè) X1，X2，…Xp轉(zhuǎn)化為 P 個(gè)新變量，Z1，Z2，…Zp若新變量 Z1，Z2，…Zp滿足下列條件：

1）Zi和 Zj相互獨(dú)立，i≠j，i，j=1，2，…，p；

2）Var（Z1）≥Var（Z2）≥…≥Var（Zp）；

則新變量 Z1，Z2，…Zp為 X1，X2，…Xp的 P 個(gè)主成分，且Z1，Z2，…Zp線性無關(guān)。

（2）主成分回歸

實(shí)際問題中不同的變量經(jīng)常具有不同的量綱，變量的量綱不同會(huì)使分析結(jié)果不合理，將變量進(jìn)行標(biāo)準(zhǔn)化處理可避免這種不合理的影響。

對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化首先要得出樣本標(biāo)準(zhǔn)差和樣本均值，記sj為xj樣本標(biāo)準(zhǔn)差，即是 xj的樣本均值，即原始數(shù)據(jù)的標(biāo)準(zhǔn)化變換為：

標(biāo)準(zhǔn)化后的數(shù)據(jù)矩陣為：

標(biāo)準(zhǔn)化后，X的相關(guān)系數(shù)陣也就是X的協(xié)方差陣（半正定矩陣）：

其中：

采用Lagrange乘子法求解，可以求得：

其中：λ1≥λ2≥…≥λp≥0 為 R 的特征值，a1，a2，…ap是相對(duì)應(yīng)的單位正交特征向量，ap=（a1pa2p，…app）T。

主成分回歸可以得到p個(gè)主成分，這p個(gè)主成分之間互相獨(dú)立，且方差呈遞減趨勢(shì)，所包含的自變量的信息也是遞減的。即主成分對(duì)因變量的貢獻(xiàn)率是遞減的，第i個(gè)主成分Zi的貢獻(xiàn)率可以用來表示。

在實(shí)際問題的分析時(shí)，由于主成分的貢獻(xiàn)率是遞減的，后面的主成分貢獻(xiàn)率有時(shí)會(huì)非常小，所以一般不選取P個(gè)主成分，而是根據(jù)累計(jì)貢獻(xiàn)率來確定主成分個(gè)數(shù)，即：前m個(gè)主成分的累計(jì)貢獻(xiàn)率達(dá)到0.85時(shí)，選取前m個(gè)主成分進(jìn)行回歸。則原始回歸問題轉(zhuǎn)化為以下回歸問題：

其中：E（εt）=0，Var（εt）=σ2，Cov（εi，εi）=0，（i≠j）

回歸模型的矩陣形式為：

采用最小二乘法估計(jì)參數(shù)矩陣B，根據(jù)式（1）和式（3）可以估計(jì)因變量矩陣Y與自變量矩陣之間X的回歸系數(shù)矩陣。

由此可見，主成分回歸模型是對(duì)普通的最小二乘估計(jì)的改進(jìn)，首先選取主成分，克服自變量間的多重共線性，然后對(duì)所選的主成分進(jìn)行線性回歸，進(jìn)而得到主成分回歸方程。

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

3.1 研究區(qū)概況

王家壩站系淮河上游總控制站，集水面積30630km2。上游干流河長(zhǎng)360km，河道比降0.5‰，年平均降水量800～1200mm，且降水年際變化大，時(shí)空分布不均勻。年降水量的60%集中在5～8月，以6、7兩月暴雨次數(shù)較多。產(chǎn)生暴雨的主要天氣系統(tǒng)是西南低渦、切變線、低壓槽和臺(tái)風(fēng)等。淮干上游及淮南山區(qū)一般是王家壩洪水的主要來源區(qū)。本項(xiàng)目主要針對(duì)息縣、潢川、班臺(tái)至王家壩區(qū)間（面積為7110km2）開展研究，見圖 2。

3.2 概率預(yù)報(bào)結(jié)果

本研究以新安江模型為確定性預(yù)報(bào)模型，在此基礎(chǔ)上，采用PCA-HUP模型分別對(duì)上述模型預(yù)報(bào)的可靠度進(jìn)行定量，并在淮河干流主要控制斷面王家壩進(jìn)行了實(shí)際應(yīng)用，采取滾動(dòng)預(yù)報(bào)方式實(shí)現(xiàn)了洪水概率預(yù)報(bào)。

將新安江模型的預(yù)報(bào)結(jié)果與實(shí)測(cè)資料輸入PCA-HUP模型中，其中20場(chǎng)洪水用于相關(guān)參數(shù)的率定，8場(chǎng)洪水用于模型驗(yàn)證。模型相關(guān)參數(shù)見表1。以置信度為90%（亦可采用其他置信度值）的預(yù)報(bào)區(qū)間為例，對(duì)概率預(yù)報(bào)結(jié)果進(jìn)行評(píng)估，同時(shí)，對(duì)流量分布函數(shù)的中位數(shù)Q50進(jìn)行分位數(shù)評(píng)價(jià)，率定期模型的模擬精度見表2。

由表2可知，PCA-HUP模型率定期模擬結(jié)果：預(yù)報(bào)區(qū)間（置信度為90%）覆蓋率較高，且離散度在0.2以內(nèi)。此外，將每一時(shí)刻預(yù)報(bào)量概率分布的中位數(shù)預(yù)報(bào)與實(shí)測(cè)流量進(jìn)行比較，確定性系數(shù)接近于1，洪峰誤差在1%以內(nèi)，說明中位數(shù)預(yù)報(bào)的精度非常高，且從不同預(yù)見期的模擬過程線中可以看出，隨著預(yù)見期的逐漸增大，區(qū)間離散度呈現(xiàn)出遞增的趨勢(shì)。

采用驗(yàn)證期8場(chǎng)洪水對(duì)概率預(yù)報(bào)模型進(jìn)行檢驗(yàn)，推求預(yù)報(bào)流量的概率分布，實(shí)現(xiàn)王家壩斷面的洪水概率預(yù)報(bào)。預(yù)報(bào)精度統(tǒng)計(jì)見表3，預(yù)報(bào)流量過程線以其中兩場(chǎng)為例，預(yù)報(bào)流量過程線如圖3～圖4所示。

由表3和2場(chǎng)洪水概率預(yù)報(bào)過程線可知，PCA-HUP模型（以新安江模型為確定性預(yù)報(bào)模型）提供的概率預(yù)報(bào)結(jié)果：預(yù)報(bào)區(qū)間（置信度為90%）覆蓋率在87%以上，且離散度在0.2以內(nèi)，說明在相對(duì)較小的區(qū)間寬度內(nèi)，預(yù)報(bào)區(qū)間仍然能夠覆蓋絕大多數(shù)實(shí)測(cè)數(shù)據(jù)，說明概率預(yù)報(bào)精度較高。此外，將每一時(shí)刻預(yù)報(bào)量概率分布的中位數(shù)預(yù)報(bào)與實(shí)測(cè)流量進(jìn)行比較，確定性系數(shù)接近于1，洪峰誤差在1%以內(nèi)，說明中位數(shù)預(yù)報(bào)的精度非常高，明顯高于新安江模型預(yù)報(bào)結(jié)果，充分體現(xiàn)了貝葉斯修正原理。從不同預(yù)見期的預(yù)報(bào)過程線中可以看出，隨著預(yù)見期的逐漸增大，區(qū)間寬度呈現(xiàn)出遞增的趨勢(shì)。

表2 PCA-HUP模型率定期模擬精度統(tǒng)計(jì)表（王家壩，Δt=2h）

圖2 王家壩斷面控制區(qū)域圖

表1 PCA-HUP模型參數(shù)表（王家壩）

圖3 a 19920505號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=2h）

圖3 b 19920505號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=6h）

圖3 c 19920505號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=12h）

圖4 a 19950707號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=2h）

圖4 b 19950707號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=6h）

圖4 c 19950707號(hào)洪水預(yù)報(bào)過程線圖（王家壩，Δt=12h）

表3 PCA-HUP模型驗(yàn)證期概率預(yù)報(bào)精度統(tǒng)計(jì)表（王家壩，Δt=2h）

4 結(jié)論

本文在貝葉斯預(yù)報(bào)系統(tǒng)（BFS）的基礎(chǔ)上，對(duì)其水文不確定性處理器（HUP）進(jìn)行改進(jìn)，提出改進(jìn)的貝葉斯概率預(yù)報(bào)系統(tǒng)（PCA-HUP），PCA-HUP模型是在正態(tài)空間中對(duì)預(yù)報(bào)誤差進(jìn)行分析，構(gòu)造線性似然函數(shù)，推求預(yù)報(bào)流量的后驗(yàn)分布，實(shí)現(xiàn)洪水概率預(yù)報(bào)。以淮河王家壩斷面的洪水預(yù)報(bào)為示例進(jìn)行了應(yīng)用研究，表明該方法不僅可以提供不同置信度的區(qū)間預(yù)報(bào)結(jié)果，還可以獲得精度更高的定值預(yù)報(bào)結(jié)果。

該方法可以與任何確定性的水文預(yù)報(bào)模型相耦合。作為示例，文中對(duì)相對(duì)誤差進(jìn)行了正態(tài)假設(shè)，但研究方法不局限于此，同樣適用于誤差服從其他分布的情況。限于研究區(qū)的預(yù)報(bào)誤差特征，本文僅考慮了誤差均值隨流量量級(jí)的變化規(guī)律，相同思路亦適用于誤差方差等分布參數(shù)隨流量變化的情況■