基于運行風(fēng)險度與馬爾科夫鏈的大壩運行壽命評估

范振東

（國家能源局大壩安全監(jiān)察中心，浙江杭州，311122）

1 概述

大壩運行壽命評估是一個多學(xué)科交叉的復(fù)雜技術(shù)問題，同時受大壩龐大的結(jié)構(gòu)體系、復(fù)雜的地質(zhì)條件和運行環(huán)境、材料老化等因素影響，很難對大壩運行壽命進(jìn)行量化評估[1-2]。然而，對大壩運行壽命進(jìn)行評估和預(yù)測可以為制定大壩除險加固決策提供參考依據(jù)，從而在全壽命周期內(nèi)基于保障大壩安全運行的前提下控制除險加固成本。為此，國內(nèi)許多專家學(xué)者對如何評估大壩運行壽命做了深入研究。如蘇懷智[3-4]等提出了基于時變可靠度的大壩使用壽命評估方法，但該方法需要事先假定壩體材料參數(shù)時變特性，具有較強(qiáng)的人為主觀性。朱伯芳[5]從混凝土碳化、凍融破壞等方面考慮大壩使用壽命，然而由于大壩體積龐大、碳化速率極為緩慢，由此計算出來的使用壽命很長，與實際不符。宋恩來[6]從大壩原型監(jiān)測數(shù)據(jù)出發(fā)，反演壩體混凝土綜合彈性模量的時變規(guī)律，但該方法沒有建立綜合彈模和運行壽命的對應(yīng)關(guān)系，且計算工作量大。

眾所周知，大壩的運行性態(tài)可以通過原型監(jiān)測資料反映，因此，大壩的原型監(jiān)測資料也隱含了大壩運行壽命的客觀信息。筆者從大壩監(jiān)測系統(tǒng)實測資料出發(fā)，基于大壩運行風(fēng)險度，運用馬爾科夫鏈方法對大壩運行壽命評估進(jìn)行初步探索，為預(yù)測大壩運行壽命探索一種新思路。

2 大壩運行風(fēng)險度評估體系和評判準(zhǔn)則

2.1 風(fēng)險度定義

大壩的運行風(fēng)險率RI也可稱為大壩失事的總概率Pf，其定義可用下列數(shù)學(xué)表達(dá)式表示：

式中：S為相互獨立隨機(jī)變量組成的大壩廣義抗力；F為相互獨立隨機(jī)變量組成的大壩廣義荷載；fSF(S,F)為荷載F和抗力S的聯(lián)合概率密度函數(shù)。

但多數(shù)情況下，受資料缺乏、求解難度大等限制，尚無法根據(jù)式（1）直接求解大壩運行風(fēng)險率RI，因此吳中如[7]提出了根據(jù)各種風(fēng)險因素實測資料來確定大壩運行風(fēng)險率的方法，并引入了“風(fēng)險度”的概念來代替風(fēng)險率，用于表征測點實測成果反映的相應(yīng)點位處物理量潛在的風(fēng)險程度，考慮已用的符號，仍用符號RI來表示風(fēng)險度。

2.2 大壩運行風(fēng)險度評估體系

大壩運行風(fēng)險度評估[7]的總目標(biāo)為估計大壩運行風(fēng)險度（A層），按風(fēng)險部位劃分，可將準(zhǔn)則層（B層）分為壩體風(fēng)險度、壩基風(fēng)險度和近壩區(qū)風(fēng)險度，各準(zhǔn)則層對應(yīng)的指標(biāo)則為影響各準(zhǔn)則層風(fēng)險度的風(fēng)險因素（C層），監(jiān)測項目的實測資料可以反映這些風(fēng)險因素的特性和變化規(guī)律，例如，變形測值反映壩體、壩基和近壩區(qū)邊坡的穩(wěn)定，相同荷載條件下，應(yīng)力應(yīng)變測值變化反映壩體材料力學(xué)特性的變化等。大壩運行風(fēng)險度評估體系如圖1所示。

圖1 大壩運行風(fēng)險度評估體系Fig.1 Dam operation risk assessment system

如圖1所示，影響大壩安全度的風(fēng)險因素眾多，而這些因素有些是主要因素，有些是次要因素，如果對風(fēng)險的嚴(yán)重程度進(jìn)行初步判斷，將次要因素從風(fēng)險因素集中剔除，能大大減少大壩運行風(fēng)險度分析計算工作。

