水庫抗暴雨能力概率預報的降雨情景構建方法

吳志勇，汪瑛琪，何 海，李 源，孫昭敏

(河海大學水文水資源學院，江蘇 南京 210098)

我國現有中小型水庫9萬多座。中小型水庫集水面積小，洪水匯集快，預報預見期短[1-5]，安全運行風險突出。水庫抗暴雨能力預報是實現中小型水庫安全運行的有效方法。該方法根據水庫集水區下墊面和水庫調度的動態情況，通過預演不同降水總量下水庫的水位變化過程，實現對預報水庫剩余防洪庫容所能容納的降水量的預報[6-10]。為防汛部門精準實施水庫調度方案，及時發布水庫預警，支撐水庫預報、預警、預演和預案工作，對保障我國中小型水庫防洪安全具有重要意義[10-12]。

當前水庫抗暴雨能力的計算主要通過設定不同等級降水總量，結合水文模型與水庫調度方案，得到不同降水總量與水庫最高水位的關系曲線；在此基礎上，依據水庫設計水位，反推水庫可容納降水總量，即該水庫當前狀態的抗暴雨能力[11-14]。該方法已在廣東、福建、江西和遼寧等地區廣泛應用。然而，傳統抗暴雨能力預報方法中降雨輸入的時程分配大多基于一場典型暴雨[14]，通過同倍比縮放得到不同降水總量下的降雨輸入，未能考慮到降水總量、降雨歷時和降雨時程分配之間的相關性及對產匯流的影響，使得水庫抗暴雨能力計算結果具有較大的不確定性。為此，本文在水庫抗暴雨能力計算中，引入Copula函數，構建考慮降水總量、降雨歷時和降雨時程分配相關性的降雨情景，生成同一降水總量下不同降雨歷時及時程分配概率的降雨情景，將其應用到水庫抗暴雨能力計算中，以實現水庫抗暴雨能力的概率預報，提高水庫抗暴雨能力的可應用性。

1 研究區概況與方法

1.1 研究區概況

選取位于江蘇省新沂市東部淋頭河上游的阿湖水庫以上流域作為研究區域(圖1)。該水庫是一座以防洪為主的中型水庫，總庫容4 032萬 m3，汛限水位25.50 m，設計水位27.73 m。水庫集水面積193 km2，位于南北氣候過渡帶和海陸交接帶，多年平均降水量為875.0 mm，年最大降水量為1 318.0 mm，汛期降水占年降水量60%左右。選取阿湖水庫水文站點2010—2020年汛期(5—9月)的逐時降水量資料開展研究。

圖1 阿湖水庫以上流域及水系分布

1.2 研究方法

水庫抗暴雨能力是指在當前水庫集水區下墊面狀態和水庫調度方式下，在一段時間內水庫剩余防洪庫容所能容納的降水量。設定不同等級的降水總量，輸入水文模型進行產匯流計算和水庫調洪演算，建立不同等級降雨與水庫最高水位的關系，通過設計水位反推水庫當前抗暴雨能力。

傳統水庫抗暴雨能力計算方法中，各級降水總量下僅考慮一種確定的降雨時程分配方案，例如典型暴雨同倍比縮放或均勻分配，忽略了降水總量、降雨歷時和降雨時程分配的相關性，以及在同一降水總量等級下不同歷時降雨情景的不確定性。

由于在降水總量一定的情況下，不同歷時及雨型對產匯流有較大影響，進而影響水庫抗暴雨能力的計算結果。因此，本文提出考慮不同降水總量下不同降雨歷時概率的降雨情景構建方法，并基于暴雨洪水同頻率假設，實現水庫抗暴雨能力的概率預報，提高水庫抗暴雨能力的可用性。該方法主要包括3個步驟：①從汛期連續的時段降雨資料中，劃分出獨立降雨事件，提取每場降雨事件的降雨特性(降水總量、降雨歷時和降雨時程分配)；②基于Copula函數構建降水總量和降雨歷時的聯合概率分布；③確定降水總量與降雨歷時分級下的代表雨型。

1.2.1降雨場次挑選

本文通過設定降雨閾值(最小降水總量)與降雨間歇時間[15](兩場降雨事件之間的無雨時間間隔)，從連續的時段降雨資料中劃分出獨立降雨事件開展研究。

