基于蒙特卡羅抽樣和最大-最小法的地震易損度算法研究

喻章程，肖 軍，李肇華，詹文輝

(1.上海核工程研究設計院有限公司，上海 200233；2.生態(tài)環(huán)境部 核與輻射安全中心，北京 100082)

福島核事故后，核電行業(yè)逐步重視地震對核電廠造成的影響。地震概率安全評價(PSA)是使用一套概率論的邏輯方法對地震風險進行評估[1-3]，定量地給出地震風險，其關鍵技術環(huán)節(jié)包括地震危險性分析、地震易損度評估、系統(tǒng)響應分析。地震危險性分析的目的是獲取廠址不同強度地震動(如峰值地面加速度，PGA)的發(fā)生頻率，確定地震始發(fā)事件發(fā)生頻率[4-6]。地震易損度評估的目的是評估重要構筑物和設備由于地震失效的條件概率[7-9]。系統(tǒng)響應分析環(huán)節(jié)可分為地震PSA模型建立和地震風險定量化兩部分。地震風險定量化需要整合地震危險性曲線和地震易損度曲線，以定量化分析地震事件序列的發(fā)生頻率。

易損度計算是地震風險定量化的重要組成部分。文獻[1]介紹了3種地震易損度算法，分別為蒙特卡羅模擬方法、拉丁超立方抽樣方法、離散概率分布方法。然而并未提供多個設備組合地震易損度的算法，本文探索并研究一種基于蒙特卡羅(MC)抽樣和最大-最小法計算單個基本事件、最小割集(MCS)和事故序列易損度的方法，為工程應用中多個設備組合地震易損度計算提供另一種可行的算法。

1 MC基本理論

MC方法是一種隨機抽樣技術[10]，在核物理領域有廣泛應用。其基本思想是：當所求問題的解是某個事件的概率，或是某個隨機變量的數(shù)學期望，或是與概率、數(shù)學期望有關的量時，通過某種試驗方法，得出該事件發(fā)生的頻率，或該隨機變量若干具體觀察值的算術平均值，通過它得到問題的解。

以常見的射擊問題為例，設r為射擊運動員彈著點到靶心的距離，g(r)為擊中r處相應的得分數(shù)(環(huán)數(shù))，f(r)為該運動員彈著點的分布密度函數(shù)，反映運動員的射擊水平。假設射擊運動員的射擊環(huán)數(shù)0～10環(huán)對應的概率依次為0、0、0、0、0、0、0、0.1、0.1、0.3、0.5，則進行隨機試驗(射擊)的方法為選取1個隨機數(shù)ξ，根據(jù)ξ值可判斷得到命中環(huán)數(shù)。這樣就進行了1次隨機試驗(射擊)，得到了1次成績g(r)，進行N次試驗后，得到該運動員射擊成績的近似值：

(1)

由上例可看出，MC方法常以一個“概率模型”為基礎，按照它所描述的過程，使用由已知分布抽樣的方法，得到部分試驗結果的觀察值，從而求得問題的近似解。

2 單個基本事件易損度計算

地震易損度指在給定地震動參數(shù)(如峰值地面加速度或峰值譜加速度)的情況下設備或構筑物的條件失效概率[2]。地震易損度分析認為地面加速度是滿足一定概率分布的隨機變量，在分布參數(shù)、分布形狀和設備的失效模式方面存在一定的不確定性。在給定設備失效模式和分布參數(shù)的情況下，能獲得1條表示隨地面峰值加速度變化的條件失效概率曲線。因此對于不同的參數(shù)假設，能獲得不同的易損度曲線。一種合理的方式是通過1簇設備易損度曲線來體現(xiàn)概率分布的不確定性。地震PSA通常假設設備和構筑物易損度服從雙對數(shù)正態(tài)分布，采用下式[2]計算：

(2)

