混凝土壩洪水漫頂生命風(fēng)險評價

李宗樾，余 雷，胡澤林，程 井

(1.貴州省大壩安全監(jiān)測中心，貴州 貴陽 550002；2.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098)

截至2019 年底，中國已建成各類水庫超過98 000 座，是目前世界修建水庫最多的國家[1-3]。水庫在發(fā)揮巨大社會效益和經(jīng)濟(jì)效益的同時，也存在不可規(guī)避的風(fēng)險[4]。為了盡可能減少水庫大壩失事帶來的損失，大壩風(fēng)險管理顯得格外重要。我國自從1954 年有較為完整的潰壩紀(jì)錄以來，隨著水庫的運(yùn)行管理等各項(xiàng)工作的不斷完善，大中型水庫的潰壩事件隨時段劇減。1954—2018 年水庫共潰壩3 541 座[5]，其中1954—1982 年共潰壩3 115 座，年均潰壩107.4 座；1983—1999 年共潰壩332 座，年均潰壩19.5 座；2000—2018 年共潰壩84 座，年均潰壩4.4 座，年均潰壩概率降至0.45×10-4[4-8]。

導(dǎo)致潰壩的原因主要包括洪水漫頂、大壩質(zhì)量問題及管理不當(dāng)?shù)?，其中洪水漫頂?0%以上[6]。特別是近年來，全球氣候變化導(dǎo)致的極端天氣頻發(fā)，超標(biāo)洪水頻率明顯上升，給水庫大壩帶來了嚴(yán)峻挑戰(zhàn)[9]，如2020 年烏茲別克斯坦的薩爾多巴水庫潰壩，以及美國密歇根州的伊登維爾大壩及桑福德大壩潰壩[5]，這些潰壩事故造成嚴(yán)重的經(jīng)濟(jì)和生命損失。如何合理評估庫壩的各類風(fēng)險，成為大壩運(yùn)行維護(hù)的重要內(nèi)容[4]。

優(yōu)質(zhì)實(shí)體混凝土壩可以長期服役，但在遭遇超標(biāo)洪水且泄流能力不足時，仍存在較大的洪水漫頂及潰壩風(fēng)險[10-11]。壩前最高水位與入庫洪水息息相關(guān)，而傳統(tǒng)的設(shè)計(jì)洪水頻率分析通常為單一變量分析，無法全面反映洪峰與洪量之間的相關(guān)關(guān)系?；贑opula 函數(shù)的多變量設(shè)計(jì)洪水頻率分析能更全面地描述洪水的內(nèi)在規(guī)律[12-16]。將多變量洪水分析引入洪水漫頂生命風(fēng)險評價，更貼合洪水發(fā)生及漫頂?shù)膶?shí)際情況。本文以貴州某水電站為例，研究混凝土壩遭遇特大洪水并發(fā)生漫頂情況下的生命風(fēng)險評價方法，基于Copula 函數(shù)構(gòu)建洪峰與洪量關(guān)系，并采用MC 抽樣方法通過三點(diǎn)式變倍比放大法得到年最大入庫洪水過程線概率序列；依據(jù)隨機(jī)生成的洪水過程線模擬洪水淹進(jìn)過程，進(jìn)行生命風(fēng)險評估，并判定風(fēng)險的可接受程度。

1 混凝土大壩安全評價的過程及原理

1.1 擬定年最大入庫洪水過程線概率序列

三點(diǎn)式變倍比放大法，是一種同時考慮峰量、利用典型洪水過程線推求設(shè)計(jì)洪水過程線的放大方法，對多變量洪水分析中用到的變倍比放大法[13-15]的前后收斂性問題進(jìn)行了改善，見式（1）～（4）。與以往洪水放大法采用單變量進(jìn)行分析不同，該方法基于Copula 函數(shù)進(jìn)行多變量洪水分析，可同時考慮洪峰、洪量及典型洪水過程線形狀對設(shè)計(jì)洪水過程線的影響。

式中：Qmax和QDmax為模擬洪水過程線和典型洪水過程線的最大洪峰流量；Q和QD為模擬洪水過程線和典型洪水過程線某一時段的平均流量；Qt和QD(t)為t時刻模擬洪水過程線和典型洪水過程線對應(yīng)流量；QD(0)和QD(T))為典型洪水過程線的初始流量和結(jié)束流量；tQDmax為典型洪水過程線最大洪峰流量相應(yīng)時刻。