2.3 大壩運行風(fēng)險度評判準(zhǔn)則

水利部2003年頒布的《水庫大壩安全鑒定辦法》將大壩分為一類壩、二類壩和三類壩，其中一類壩安全可靠，大壩工作狀態(tài)正常，在設(shè)計工況下能正常運行；二類壩類似于病壩，需要在一定控制條件下運行；三類壩類似于險壩，工程存在較嚴(yán)重安全隱患，不能按設(shè)計正常運行。參考文獻(xiàn)[7]將大壩運行風(fēng)險度分為四個等級，相應(yīng)的評語依次為N1（微風(fēng)險）、N2（低風(fēng)險）、N3（中風(fēng)險）和N4（高風(fēng)險），其中取一類壩中的20%劃歸為“微風(fēng)險”，其余的劃歸為“低風(fēng)險”；將二類壩全部劃歸為“中風(fēng)險”；將三類壩全部劃歸為“高風(fēng)險”。

各類事件的先驗概率P(Ni)可以通過我國水利系統(tǒng)現(xiàn)有的大壩安全分類統(tǒng)計情況進(jìn)行確定，表1統(tǒng)計的為某一年水利系統(tǒng)80 427座大壩安全分類情況[7]。

表1 某年全國水利系統(tǒng)水庫大壩安全分類Table 1 Safety classification of reservoir dams in China's water conservancy system in a given year

通過計算P(N1)=0.06，P(N2)=0.25，P(N3)=0.33，P(N4)=0.36。取大壩運行風(fēng)險度值區(qū)間為[0,1]，則根據(jù)各類事件的先驗概率，可以確定各風(fēng)險等級對應(yīng)的大壩運行風(fēng)險度區(qū)間如下：

式（2）表示若某一風(fēng)險因素的風(fēng)險度值在[0,0.06]區(qū)間時，該因素的特性為“微風(fēng)險”，以此類推。當(dāng)計算出風(fēng)險因素的風(fēng)險度值，即可通過評價體系自下而上逐層確定大壩運行的綜合風(fēng)險度。

3 基于實測資料的大壩運行風(fēng)險度分析

大壩是一個由不同部位組合而成的復(fù)雜體系，其綜合運行風(fēng)險度實質(zhì)上是多個組成部位風(fēng)險度的綜合反映。各個部位的風(fēng)險度可通過不同項目各個測點的實測資料分析得到，因此大壩綜合運行風(fēng)險度確定的過程可分解為大壩運行單項、局部和綜合風(fēng)險度分析。

文獻(xiàn)[7]通過數(shù)值趨勢和時效趨勢兩方面建立了單個測點的單項運行風(fēng)險度，具體計算方法如下。

圖2為單測點單項測值風(fēng)險度的評價級別劃分方案，圖中，y為效應(yīng)量的實測值，如變形、滲流等；y為根據(jù)監(jiān)控模型求得的效應(yīng)量的擬合值，y=f(x)，x指影響效應(yīng)量的自變量因子；S為監(jiān)控模型的均方差；ymax、ymin為擬定的安全閾值上限值和下限值，但對于有些效應(yīng)量（如滲流量等），沒有下限值。

圖2中，把大壩效應(yīng)量測值y分為4個區(qū)域，如下式：

時效分量可以充分反映大壩結(jié)構(gòu)風(fēng)險的變化趨勢，圖3為單測點測值風(fēng)險度趨勢表現(xiàn)的評價級別劃分方案，圖中將測值的時效分量劃分為4種形式：

圖2 單測點測值風(fēng)險度數(shù)值表現(xiàn)評價級別劃分Fig.2 Single point operation risk division by numerical performance

圖3 單測點測值風(fēng)險度趨勢表現(xiàn)評價級別劃分Fig.3 Single point operation risk division by trend performance

（2）時效分量初期增長較快，后期趨于穩(wěn)定，如曲線B所示，。這種情況最為常見，表示大壩的運行情況良好，與之對應(yīng)的風(fēng)險度評估等級為“低風(fēng)險N2”。

