基于蒙特卡罗抽样和最大-最小法的地震易损度算法研究

2022-07-29喻章程李肇华詹文辉

原子能科学技术 2022年7期

喻章程，肖军，李肇华，詹文辉

(1.上海核工程研究设计院有限公司，上海 200233；2.生态环境部核与辐射安全中心，北京 100082)

福岛核事故后，核电行业逐步重视地震对核电厂造成的影响。地震概率安全评价(PSA)是使用一套概率论的逻辑方法对地震风险进行评估[1-3]，定量地给出地震风险，其关键技术环节包括地震危险性分析、地震易损度评估、系统响应分析。地震危险性分析的目的是获取厂址不同强度地震动(如峰值地面加速度，PGA)的发生频率，确定地震始发事件发生频率[4-6]。地震易损度评估的目的是评估重要构筑物和设备由于地震失效的条件概率[7-9]。系统响应分析环节可分为地震PSA模型建立和地震风险定量化两部分。地震风险定量化需要整合地震危险性曲线和地震易损度曲线，以定量化分析地震事件序列的发生频率。

易损度计算是地震风险定量化的重要组成部分。文献[1]介绍了3种地震易损度算法，分别为蒙特卡罗模拟方法、拉丁超立方抽样方法、离散概率分布方法。然而并未提供多个设备组合地震易损度的算法，本文探索并研究一种基于蒙特卡罗(MC)抽样和最大-最小法计算单个基本事件、最小割集(MCS)和事故序列易损度的方法，为工程应用中多个设备组合地震易损度计算提供另一种可行的算法。

1 MC基本理论

MC方法是一种随机抽样技术[10]，在核物理领域有广泛应用。其基本思想是：当所求问题的解是某个事件的概率，或是某个随机变量的数学期望，或是与概率、数学期望有关的量时，通过某种试验方法，得出该事件发生的频率，或该随机变量若干具体观察值的算术平均值，通过它得到问题的解。

以常见的射击问题为例，设r为射击运动员弹着点到靶心的距离，g(r)为击中r处相应的得分数(环数)，f(r)为该运动员弹着点的分布密度函数，反映运动员的射击水平。假设射击运动员的射击环数0～10环对应的概率依次为0、0、0、0、0、0、0、0.1、0.1、0.3、0.5，则进行随机试验(射击)的方法为选取1个随机数ξ，根据ξ值可判断得到命中环数。这样就进行了1次随机试验(射击)，得到了1次成绩g(r)，进行N次试验后，得到该运动员射击成绩的近似值：

(1)

由上例可看出，MC方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，从而求得问题的近似解。

2 单个基本事件易损度计算

地震易损度指在给定地震动参数(如峰值地面加速度或峰值谱加速度)的情况下设备或构筑物的条件失效概率[2]。地震易损度分析认为地面加速度是满足一定概率分布的随机变量，在分布参数、分布形状和设备的失效模式方面存在一定的不确定性。在给定设备失效模式和分布参数的情况下，能获得1条表示随地面峰值加速度变化的条件失效概率曲线。因此对于不同的参数假设，能获得不同的易损度曲线。一种合理的方式是通过1簇设备易损度曲线来体现概率分布的不确定性。地震PSA通常假设设备和构筑物易损度服从双对数正态分布，采用下式[2]计算：

(2)

其中：F为失效概率；Am为易损度中值；βR和βU分别为Am的随机不确定性参数和认知不确定性参数；Q为置信度；Φ-1(·)为标准正态分布的反函数。

构筑物和设备易损度还可采用均值表示：

(3)

(4)

(5)

其中，ξ1、ξ2、ξ3、ξ4为(0，1)之间的随机数。采用式(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的易损度三参数分别为Am=1.6g、βR=0.2、βU=0.3。B2的易损度三参数分别为Am=1.2g、βR=0.2、βU=0.4。

采用MC计算最小割集概率时，需首先给出抽样方法。割集{B1，B2}的物理意义是发生某一地震动水平(地面峰值加速度，PGA)时，B1和B2同时失效则最小割集发生。根据以上物理意义，结合最大-最小法[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区间内，B1和B2的均值易损度，区间间隔为0.01g；然后在每个区间间隔内将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的升高而增大，其成功逻辑的发生概率逐渐接近0，因此会对最终结果产生影响。从图5可看出，在0～1g区间内，由于B2的地震失效条件概率很低，因此最小割集易损度曲线与B1的易损度曲线(即忽略B2成功概率的最小割集)基本相同；当进入PGA大于1.0g的区域时，B2成功概率对最小割集结果的作用明显增加，最小割集条件失效概率在越过峰值后逐渐降低到0。由此可得，一般情况下，地震失效模式的成功概率在PGA较大时会显著降低最小割集的条件失效概率。

类型3通常出现在故障树/事件树模型含有地震失效模式的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的最小割集，则可根据前文介绍的定量化方法分别进行计算。如果事故序列包括多种类型最小割集，则处理方式会更复杂。

对于类型1和类型2组成的最小割集集合，采用本算法计算时，首先需将割集集合分为2个子集，第1个子集包含的最小割集全部属于类型1，第2个子集包含的最小割集全部属于类型2。

采用本算法可很容易地对第1个子集进行定量化。子集中的所有最小割集均为“或”的关系，其物理意义是只要地震动水平达到所有最小割集的最低抗震能力，则事故发生。因此抗震能力最低的最小割集即为该子集的抗震能力，即A=min(MCS1，MCS2，…，MCSn)，其中MCSi表示每个最小割集的抗震能力。MCSi可根据前文讨论的方法进行抽样，即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个最小割集失效概率。

根据以上分析，对于仅包括地震失效的事故序列(特别是对于采用地震前置树的建模方式，所有地震失效仅出现在地震前置树中，不会出现在事故缓解故障树中)，当事故序列较复杂时，难以采用解析方式有效获得事故序列的易损度，一般需采用近似算法计算事故序列的易损度，如MCUB等。但若采用本文介绍的方法计算其易损度时，不需再采用MCUB或其他割集定量化算法，仅通过比较基本事件和最小割集的抗震能力即可获得事故序列的抗震能力，进而获得事故序列的易损度曲线，且计算结果精度高。