1.2.2降水總量與降雨歷時的聯合分布

Copula函數是一種描述變量間相關性的連接函數，已被眾多學者用來描述降水總量與降雨歷時的相依結構[16-20]。參照Sadegh等[21]的相關研究，本文選用26種Copula函數作為備選函數，求解降水總量和降雨歷時的聯合分布，以得到同一降水總量下不同降雨歷時的概率。主要步驟如下：

a.選擇常用的伽馬(Gamma)分布、廣義極值(generalized extreme value，GEA)分布、廣義帕累托(generalized Pareto，GPA)分布、對數邏輯(Loglogistic)分布、對數正態(Lognormal)分布和威布爾(Weibull)分布作為降水總量與降雨歷時的備選邊緣分布。采用極大似然估計(maximum likelihood estimate，MLE)方法估計邊緣分布的參數，依據貝葉斯信息準則(Bayesian information criterion，BIC)和赤池信息準則(Akaike information criterion，AIC)的最小統計量選擇降水總量R與降雨歷時T的最優邊緣分布FR(R)和FT(T)，采用Kolmogorov-Smirnov(K-S)檢驗方法進行擬合優度檢驗。

b.將優選的邊緣分布帶入Copula備選函數中，分別采用MLE和基于拉丁超立方體抽樣的馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)方法[22]估計Copula參數，基于極大似然值(maximum likelihood values，MLV)最大、AIC和BIC最小準則比較兩種參數估計方法，并確定最優Copula函數F(R,T)=C(FR(R),FT(T))，F(R,T)表示降水總量R和降雨歷時T的聯合分布函數。通過均方根誤差RMSE、納什效率系數NSE，并結合Q-Q圖(quantile-quantile plot)檢驗最優的Copula函數。RMSE、NSE和Q-Q圖表示構建Copula聯結的雙變量概率與經驗聯合概率的擬合程度，當擬合效果好時，RMSE為0，NSE為1，Q-Q圖中散點落在45°對角線附近。

1.2.3降水總量和降雨歷時分級下代表雨型的確定

一場降雨事件中，降雨強度是隨時間變化的，即存在降雨時程分配。不同降雨事件具有不同的降雨時程分配，且受降水總量和降雨歷時的影響。由于降雨歷時可以反映降雨的集中程度，從某種程度上亦能代表降雨的時程分配。當降水總量一定時，若降雨歷時短，則降雨集中且強度大；若降雨歷時長，則降雨分散且均勻。因此，基于已劃分的降雨場次，對降雨時程分配進行分類，并結合降水總量和歷時的分級，采用每種分級下頻率最高的雨型作為代表雨型。結合第1.2.2節中生成的降水總量與降雨歷時，即可生成完整的降雨情景。

為了消除降雨事件中其他降雨特性對時程分配分類的影響，使降雨隨時間的變化是唯一影響不同時程分配的因素，首先要對降雨時程分配進行無量綱化處理[23]：

τ=t/T

(1)

dτ=Dt/R

(2)

式中：τ為時刻t的無量綱時間，τ∈(0,1]；T為一場降雨事件的降雨歷時；dτ為無量綱累計降水總量，即降水總量隨時間的累計百分比，dτ∈(0,1]；Dt為時刻t的累計降水總量。

其次，將每場降雨事件的無量綱降水總量累計曲線均分成K個時段(K值為研究區平均降雨歷時)。此時，無量綱時間τm=m/K(m=1,2,…,K)對應的無量綱累計降水總量值為dm。將每場降雨的K個dm值作為其K個屬性輸入到聚類方法中，以雨型劃分。采用基于歐幾里得距離進行分類的K均值聚類算法計算不同時程分配的K個屬性值間的歐幾里得距離，并通過迭代算法最終得到不同類別的時程分配。同一類別內的雨型數據相似度較高，而不同類別間的數據相似度較低。最后，采用列聯法統計分析在降水總量和歷時分級下各雨型的頻率，頻率最高的雨型為該分級下的代表雨型。將代表雨型與第1.2.2節中的降水總量、降雨歷時及其概率相結合，即可得到在同一降水總量下不同概率的降雨情景。