式（3）中，A、B、C、D區(qū)分別對應(yīng)該測點測值風(fēng)險P1特性為微風(fēng)險、低風(fēng)險、中風(fēng)險和高風(fēng)險，圖3中，A、B、C、D四條趨勢線分別對應(yīng)該測點趨勢風(fēng)險P2特性為微風(fēng)險、低風(fēng)險、中風(fēng)險和高風(fēng)險。根據(jù)不同項目，不同測點的實測數(shù)據(jù)和上述劃分準(zhǔn)則，可以建立一個模糊的單因素評價向量R：

式中：rij為該風(fēng)險元素Pi對評價向量Nj的模糊隸屬度，其中，r1j等于落在各個區(qū)間內(nèi)的測值個數(shù)占總測值個數(shù)的比。當(dāng)測點的時效分量與相應(yīng)于圖3中某一條曲線Nj變化一致時，r2j=1，而取r2k=0（k≠j）；當(dāng)測點的時效分量介于相應(yīng)于Nj和Nk的兩種曲線之間時，取r2j＞0，r2k＞0且r2j+r2k=1，其余r2m=0（m≠j,k）。

為了獲得數(shù)值表現(xiàn)和趨勢表現(xiàn)的綜合風(fēng)險度，需要構(gòu)造兩者的權(quán)重向量W=[ω1ω2]，權(quán)重向量的取值根據(jù)專家經(jīng)驗確定，最后得到測值綜合評價指標(biāo)向量Z：

根據(jù)模糊隸屬向量Z，可根據(jù)式（6）計算該測值的風(fēng)險度：

式中：RIE為某風(fēng)險因素的單項風(fēng)險度；zi為該風(fēng)險因素模糊隸屬向量Z中的第i個級別對應(yīng)的模糊隸屬度；ai、bi分別為zi對應(yīng)的風(fēng)險等級所屬區(qū)間的下限和上限。

在確定了大壩的全部單項風(fēng)險度值之后，可以根據(jù)文獻(xiàn)[7]給出的公式計算大壩局部運行風(fēng)險度RIB和大壩綜合運行風(fēng)險度RI。

4 基于Markov鏈的大壩時變風(fēng)險度模型構(gòu)建

4.1 Markov基本理論

若當(dāng)前隨機(jī)變量取值是通過它的前一個變量決定，而與其之后的變量無關(guān)，則該隨機(jī)過程稱為馬爾科夫過程。隨著馬爾科夫理論的完善與推廣，目前已被廣泛用于城市人口、水質(zhì)分析、語音和人臉識別等眾多領(lǐng)域的預(yù)測研究中[8-10]。

一個時間與空間都是離散的馬爾科夫過程叫做馬爾科夫鏈[11]。離散的馬爾科夫鏈用數(shù)學(xué)語言可表達(dá)如下：

設(shè)一隨機(jī)過程{Xt,t≥ 0}，如果Xt∈S(t=1,2,···)，其中S為一有限集合的狀態(tài)空間，且對于任意的i，j，i0，i1，···，in+1，S都有：

式中：i，j，i0，i1，…，in+1為系統(tǒng)所處的狀態(tài)；P為轉(zhuǎn)移概率矩陣；t為時間；X為離散的馬爾科夫鏈。

馬爾科夫鏈最重要且最基本變量是轉(zhuǎn)移概率矩陣P，其表達(dá)式如下：

式中：Pij(t)表示當(dāng)系統(tǒng)t時刻處于i狀態(tài)時，系統(tǒng)t+1時刻處于j狀態(tài)的概率，；m為系統(tǒng)所具有的狀態(tài)個數(shù)。

齊次馬爾科夫鏈的狀態(tài)矩陣與時間無關(guān)，系統(tǒng)任一時刻的所處狀態(tài)均可根據(jù)式（9）求得：

式中：Xt為系統(tǒng)在時間t所處的狀態(tài)；X0為系統(tǒng)的初始狀態(tài)。

4.2 轉(zhuǎn)移概率矩陣的確定

轉(zhuǎn)移概率矩陣是馬爾科夫鏈模型最核心的部分，目前常用的確定轉(zhuǎn)移概率矩陣的方法有經(jīng)驗判斷法、統(tǒng)計分析法、回歸分析法和基于逆陣的方法[12]。除了經(jīng)驗判斷法之外，剩余三種方法對監(jiān)測數(shù)據(jù)的要求均較高，計算工作量也較大，因此，為了簡化計算工作，通常需要做出某些合理的假定條件，從而簡化轉(zhuǎn)移概率矩陣。