其中：F為失效概率；Am為易損度中值；βR和βU分別為Am的隨機不確定性參數(shù)和認知不確定性參數(shù)；Q為置信度；Φ-1(·)為標準正態(tài)分布的反函數(shù)。

構筑物和設備易損度還可采用均值表示：

(3)

(4)

(5)

其中，ξ1、ξ2、ξ3、ξ4為(0，1)之間的隨機數(shù)。采用式(4)、(5)可求出以置信度表示的易損度曲線。

圖1 單個基本事件易損度抽樣結果與理論值對比Fig.1 Comparison of fragility sampling result and theoretical value for basic event

從圖1可看出，MC抽樣曲線與理論曲線符合得很好。采用MC方法計算得到的高置信度低概率失效(HCLPF)值為0.71g，與理論值0.70g一致。另外也可采用三維曲面圖直觀表示設備和構筑物的易損度曲線簇。該設備易損度累積概率分布示于圖2。

圖2 易損度累積概率分布Fig.2 Cumulative probability distribution of fragility

3 最小割集易損度計算

由于地震PSA定量化的特殊性，需考慮以下3種最小割集類型：1) 純地震失效最小割集，包含的基本事件均為不同設備的地震失效模式，且不含成功邏輯基本事件(類型1)；2) 包含非事件的最小割集，包含的基本事件均為不同設備的地震失效模式，且至少有1個基本事件為成功邏輯(類型2)；3) 地震失效和隨機失效混合割集，包含的基本事件含有不同設備的地震失效模式和隨機失效模式(類型3)。

對于類型1，假設最小割集形式為{B1，B2}，其中B1和B2都是地震導致失效的基本事件，可代表設備或構筑物。B1的易損度三參數(shù)分別為Am=1.6g、βR=0.2、βU=0.3。B2的易損度三參數(shù)分別為Am=1.2g、βR=0.2、βU=0.4。

采用MC計算最小割集概率時，需首先給出抽樣方法。割集{B1，B2}的物理意義是發(fā)生某一地震動水平(地面峰值加速度，PGA)時，B1和B2同時失效則最小割集發(fā)生。根據(jù)以上物理意義，結合最大-最小法[11]，可制定圖3所示抽樣方法。其中，A為一次抽樣的地震動水平(樣本)。

圖3 割集{B1，B2}的抽樣方法Fig.3 Sampling method for MCS {B1, B2}

圖1所示的抽樣過程可簡化為A=max(AB1，AB2)，即1次MC抽樣的地震動水平為AB1和AB2的最大值。該抽樣方法同樣可擴展到最小割集含有3個以上基本事件的情況，即A=max(AB1，AB2，…，ABn)。將所有樣本按升序排列即可得到此割集的易損度概率密度分布。本文以割集{B1，B2}為例計算該割集的均值易損度曲線，如圖4所示。

圖4中也示出了割集{B1，B2}的理論曲線，計算方法為首先計算0～5gPGA區(qū)間內，B1和B2的均值易損度，區(qū)間間隔為0.01g；然后在每個區(qū)間間隔內將B1和B2的易損度相乘，即得到理論結果。從圖4可看出，割集{B1，B2}的MC抽樣結果與理論結果基本一致。

圖4 割集{B1，B2}抽樣結果Fig.4 Sampling result for MCS {B1, B2}

