凍 芳 芳,劉 沁 薇,王 淑 梅,劉 建 偉
(1.江西省水文局,江西 南昌 330000; 2.江西省上饒市水文局,江西 上饒 335100)
水文模型應用于洪水預報時,除了水文過程本身固有的自然不確定性外,還存在模型輸入、模型參數和模型結構的不確定性,這些不確定性的存在必將導致洪水預報成果存在不確定性[1-3]。關于如何定量描述洪水預報成果不確定性的研究一直是水文領域的熱點。縱觀國內外關于洪水預報不確定性分析的研究,大體上可以分為基于不確定性要素耦合途徑和基于預報總誤差分析途徑兩類[4]。不確定性要素耦合途徑的實質是分析預報過程中各環節的主要不確定性因子,估計其概率分布,再將這些不確定性耦合到洪水預報模型中,從而量化預報成果的不確定性[5-8];預報總誤差分析途徑不直接處理模型輸入數據、模型結構和參數的不確定性,代之以處理其綜合誤差(即終端輸出誤差),算法上不依賴于具體的水文預報模型,可以與任何模型相耦合,最終以后驗分布的形式直接提供模型輸出變量概率分布的估計,量化洪水預報成果的不確定性[4,9-10]。
水文不確定性處理器(HUP)是預報總誤差分析方法的一種,被廣泛應用于洪水預報不確定性分析中[11-14]。針對采用最小二乘法估計HUP模型參數過程中存在的多重共線性問題,Yao等[15]提出采用主成分分析技術(PCA)對傳統HUP模型進行改進,避免了HUP模型中似然函數求解過程中的多元共線性問題。
本文以江西省信江流域梅港水文站為研究對象,采用PCA-HUP模型分析洪水預報不確定性,在提供類似于傳統定值預報結果的同時,可實現預報成果的不確定性定量評估,為進一步提升信江流域防洪減災決策的科學水平提供技術支撐。
設H0為預報時刻已知的實測流量,Hn,Sn(n=0,1,2,…,N)分別為待預報的實測流量、確定性水文模型的預報流量,估計值分別用hn,sn表示,N為預見期的長度。根據貝葉斯原理及任意的觀測值H0=h0,在Sn=sn的條件下,可求得Hn的后驗密度函數[10-11]:
(1)
式中:hn為Hn的估計值;sn為Sn的估計值;h0為預報時刻已知的實測流量;fn為Hn的似然函數;gn為先驗密度函數。
預報量后驗分布函數的形式復雜,為了簡化計算及推求預報量后驗分布的解析解,一般先通過正態分位數轉換將非正態的流量序列分布轉換為亞高斯分布,在正態空間中對先驗分布式和似然函數式進行線性假設,并采用最小二乘法對相關參數進行估計[11]。
令Q表示標準正態分布,q表示相應的標準正態密度函數,則Hn與Sn轉換后的正態分位數可分別表示為
Wn=Q-1(Γ(Hn)),n=0,1,2,…,N
(2)
(3)

在獲得Hn和Sn的正態分位數Wn和Xn后,通過假定轉化空間中的流量過程服從一階馬爾可夫過程的正態-線性關系,構造了如下的先驗分布函數:
Wn=cWn-1+G
(4)
式中:c為參數,G為不依賴于Wn-1的殘差系列,且服從N(0,1-c2)的正態分布。
通過假定正態空間中的各變量Xn,Wn,W0服從正態-線性關系,構造了如下的似然函數:
Xn=anWn+dnW0+bn+Θn
(5)

在正態空間中,進一步推出Wn的后驗密度函數:
(6)
(7)
通過下式可獲得原始空間中預報流量的分布函數為
(8)
根據公式(8)預報流量的分布函數,可計算任一分位點的預報值,如50%概率對應的中位數預報及90%置信度對應的預報區間(5%概率和95%概率對應的分位數組成的區間)。
由先驗分布(式(4))可以看出,似然函數(式(5))中的自變量Wn和W0之間存在明確的線性關系,會導致公式(5)存在多重共線性問題。采用傳統最小二乘法進行參數估計,會使得估計的回歸系數不唯一,導致回歸方程不穩定。因此,Yao等[15]采用主成分分析技術(PCA)對傳統HUP模型參數估計方式進行了改進。
主成分回歸的基本思想是對原始回歸變量進行主成分分析,將線性相關的自變量轉化為線性無關的新綜合變量,采用新綜合變量建立模型回歸方程。針對變量Wn和W0。考慮它的線性變換:
(9)