結(jié)合工程結(jié)構(gòu)實際情況，做如下假定：（1）在正常運行條件下，不考慮除險加固措施影響，結(jié)構(gòu)在壽命周期內(nèi)的轉(zhuǎn)移概率不變；（2）結(jié)構(gòu)任意時刻所處的風(fēng)險等級只有兩種情況，維持當(dāng)前風(fēng)險等級狀態(tài)和向更高風(fēng)險等級狀態(tài)轉(zhuǎn)移。假定（1）認(rèn)為結(jié)構(gòu)的運行風(fēng)險度在壽命期內(nèi)連續(xù)均勻變化，這一假定不考慮轉(zhuǎn)移概率矩陣與時間的關(guān)系，故建立的馬爾科夫鏈為齊次的，可以大大減少計算工作量。假定（2）認(rèn)為結(jié)構(gòu)只能朝更高風(fēng)險等級發(fā)展，由于不考慮除險加固的影響，大壩的性能隨時間將逐漸劣化?；谏鲜黾俣l件的馬爾科夫鏈如圖4所示。

圖4 大壩運行風(fēng)險的馬爾科夫鏈模型Fig.4 Dam operation risk model based on Markov chain

由此可得轉(zhuǎn)移概率矩陣為

式中：Px為大壩結(jié)構(gòu)運行風(fēng)險從當(dāng)前狀態(tài)向下一狀態(tài)轉(zhuǎn)移的概率。

由式（10）可以看出，簡化后的概率矩陣只含有一個未知數(shù)Px，如果已知t1、t2（t2＞t1）時刻大壩運行風(fēng)險狀態(tài)的分布向量Xt1、Xt2，則根據(jù)式（11）可以反求出變量Px，進(jìn)而計算狀態(tài)轉(zhuǎn)移矩陣P。

4.3 大壩運行風(fēng)險度轉(zhuǎn)異判據(jù)

根據(jù)大壩運行風(fēng)險度確定大壩工程從正常轉(zhuǎn)為異常的征兆是一個涉及多方面的復(fù)雜技術(shù)問題。在水利工程結(jié)構(gòu)運行風(fēng)險度分析的基礎(chǔ)上，結(jié)合現(xiàn)有規(guī)范，嘗試建立大壩運行風(fēng)險度轉(zhuǎn)異判據(jù)。

大壩運行風(fēng)險度評判準(zhǔn)則將大壩運行風(fēng)險度分為4個等級，各運行風(fēng)險度等級對應(yīng)的運行風(fēng)險度值見表2。

表2 大壩運行風(fēng)險度等級劃分表Table 2 Dam operation risk classification

根據(jù)《水利大壩安全鑒定辦法》，初建大壩運行風(fēng)險度轉(zhuǎn)異判據(jù)如下：

（1）當(dāng)RI≤0.31時，對應(yīng)《水利大壩安全鑒定辦法》，該大壩處于“一類壩”，認(rèn)為大壩運行正常，只需采取日常維護(hù)措施即可；

（2）當(dāng)0.31＜RI≤0.64時，對應(yīng)《水利大壩安全鑒定辦法》，該大壩處于“二類壩”，認(rèn)為大壩從正常轉(zhuǎn)為帶“病”運行，此時大壩需要進(jìn)行維修處理；

（3）當(dāng)RI＞0.64時，對應(yīng)《水利大壩安全鑒定辦法》，該大壩處于“三類壩”，認(rèn)為大壩由“病壩”轉(zhuǎn)異為“險壩”，此時需要對大壩進(jìn)行加固處理。

實際運行情況中，當(dāng)某一大壩被判定為“三類壩”時，經(jīng)過除險加固后大壩還能夠運行。吳中如[7]將三類壩的80%劃歸為高風(fēng)險，剩余的20%劃歸為惡性風(fēng)險，即高風(fēng)險對應(yīng)的運行風(fēng)險度區(qū)間為（0.64,0.93]，惡性風(fēng)險對應(yīng)的運行風(fēng)險度區(qū)間為（0.93,1.00]，因此，認(rèn)為當(dāng)大壩的運行風(fēng)險度RI＞0.93時，大壩處于風(fēng)險度極高的狀態(tài)，不適宜繼續(xù)運行。