Copula 函數(shù)可以連接多個隨機(jī)變量的邊緣分布來構(gòu)造聯(lián)合分布。近年來在水文領(lǐng)域經(jīng)常使用的3 種單參數(shù)二維Archimedean Copula 函數(shù)[16]見表1?；贕umbel-Hougaard 函數(shù)構(gòu)建的洪峰與洪量關(guān)系，通過MC（Monte Carlo method）法隨機(jī)生成較大組數(shù)的洪峰與洪量。利用生成的洪峰、洪量及典型洪水過程線通過三點(diǎn)式變倍比放大法得到年最大入庫洪水過程線。就水庫大壩工程防洪安全而言，壩前最高水位起著決定性作用，將壩前最高水位的重現(xiàn)期看作設(shè)計(jì)洪水過程線的重現(xiàn)期。用年最大入庫洪水過程線序列進(jìn)行調(diào)洪演算得到壩前最高水位序列，進(jìn)一步可得到壩前最高水位的累積頻率。將壩前最高水位與年最大入庫洪水過程線相對應(yīng)，便可得到年最大入庫洪水過程線的累積頻率，即年最大入庫洪水過程線發(fā)生概率。

表1 Archimedean Copula 函數(shù)參數(shù)θ 與Kendall 相關(guān)系數(shù)τ 的關(guān)系Tab.1 Relationship between Archimedean Copula function parameter and Kendall rank correlation coefficient

1.2 大壩洪水漫頂生命損失評估

文獻(xiàn)[17]通過對8 座水庫潰壩洪水嚴(yán)重性的分析，并結(jié)合我國歷史潰壩生命損失的主要規(guī)律，提出適合我國水庫潰壩實(shí)際情況的生命損失模型。模型主要考慮3 個部分：風(fēng)險人口PAR、風(fēng)險人口死亡率f及潰壩對下游造成的嚴(yán)重性系數(shù)α，LOL=RARfα，RAR=mPSP，α=m1+bm2，其中：mP和SP分別為下游受災(zāi)地區(qū)的人口密度和淹沒面積；m1和m2分別為生命損失的直接和間接影響因素的嚴(yán)重性影響因子，文獻(xiàn)[17]定性給出了影響因素的嚴(yán)重性程度判別建議標(biāo)準(zhǔn)，m1/m2≤1；b為生命損失間接影響因素的嚴(yán)重性程度影響因子系數(shù)，0＜b＜1，采用0.25；f為風(fēng)險人口死亡率，主要考慮潰壩洪水嚴(yán)重性、警報(bào)時間和風(fēng)險人口對潰壩洪水嚴(yán)重性理解程度3 個影響因素，文獻(xiàn)[17]給出了適合我國水庫潰壩的風(fēng)險人口死亡率建議表。

利用1.1 節(jié)得到的年最大入庫洪水過程線序列作為HEC-RAS 軟件的基礎(chǔ)數(shù)據(jù)，進(jìn)行大壩洪水漫頂模擬演進(jìn)，并對洪水漫頂后對下游造成的生命損失進(jìn)行評估。

1.3 生命風(fēng)險標(biāo)準(zhǔn)制定

ALARP (As Low As Reasonably Practicable)是目前國際公認(rèn)的風(fēng)險準(zhǔn)則[6,18]，依據(jù)可接受風(fēng)險水平與可容忍風(fēng)險水平將風(fēng)險劃分為可接受風(fēng)險、可容忍風(fēng)險和不可容忍風(fēng)險3 個區(qū)域，這兩個風(fēng)險水平即為所謂的風(fēng)險標(biāo)準(zhǔn)。F-N線法是確定風(fēng)險標(biāo)準(zhǔn)的主要方法之一，N為死亡人數(shù)，F(xiàn)為N的累計(jì)頻率限制曲線。式（5）表示F的超過概率與死亡人數(shù)N之間關(guān)系，通過在雙對數(shù)坐標(biāo)系中確定起點(diǎn)位置、斜率和極值線來描述F-N標(biāo)準(zhǔn)線。

式中：FN(x)為年死亡人數(shù)小于x的概率分布函數(shù)；C為常數(shù)；n為標(biāo)準(zhǔn)線的斜率。

生命風(fēng)險標(biāo)準(zhǔn)的構(gòu)建與社會、環(huán)境及當(dāng)?shù)氐恼巍⒔?jīng)濟(jì)、文化等因素有關(guān)，不同國家地區(qū)考慮的側(cè)重點(diǎn)有所不同。2006 年李雷等[6]遵循ALARP 準(zhǔn)則，依據(jù)我國年均潰壩率及潰壩人數(shù)的估算，提出我國生命風(fēng)險標(biāo)準(zhǔn)的F-N線建議圖。2015 年李宗坤等[18]結(jié)合我國水庫大壩安全現(xiàn)狀及其他行業(yè)的風(fēng)險標(biāo)準(zhǔn)情況，對F-N標(biāo)準(zhǔn)線的參數(shù)取值進(jìn)行了更新：建議對水庫大壩可容忍風(fēng)險標(biāo)準(zhǔn)的C值取10-2，可接受風(fēng)險標(biāo)準(zhǔn)的C值取10-3；將大中型水庫歸屬于風(fēng)險厭惡型（n值取2），小型水庫歸屬于風(fēng)險中立型（n值取1）；同時考慮了中國大壩事故概率和中國可靠度標(biāo)準(zhǔn)分別進(jìn)行極值線的計(jì)算，并取較小值作為標(biāo)準(zhǔn)，見圖1。