Xn=bn+b1Z1+b2Z2
(10)
由于各個主成分之間是相互獨立的,為此可采用最小二乘法估計公式(10) 中的參數bn,b1和b2,進而克服了自變量間的多重共線性問題。
將公式(9)帶入到公式(10)中,公式(5)可進一步表示為
(11)
即,an=b1a11+b2a12,dn=b1a21+b2a22
梅港水文站為信江流域主要控制站,位于信江下游余干縣梅港鄉。梅港水文站集水面積15 535 km2,占信江總流域面積的97.4%。1952年4月設站,有流量、水位、降雨量、蒸發量、含沙量等測驗項目。實測最高水位為1998年的29.84 m(吳淞基面),實測最大流量為2010年的13 800 m3/s。本研究中,采用2012~2019年共10場洪水資料對PCA-HUP模型進行率定驗證,其中2012~2017年的8場洪水用于模型率定,而2019年的2場洪水用于模型驗證。場次洪水的確定性預報值是通過新安江模型預報獲得。
圖1給出了實測系列和模型預報系列的經驗累積分布函數及三參數對數威布爾分布函數的擬合效果圖。從圖1可以看出,三參數對數威布爾分布函數能很好地擬合實測系列和模型預報系列。

圖1 實測系列和模型預報系列的擬合效果Fig.1 Fitting effect of the measured series and prediction model series
表1給出了基于PCA-HUP模型的率定期8場洪水的分析結果,并與原始預報成果進行了對比分析。其中,PCA-HUP模型的確定性系數和洪峰相對誤差指標值是基于PCA-HUP模型提供的中位數(Q50)預報結果計算獲得,而離散度和覆蓋率指標值是基于PCA-HUP模型提供的90%置信度下的預報區間計算獲得。

表1 基于PCA-HUP模型的預報評估結果(率定期)Tab.1 Evaluation results of flood forecast based on PCA-HUP Model(rate period)
從表1中的中位數(Q50)預報結果來看,就確定性系數指標而言,8場洪水的確定性系數指標值均為0.99,高于原始預報的平均確定性系數為0.92;就洪峰相對誤差指標而言,8場洪峰的相對誤差絕對值的平均值為0.40%,遠小于原始預報的6.23%。在分析的8場洪水中,經PCA-HUP模型處理后,8場洪水的確定性系數均增大,且洪峰相對誤差均減小,表明經PCA-HUP模型處理后,洪水預報精度整體上有明顯的提高。從90%置信度下的預報區間來看,就離散度指標而言,其值在0.20~0.21之間,平均值為0.20;就區間覆蓋率指標而言,其值在91%~100%之間,平均值達97%。在離散度較小情況下,保證了較高的覆蓋率。
圖2給出了率定期2場洪水的實測值、原始預報、概率預報模型的中位數預報(Q50)及90%置信度下的區間預報成果。從圖2中可以看出,90%置信度的預報區間很好地描述了洪水預報的不確定性,且覆蓋了絕大多數的實測點據。

圖2 率定期2場洪水的PCA-HUP模型預報成果Fig.2 Forecast results of PCA-HUP model of two floods in rate period
表2給出了基于PCA-HUP模型的驗證期2場洪水的分析結果,并與原始預報成果進行了對比分析。從表中可以看出,就中位數(Q50)的預報結果而言,2場洪水的確定性系數值均為0.99,高于原始預報的平均確定性系數0.91;而2場洪峰的相對誤差絕對值的平均值為0.28%,遠小于原始預報的7.15%。就90%置信度下的預報區間而言,2場洪水的離散度均值為0.20;而區間覆蓋率的均值為0.98。與率定期結果相似,在離散度較小情況下,保證了較高的覆蓋率。

表2 基于PCA-HUP模型的預報評估結果(驗證期)Tab.2 Evaluation results of flood forecast based on PCA-HUP Model(verification period)
圖3給出了驗證期2場洪水的實測值、原始預報、概率預報模型的中位數預報(Q50)及90%置信度下的區間預報成果。從圖中可以看出,90%置信度的預報區間同樣很好地描述了驗證期2場洪水的預報不確定性,且覆蓋了絕大多數的實測點據。

圖3 驗證期2場洪水的PCA-HUP模型預報成果Fig.3 Forecast results of PCA-HUP model of two verification periodic floods
以信江流域梅港站為例,采用PCA-HUP模型對洪水預報不確定性進行了分析,在提供中位數(Q50)定值預報結果的同時,也提供了90%置信度下的預報區間成果。主要結論如下:
(1) 基于確定性系數指標和洪峰相對誤差指標,評估了原始確定性模型預報精度和PCA-HUP模型的中位數(Q50)預報精度。結果表明:無論是在率定期還是在驗證期,PCA-HUP模型的確定性系數平均值均高于原始預報的確定性系數平均值。經PCA-HUP模型處理后,洪水預報精度整體上有明顯的提高。
(2) 基于離散度和覆蓋率指標,評估了PCA-HUP模型計算的90%置信度預報區間的可靠性。結果表明:率定期和驗證期10場洪水的覆蓋率均在90%以上,且離散度均在0.20左右,具有較高的可靠性。