4.4 大壩運行風(fēng)險度轉(zhuǎn)異診斷過程

不考慮除險加固效應(yīng)的大壩運行風(fēng)險度轉(zhuǎn)異診斷基本步驟如下：

（1）選取若干段大壩原型監(jiān)測資料，一般以每5年的監(jiān)測資料為一個評估時段；

（2）根據(jù)大壩原型監(jiān)測資料建立統(tǒng)計模型，并分離時效分量；

（3）擬定大壩安全閾值，計算大壩單項、局部和綜合運行風(fēng)險度；

（4）根據(jù)大壩單項運行風(fēng)險度，由式（11）計算狀態(tài)轉(zhuǎn)移矩陣P，建立單項運行風(fēng)險時變模型，進(jìn)而建立局部運行風(fēng)險度和綜合運行風(fēng)險度時變模型；

（5）根據(jù)大壩綜合運行風(fēng)險度時變模型以及大壩運行風(fēng)險度轉(zhuǎn)異判據(jù)，評估大壩在全壽命周期內(nèi)的風(fēng)險度等級。

5 實例分析

以某混凝土大壩為例，分析該壩的運行風(fēng)險度和運行壽命。該混凝土壩位于第二松花江干流上，于1937年修建，1951～1953年進(jìn)行了擴(kuò)建和改建，1953年全部建成。樞紐工程由混凝土重力壩、溢流壩、壩后式廠房及變電站、左岸泄洪洞（后改建成三期引水洞）及三期發(fā)電廠房組成。攔河壩最大壩高91.7 m，壩頂長度1 080 m，壩頂高程267.70 m，共有60個壩段。目前電站總裝機(jī)容量為1 002.5 MW，多年平均發(fā)電量19.68×108kW·h。水庫正常蓄水位263.50 m，汛限水位260.50 m，死水位242.00 m，校核洪水位267.70 m，總庫容109.88×108m3，為多年調(diào)節(jié)水庫。

5.1 大壩運行風(fēng)險辨識

初步確定影響大壩運行的風(fēng)險因素主要包括變形、強(qiáng)度、老化、滲流等，如圖5所示。

圖5 某大壩運行風(fēng)險評價粗集Fig.5 Rough set of operation risk assessment for a dam

參考類似工程，相應(yīng)于圖5的四個判斷矩陣為：

采用文獻(xiàn)[13]的方法對大壩運行風(fēng)險因子進(jìn)行權(quán)值排序，計算結(jié)果見表3。

從表3得到的各風(fēng)險因素總排序權(quán)值來看，影響該大壩安全運行的風(fēng)險因素主要包括壩體變形（ω1=0.416）、壩體強(qiáng)度（ω3=0.170）和壩基滲流（ω3=0.169）。

5.2 基于實測資料的大壩運行風(fēng)險度計算分析

根據(jù)該大壩的實際情況，其安全監(jiān)控重點壩段包括14號、22號和35號壩段，其中，14號壩段壩基存在斷層，破碎帶寬度超過5 m，此外還存在壩頂向上游時效變形、壩頂上抬較明顯、壩體裂縫較嚴(yán)重等問題；22號壩段存在壩基揚壓力大、壩基漏水量較大、壩體混凝土質(zhì)量較差等問題；35號壩段是該大壩監(jiān)控的核心壩段，其重要性遠(yuǎn)大于其他壩段，而35號壩段處于壩址區(qū)最大的斷層破碎帶中心區(qū)域上，建基面物理力學(xué)參數(shù)偏低，抗滑承載力不足，壩基穩(wěn)定安全系數(shù)低，此外，35號壩段還存在壩頂位移大、壩頂向上游時效變形、壩體混凝土強(qiáng)度低、壩基揚壓力大等問題。

表3 某大壩運行風(fēng)險因素總排序權(quán)值計算結(jié)果Table 3 Weights of operation risk factors for a dam