圖1 大中型及小型水庫社會生命風(fēng)險標(biāo)準(zhǔn)Fig.1 Social life risk criteria for large and medium-sized reservoirs and small reservoirs

2 工程應(yīng)用

某水電站水庫位于貴州省黎平縣境內(nèi)，工程等別為Ⅲ等，主要建筑物為3 級。水電站原大壩為拋物線雙曲拱壩，壩頂高程288.8 m，最大壩高65.8 m，壩頂設(shè)置0.9 m 高防浪墻。2016 年因調(diào)度不當(dāng)發(fā)生漫壩事故，壩前水位達(dá)到歷史最高洪水位290.1 m，高于防浪墻頂0.4 m，造成較大經(jīng)濟(jì)損失和社會影響。2019 年針對工程存在的抗洪能力不足問題進(jìn)行了除險加固，拆除原有溢流堰頂胸墻，對大壩壩頂上游面防浪墻進(jìn)行加高處理，增強(qiáng)大壩的抗洪能力。

2.1 擬定洪水過程線序列

根據(jù)1964—2015 年的實(shí)測水文數(shù)據(jù)，結(jié)合調(diào)查所得歷史洪水，推得該水電站的設(shè)計(jì)洪峰及洪量。采用P-Ⅲ分布曲線分別描述洪峰流量邊緣分布u和最大3 日洪量邊緣分布v。構(gòu)建P-Ⅲ分布的相關(guān)系數(shù)結(jié)果為：最大洪峰流量q為693 m3/s，變差系數(shù)Cv為0.68，偏差系數(shù)與變差系數(shù)的比值Cs/Cv為4.5；最大3 日洪量w為0.34 億m3，Cv為0.65，Cs/Cv為4.0。

利用表1 的Gumbel-Hougaard Copula 函數(shù)構(gòu)建洪峰流量與最大3 日洪量二者邊緣分布的聯(lián)合分布C(u,v)，其中τ=0.92?；贛C 法在（0,1）區(qū)間隨機(jī)生成10 萬組相互獨(dú)立且均勻分布的u與條件概率Pv|u，利用式（6）得到10 萬組符合C(u,v)分布的(u,v)見圖2。利用反函數(shù)u-1、v-1推得10 萬組(q,w)，隨機(jī)生成的q與w的P-III 分布見圖3～4。利用三點(diǎn)式變倍比放大法式（1）～（4）對典型洪水過程線（圖5）進(jìn)行放大模擬得到10 萬條年最大入庫洪水過程線?；谒心曜畲笕霂旌樗^程線對除險加固后的水庫進(jìn)行調(diào)洪演算，得到壩前最高水位累積頻率曲線如圖6，相應(yīng)可得年最大入庫洪水過程線發(fā)生概率。

圖2 (u,v)抽樣分布Fig.2 Sampling of (u,v)

圖3 隨機(jī)生成的洪峰q 的P-III 分布Fig.3 P-III distribution of randomly generated flood peak q

圖4 隨機(jī)生成的洪量w 的P-III 分布Fig.4 P-III distribution of the randomly generated flood volume w

圖5 典型洪水過程線Fig.5 Typical flood process line

圖6 壩前最高水位累積頻率曲線Fig.6 Cumulative frequency curve of highest water levels in front of the dam

式中：Pv|u為條件概率；C(u,v)為雙變量Copula 函數(shù)。

2.2 大壩洪水漫頂生命風(fēng)險評價

工程所處位置偏僻，受工程規(guī)模及勘測精度等客觀條件限制，對下游3 個村鎮(zhèn)的高程取平均海拔高程，人口采用均勻分布，具體信息見表2。

表2 各村鎮(zhèn)信息Tab.2 Information of villages and towns

利用HEC-RAS 軟件對上游水庫、拱壩及下游可能的淹沒區(qū)域進(jìn)行模擬。上游水庫初始水位設(shè)置為280 m（溢流堰堰頂高程），下游區(qū)域以河床起算，模擬時間設(shè)置為2050 年6 月10 日8:00—2050 年6 月13 日11:00，入庫洪水的模擬在隨機(jī)生成的年最大入庫洪水過程線序列中選取1 條輸入。