2 實例驗證

2.1 降雨場次提取

利用阿湖水庫報汛站汛期降雨數據提取降雨事件的準則為：①場次降水總量大于或等于20.0 mm，若該降水總量選取得過大，會導致降雨場次的樣本量較少，若閾值較小，可能會出現一種天氣狀況下的降雨被分為兩場降雨；②無雨持續時間大于或等于6 h[24]，若無雨持續時間小于6 h，則認為是一場降雨事件，反之則為兩場降雨事件。

根據該原則，阿湖水庫以上流域共挑選出83場降雨，流域平均每年約發生7.5場降水總量大于20.0 mm的降雨事件，每場降雨事件的平均降水總量達到44.9 mm，場次降水總量最高達144.4 mm，平均降雨歷時為16.6 h，最大降雨歷時為48 h，其中最大降雨歷時與最大降水總量相對應。降水總量和降雨歷時通過了95%的相關性顯著檢驗，表明該地區適合采用Copula函數建立降水總量與降雨歷時之間的聯合分布。

2.2 降水總量與降雨歷時的聯合分布

采用二維Copula函數構建降水總量與降雨歷時的聯合分布。不同邊緣分布下降水總量與降雨歷時的BIC和AIC指標值如表1所示，依據BIC和AIC最小原則優選出的邊緣分布分別為GPA分布和Weibull分布，且兩種邊緣分布函數均通過了K-S擬合優度檢驗。降水總量邊緣分布中的形狀、尺度和位置參數分別為k=0.43，σ=23.87，θ=20.0；降雨歷時邊緣分布中的尺度和形狀參數分別為A=18.45，B=1.55。經驗分布與理論分布的擬合結果如圖2(a)和(b)所示，可見經驗頻率與理論頻率非常接近，擬合效果較好。因此，對阿湖水庫地區降水總量和降雨歷時的模擬適宜采用GPA分布和Weibull分布。

表1 降水總量與降雨歷時不同邊緣分布下的BIC和AIC指標值

(a) 降水總量邊緣分布 (b) 降雨歷時邊緣分布 (c) 聯合概率Q-Q圖

分別采用MLE和MCMC[22]方法估計Copula函數的參數，表2為依據MLV最大、AIC和BIC最小準則優選出的Copula備選函數。MCMC算法下的參數估計優化了原本在MLE中表現不好的函數(Roch-Alegre、t和Gaussian)，使該函數更適用于表示變量間的聯合分布。雖然MLE法可以高效地估計參數，但對于有些Copula函數形式，MLE法容易陷入局部最優解，而MCMC算法下的參數尋優具有多個初始鏈，可以在全局尋優[22]和迭代次數足夠多的前提下，更易找到最優解。

在本研究區中，兩種參數估計算法下的最優聯合分布函數形式均是BB1，并且在兩種參數估計方法下的指標值類似。理論聯合概率與經驗聯合概率的RMSE(0.171 5)和NSE(0.995)均表現較好，結合降水總量和降雨歷時的聯合分布Q-Q圖(圖2(c))，發現圖中所有散點均落在45°對角線附近，由此可知，BB1形式可以較好地擬合降水總量和降雨歷時的聯合分布，最終其聯合分布和條件分布為

C(u,v)={1+[(u-θ1-1)θ2+

(v-θ1-1)θ2]1/θ2}-1/θ1

(3)

CT|R(v|u)=[1+(x+y)1/θ2]-1/θ1-1·

(x+y)1/θ2-1x1-1/θ2u-θ1-1

(4)

式中：C(u,v)表示降水總量和降雨歷時的聯合分布函數；u和v分別表示降水總量與降雨歷時的邊緣分布函數,u=FR(R)，v=FT(T)；CT|R(v|u)表示在降水總量為R發生的概率u下，降雨歷時為T發生的條件概率；x=(u-θ1-1)θ2；y=(v-θ1-1)θ2；參數θ1和θ2在MCMC算法優化下的值分別為0.220 3和1.33。

2.3 代表雨型分析