根據(jù)運行風(fēng)險辨識的結(jié)果和重點壩段存在的問題，收集如下實測數(shù)據(jù)資料用于建模分析：（1）14號壩段壩頂垂直位移y1、壩基滲漏量y2資料；（2）22號壩段的壩基滲漏量y3、壩基揚壓力y4資料；（3）35號壩段的壩頂水平位移y5、壩體滲漏量y6、壩基橫向揚壓力y7資料。根據(jù)各測點2001年1月～2010年9月的觀測數(shù)據(jù)，建立統(tǒng)計模型并分離時效分量。風(fēng)險度計算第一個評價時段選取2001～2005年的監(jiān)測資料，第二個評價時段選取2006～2010年的監(jiān)測資料。數(shù)值表現(xiàn)和趨勢表現(xiàn)的權(quán)重向量W=[0.7,0.3]，2001～2005年大壩單項、局部運行風(fēng)險度計算結(jié)果見表4，2006～2010年大壩單項、局部運行風(fēng)險度計算結(jié)果見表5。

由大壩的局部運行風(fēng)險度計算公式可以計算大壩綜合運行風(fēng)險度，其中2001～2005年大壩綜合運行風(fēng)險度為0.494，大壩處于“中風(fēng)險”狀態(tài)。2006～2010年大壩綜合運行風(fēng)險度為0.634，大壩處于“中風(fēng)險”狀態(tài)。

表4 2001～2005年某大壩單項、局部運行風(fēng)險度計算結(jié)果Table 4 Calculation results of single and local operation risk degree from 2001 to 2005

表5 2006～2010年某大壩運行單項、局部風(fēng)險度計算結(jié)果Table 5 Calculation results of single and local operation risk degree from 2006 to 2010

現(xiàn)計算各因素綜合模糊向量的Markov狀態(tài)轉(zhuǎn)移矩陣，由上述兩個時段綜合模糊隸屬向量和式（11）可計算大壩各因素結(jié)構(gòu)運行風(fēng)險從當(dāng)前狀態(tài)向下一狀態(tài)轉(zhuǎn)移的概率Px，見表6。

由狀態(tài)轉(zhuǎn)移矩陣Px，可計算任意時刻該大壩的綜合模糊隸屬向量，從而可計算任意時刻該大壩的單項、局部和綜合運行風(fēng)險度。其中，各因素單項運行風(fēng)險度時變圖見圖6，各壩段局部運行風(fēng)險度時變圖見圖7，大壩綜合運行風(fēng)險度見圖8。由圖8可以看出，該大壩再運行8年之后（2013年）將處于高風(fēng)險狀態(tài)，有必要進(jìn)行除險加固處理，若不考慮除險加固效應(yīng)，該大壩在第34年后處于惡性風(fēng)險狀態(tài)，大壩不適宜繼續(xù)服役。

表6 大壩各因素運行風(fēng)險轉(zhuǎn)移概率Table 6 Transition probability of dam operation risk factors

圖6 大壩單項運行風(fēng)險度時變圖Fig.6 Time-varying diagram of dam single operation risk degree

圖7 大壩局部運行風(fēng)險度時變圖Fig.7 Time-varying diagram of dam local operation risk degree

圖8 大壩綜合運行風(fēng)險度時變圖Fig.8 Time-varying diagram of dam comprehensive operation risk degree

6 結(jié)語

大壩運行壽命評估是一個復(fù)雜的多技術(shù)問題，基于大壩原型觀測數(shù)據(jù)，利用馬爾科夫鏈理論，初步探究了大壩運行壽命的評估方法，并將該方法用于某大壩運行壽命評估中，在對該壩原型監(jiān)測資料深入剖析的基礎(chǔ)上，分析評價了該壩的運行風(fēng)險狀況和時變過程，其結(jié)果表明：

（1）根據(jù)不同監(jiān)測項目各個測點的實際監(jiān)測資料，可以有效評估大壩的運行風(fēng)險度，判斷大壩的運行性態(tài)；

（2）引入馬爾科夫理論可以較快速和方便地建立大壩運行風(fēng)險度時變模型，并通過風(fēng)險度轉(zhuǎn)異判據(jù)實現(xiàn)大壩工程全壽命周期內(nèi)的風(fēng)險度等級評定和工程除險加固時間的確定，從而為制定工程除險加固策略提供依據(jù)。