圖5 割集抽樣結果Fig.5 Sampling result for MCS {B1,

對于地震失效模式，基本事件條件失效概率隨PGA的升高而增大，其成功邏輯的發(fā)生概率逐漸接近0，因此會對最終結果產生影響。從圖5可看出，在0～1g區(qū)間內，由于B2的地震失效條件概率很低，因此最小割集易損度曲線與B1的易損度曲線(即忽略B2成功概率的最小割集)基本相同；當進入PGA大于1.0g的區(qū)域時，B2成功概率對最小割集結果的作用明顯增加，最小割集條件失效概率在越過峰值后逐漸降低到0。由此可得，一般情況下，地震失效模式的成功概率在PGA較大時會顯著降低最小割集的條件失效概率。

類型3通常出現(xiàn)在故障樹/事件樹模型含有地震失效模式的SCDF最小割集中。該模型與SET連接，且通常基于內部事件模型建立。假設最小割集形式為{B1，B2，R1}，其中R1表示某設備的隨機失效基本事件。計算時可將最小割集分為純地震失效{B1，B2}和隨機失效{R1}兩部分。首先采用類型1的方法計算{B1，B2}部分，然后將每個PGA對應的累積失效概率與R1對應的失效概率相乘，即可得到{B1，B2，R1}的累積條件失效概率。

本文假設R1的失效概率為P(R1)=1.0×10-2，計算了最小割集{B1，B2，R1}的易損度均值曲線，如圖6所示。

圖6 割集{B1，B2，R1}抽樣結果Fig.6 Sampling result for MCS {B1, B2, R1}

4 事故序列易損度計算

如果事故序列僅包括類型1、2或3的最小割集，則可根據(jù)前文介紹的定量化方法分別進行計算。如果事故序列包括多種類型最小割集，則處理方式會更復雜。

對于類型1和類型2組成的最小割集集合，采用本算法計算時，首先需將割集集合分為2個子集，第1個子集包含的最小割集全部屬于類型1，第2個子集包含的最小割集全部屬于類型2。

采用本算法可很容易地對第1個子集進行定量化。子集中的所有最小割集均為“或”的關系，其物理意義是只要地震動水平達到所有最小割集的最低抗震能力，則事故發(fā)生。因此抗震能力最低的最小割集即為該子集的抗震能力，即A=min(MCS1，MCS2，…，MCSn)，其中MCSi表示每個最小割集的抗震能力。MCSi可根據(jù)前文討論的方法進行抽樣，即AMCSi=max(AMCSiB1，AMCSiB2，…，AMCSiBn)，其中AMCSiBn表示割集MCSi包含的基本事件。由此可知，該子集的抗震能力抽樣為：A=min(max(AMCS1B1，AMCS1B2，…，AMCS1Bn)，max(AMCS2B1，AMCS2B2，…，AMCS2Bn)，…，max(AMCSnB1，AMCSnB2，…，AMCSnBn))。

上述抽樣過程與抗震裕度評價(SMA)中計算事故序列易損度時采用的最小-最大法十分類似，不同之處在于，SMA中的基本事件為95%置信度5%失效概率的點的估計值，因此無需進行多次抽樣；而地震PSA中的地震失效基本事件為隨機變量，需采用MC方法進行多次抽樣。

對于第2個子集，可采用前文介紹的方法計算其中的每個最小割集，然后利用極限近似方法(MCUB)計算該子集的易損度曲線，最后將第1、2個子集形成的易損度曲線相加即可得到類型1和類型2組成的最小割集集合的易損度曲線。

MCUB公式為：

(6)

其中：P為事故序列失效概率；P(MCSi)為第i個最小割集失效概率。

根據(jù)以上分析，對于僅包括地震失效的事故序列(特別是對于采用地震前置樹的建模方式，所有地震失效僅出現(xiàn)在地震前置樹中，不會出現(xiàn)在事故緩解故障樹中)，當事故序列較復雜時，難以采用解析方式有效獲得事故序列的易損度，一般需采用近似算法計算事故序列的易損度，如MCUB等。但若采用本文介紹的方法計算其易損度時，不需再采用MCUB或其他割集定量化算法，僅通過比較基本事件和最小割集的抗震能力即可獲得事故序列的抗震能力，進而獲得事故序列的易損度曲線，且計算結果精度高。

5 結論

本文探索并研究了一種新的地震易損度算法，該算法基于蒙特卡羅抽樣方法和最大-最小法計算單個基本事件、最小割集和事故序列割集的易損度。經對比，該方法的結果與理論值一致，為工程應用中的地震易損度計算提供了另一種可行的算法。