此次模擬的年最大入庫洪水過程序列有10 萬條，無法逐條進(jìn)行模擬，但由大壩洪水漫頂造成下游海拔相對較高地區(qū)淹沒的情況為極值事件，只需對幾個極小概率的洪水過程線的情況進(jìn)行洪水模擬演進(jìn)。圖7 為HEC-RAS 輸入概率為1.0×10-5的年最大入庫洪水過程線時（本次模擬中發(fā)生概率最小的年入庫洪水過程線）下游最大淹沒水深示意。從圖7可見，下游村落中只有最靠近大壩的SJ 村出現(xiàn)被洪水淹沒的情況。模擬的目的是為了進(jìn)行生命損失的評估，則此次只需要對能淹沒SJ 村的年入庫洪水過程線進(jìn)行模擬。

圖7 概率10-5 的洪水漫頂工況的下游最大淹沒水深示意（2050-06-10T15:43）（單位:m）Fig.7 Schematic diagram of the downstream maximum submerged depth under flood overtopping condition with probability of 10-5 (2050-06-10T15:43) (unit:m)

在SJ 村區(qū)域內(nèi)，選取具有代表性且盡可能均勻分布的5 個測點(diǎn)，將人口均分在各測點(diǎn)上，從出現(xiàn)概率為1.0×10-5的年最大入庫洪水過程線開始，概率從小到大分別輸入HEC-RAS 軟件中進(jìn)行洪水模擬演進(jìn)。對HEC-RAS 洪水模擬演進(jìn)的結(jié)果，查看各測點(diǎn)洪水淹沒水深和淹沒洪水流速（圖8～9），最終用洪水嚴(yán)重性即最大淹沒水深與流速的乘積值來反映SJ 村的整體受災(zāi)嚴(yán)重性，結(jié)果見表3。當(dāng)模擬輸入的年最大入庫洪水過程線概率大于等于7.0×10-5時，SJ 村不再出現(xiàn)洪水淹沒，則更大概率的洪水過程線無需進(jìn)行模擬。

圖8 SJ 村淹沒水深過程線Fig.8 Process line of submerged water depth in SJ village

根據(jù)文獻(xiàn)[17]建立的生命損失模型，將人口集中于5 個測點(diǎn)上，各測點(diǎn)的人數(shù)均分為230 人。SJ 村風(fēng)險人口死亡率f的取值，由SJ 村村民對潰壩洪水嚴(yán)重性理解清楚（SJ 村地處多庫壩的省份）、漫頂洪水屬于完全報(bào)警的情況（HEC-RAS 洪水模擬演進(jìn)結(jié)果顯示，從洪水開始漫頂?shù)絊J 村出現(xiàn)淹沒的時間間隔為2.0～2.5 h）、SJ 村的受災(zāi)嚴(yán)重性（表3）3 個方面，根據(jù)文獻(xiàn)[17]的風(fēng)險人口死亡率建議表進(jìn)行取值。最終得到各極小概率洪水的生命損失及死亡率，結(jié)果見表3。

圖9 SJ 村淹沒洪水流速過程線Fig.9 Process line of flood velocity in SJ village

由于模擬的最小概率洪水過程線對應(yīng)的壩前最高水位已經(jīng)達(dá)到設(shè)置水位條件的最大邊界，對于超出邊界水位，例如超出水位-庫容曲線范疇，無法探究規(guī)律，一般認(rèn)為發(fā)生概率為零。由表3 可以看出，經(jīng)除險加固后的拱壩的洪水漫頂生命風(fēng)險超過1 人的概率為2.5×10-9～8.0×10-9。該水電站為中型水庫，對比圖1 大中型水庫社會生命風(fēng)險標(biāo)準(zhǔn)，可見除險加固后的水電站洪水漫頂生命風(fēng)險處于可接受范圍。

表3 SJ 村受災(zāi)嚴(yán)重性Tab.3 Severity of the disaster in SJ village

3 結(jié)語

本文提出了一種洪水概率序列生成-漫頂洪水演進(jìn)-下游生命損失評估的混凝土壩漫頂生命風(fēng)險定量評估方法。該方法基于Copula 函數(shù)與三點(diǎn)式變倍比放大法進(jìn)行設(shè)計(jì)洪水多變量聯(lián)合分析，可考慮洪水特征量間的相互關(guān)系，得到的年最大入庫洪水過程線發(fā)生概率序列更加符合洪水內(nèi)在規(guī)律；漫頂生命風(fēng)險定量評估與符合我國現(xiàn)狀的風(fēng)險標(biāo)準(zhǔn)共同構(gòu)成了完整的混凝土壩生命風(fēng)險評價體系。以貴州省某拱壩為例進(jìn)行了分析評估，結(jié)果表明，大壩加固后洪水漫頂對下游村鎮(zhèn)造成生命損失的風(fēng)險較小，屬于社會可接受的范圍。