將83場降雨事件的降雨過程線進行無量綱化處理，根據流域平均降雨歷時(約17 h)，將降雨過程線分為17段，即K=17，并采用K均值聚類方法對其進行分類。我國雨型劃分采用的聚類數一般為3～6[25]，若聚類數太少，會忽略頻率低但對防洪較不利的雨型；若聚類數太多，則會導致聚類雨型間差異較小。本實例選擇聚類數為4和5，得到的雨型分類如圖3所示。阿湖水庫地區的雨型主要以單峰型為主，雨峰在前期且降雨主要集中在前中期的降雨事件頻率最高。聚類數為4時，R4-1、R4-2、R4-3、R4-4雨型的發生頻率依次為32.53%、28.92%、20.48%和18.07%。雨型差異主要體現在雨峰位置和降雨集中程度。R4-3的雨峰較R4-1更靠前，R4-1的降雨集中在雨峰兩側，而R4-3的降雨集中在前期，雨強在雨峰后陡然減小。R4-4的雨峰較R4-2更靠后，R4-2的降雨集中在雨峰之前，而R4-4的降雨集中在雨峰之后。相較于聚類數為4時的雨型劃分結果，聚類數為5時，R5-1、R5-2、R5-3、R5-4、R5-5雨型的發生頻率依次為33.73%、28.92%、20.78%、10.84%、6.02%。將雨峰在中(R5-2)、極端在前(R5-4)和極端在后(R5-5)細分出來。R5-2為雙峰型，兩個雨峰接近且降雨集中在中期，發生頻率較高為28.92%。雖然R5-5發生頻率僅有6.02%，但在實際防洪規劃設計及管理中該雨型情況較不利，需要單獨著重考慮這一雨型。因此，最終選擇雨型聚類數為5。

表2 兩種參數估計方法下的優選函數及指標值

(a)R4-1 (b)R4-2 (c)R4-3 (d)R4-4

(e)R5-1(f)R5-2(g)R5-3(h)R5-4(i)R5-5

本文將降水總量分為3級：20～50 mm、50～100 mm和100 mm以上，將降雨歷時分為4級：1～6 h、6～12 h、12～24 h和24 h以上，采用列聯法統計在降水總量和降雨歷時分級下各雨型的頻率(表3)。在降水總量和降雨歷時分級下，各雨型頻率和代表雨型均有較大的差別。對于20～50 mm的降雨事件，雨型多樣，隨著降雨歷時的增加，雨型從以R5-3和R5-1為主變換為以R5-1和R5-2為主，選擇最有可能的雨型(頻率最高)為代表雨型。對于50～100 mm的降雨事件，12 h以內R5-3發生頻率較高，為代表雨型；在12～24 h內R5-1和R5-2居多；當降雨歷時大于24 h時，R5-1和R5-5居多，均選擇R5-1為代表雨型。當降水總量大于100 mm時(樣本量較少)，降雨歷時以長歷時為主，當降雨歷時在12～24 h和24 h以上時，分別選擇R5-2和R5-3為代表雨型；當歷時較短時，選擇對防洪較為不利的R5-5為代表雨型。

表3 5種雨型在降水總量和歷時分級下的頻率及代表雨型

2.4 同一降水總量下不同歷時的降雨情景

以降水總量為40 mm、70 mm和110 mm為例，實際降雨情景和采用典型場次同倍比放大得到的降雨情景分別見圖4和圖5。依據式(4)，挑選并計算條件概率P(T≤t*|R=r*)分別為0.1、0.5和0.9時生成的降雨情景結果見圖6。3種條件概率分別代表在同一降水總量下，降雨歷時為短、中和長的降雨情景，也可反映降雨集中、均勻和分散時的降雨情景。在實際降雨情景中，降雨歷時隨著降水總量的增加而增加。本文生成的降雨情景，在同一條件概率下，亦能夠模擬出該特性。而采用典型場次同倍比放大得到的降雨情景，在不同降水總量下的降雨歷時均為21 h(典型降雨歷時為21 h)，在降水總量較小時(如40 mm和70 mm)，坦化了降雨過程，既不能反映實際降雨情況，也沒有考慮對防洪的不利影響。

利用本文提出的方法生成的降雨情景，總是可以找到與實際降雨情景較為相似的降雨特性或者生成對防洪較不利的降雨情景。當降水總量為40 mm和70 mm時，條件概率為0.1時的降雨情景(圖6(a)和(d))能夠反映實際降雨情景的集中程度。由于僅考慮了在降水總量和歷時分級下最有可能發生的雨型，在雨峰強度和雨峰位置上難免與實際降雨情景有出入。當降水總量為110 mm時，實際降雨情景中降雨歷時為31 h，降雨主要集中在第20～30 h內。在本文方法生成的降雨情景中，雖然對于歷時相似的降雨情景(圖6(i))，雨峰與實際雨峰有較大差異，但條件概率為0.1的降雨情景(圖6(g))歷時短，降雨集中，其雨峰強度與實際降雨情景的雨峰強度相近，且該降雨情景對防洪十分不利。

因此，本文方法生成的降雨情景考慮了降水總量、降雨歷時和降雨時程分配的相關性，可以較好地反映多種實際降雨情景或生成對防洪不利的降雨情景，并可定量計算同一降水總量下不同降雨情景的概率。

2.5 抗暴雨能力預報結果對比

以2020年7月12日8時阿湖水庫的狀況為例，該時刻庫水位為25.54 m。將傳統典型降雨同倍比縮放得到的降雨情景和本文生成的降雨情景分別輸入到水庫抗暴雨能力模型[14]中，傳統方法預報此時阿湖水庫的抗暴雨能力為160 mm，而本文方法預報的水庫抗暴雨能力分布在146～169 mm區間內(圖7)，相應概率在9.9%～93.1%之間，其中160 mm的概率為62.2%。

實測結果表明，7月12日8時至7月15日8時實際降水總量為64.4 mm，水庫水位于7月12日21時漲至最高(26.68 m)，低于設計水位(27.73 m)1.05 m，實際發生降雨小于抗暴雨能力，水庫仍有剩余庫容。根據實測降雨和入庫徑流計算此次降雨產流系數為0.58，若將7月12日8時剩余的防洪庫容(1 677萬 m3)折算成水庫抗暴雨能力為149.8 mm，該值為接近實際水庫抗暴雨能力的一種情況，落在本文的預報區間內，發生概率為11.6%。從以上分析來看，利用本文提出的降雨情景構建方法，可以實現阿湖水庫抗暴雨能力的概率預報，結果合理。

(a)降水總量40 mm、降雨歷時5 h (b)降水總量70 mm、降雨歷時13 h (c)降水總量110 mm、降雨歷時31 h

(a)降水總量40 mm (b)降水總量70 mm (c)降水總量110 mm

(a)P(t≤6 h|R=40 mm)=0.1 (b)P(t≤16 h|R=40 mm)=0.5 (c)P(t≤29 h|R=40 mm)=0.9

(d)P(t≤9 h|R=70 mm)=0.1 (e)P(t≤22 h|R=70 mm)=0.5 (f)P(t≤37 h|R=70 mm)=0.9

(g)P(t≤12 h|R=110 mm)=0.1 (h)P(t≤29 h|R=110 mm)=0.5 (i)P(t≤46 h|R=110 mm)=0.9

圖7 阿湖水庫2020年7月12日8時水庫抗暴雨能力概率預報結果

3 結 論

a． 提出了基于Copula函數構建降水總量和降雨歷時的聯合分布，結合K均值聚類法得到的代表雨型，生成同一降水總量下不同概率降雨情景的方法。相比于傳統典型暴雨放大方法的輸入，本文將概率降雨情景應用到水庫抗暴雨能力計算中，考慮了流域不同情況降雨的可能性及其對產匯流的影響，為實現中小水庫抗暴雨能力從定量預報到概率預報，拓展中小水庫抗暴雨能力預報方法提供了技術支撐。

b． 基于多樣的Copula備選函數形式并結合MCMC參數估計方法，可較好地找到描述變量聯合關系的最優Copula函數。在阿湖水庫流域，BB1函數形式可較好地反映降水總量與降雨歷時的聯合分布。

c． 雨型聚類數選取需根據流域調整。通過對比不同聚類數結果下的降雨特性(如雨峰位置和降雨集中程度)和對實際防洪規劃中不利雨型的考慮，分析選擇合適的聚類數。在阿湖水庫地區，雨型以單峰型為主，雨型聚類數為5較合適。