随着世界经济和城市建设的快速发展，变电站系统的安全性日益受到关注. 然而，进入21世纪后，全球，尤其是我国的地震强度和频度均有所增强，严峻的地震灾情对我国防震救灾工作提出了更高的要求. 在过去的30 a中，地震引起的变电站设施（结构、设备）破坏已经在世界范围内严重阻碍了社会经济发展^[1-2]. 变电站是维持社会、经济正常运转的重要基础设施之一，也是城市遭遇地震灾害的主要承灾体，探明其功能震损机理对变电站抗震可靠性研究具有重大意义. 目前，国内外学者已经对结构、设备及变电站系统开展了大量的抗震可靠性研究.

在变电站系统抗震可靠性分析方面，Hwang等^[3]通过二分试验法将系统的每一种状态以状态向量表示，系统总状态数量为2ⁿ，计算量将随着单元数量的增加而急剧增加. Volkanovski等^[4]通过故障树建立了变电站系统的基本模型，所建模型能够清楚、完整地表达系统单元之间的逻辑关系，并通过最小割法（minimum cut set，MCS）对建立的故障树进行化简^[5]. 考虑到变电站中存在的冗余性，李吉超^[6]提出基于概率的状态树方法，用以评估复杂系统的地震可靠性，状态树的建立过程实际上是在先评估过程中寻找所有的成功路径^[7]. 综上所述，目前对于变电站系统抗震可靠性的研究多聚焦于先评估过程，即考虑如何表达设备间的串并联关系，建立合理的系统模型，而忽略后评估过程中的求解效率. 邻接矩阵法作为后评估的一种典型方法，能够有效表示系统节点间的拓扑关系，实现程序化计算. 近年来已在车辆工程^[8]、电力系统、供水系统和道桥系统^[9-11]等诸多领域有所应用，然而其在变电站系统领域的相关研究仍相对欠缺.

本研究提出新的基于邻接矩阵法的变电站系统抗震可靠性分析流程，邻接矩阵法在先评估过程中建立的系统模型能够充分考虑各单元间的相互作用，而所采用的拟Warshall算法则能够大幅减少后评估过程中的计算量^[12]. 所建立的系统模型通过Quasi-Monte Carlo （QMC）模拟能够准确地计算出整个系统的失效概率. 本方法的优点是在建模过程中能够保留系统单元间的逻辑关系，结合“功能单元”概念与边权网络模型，综合考虑各供电设施的易损性，大大降低了系统网络基本模型的复杂度，使其在拟蒙特卡罗仿真中具有可计算性. 在评估过程中可以实现程序化模拟并高效得到系统的抗震可靠性分析结果.

变电站系统网络具有大量的供电设施，并且这些设施具有复杂的串并联关系. 为了研究变电站系统的抗震可靠性，首先将变电站系统抽象为一个由许多节点组成的几何图形，并通过相关的链路将这些节点连接起来. 节点代表这些供电设施，例如设备和架构，边代表忽略了其物理结构和电气参数差异的母线. 变电站系统被抽象为没有权重的有向拓扑图，以集合的形式表示如下：

(1) $ G=\left\{V\text{，}E\right\}. $

式中：非空集合 $ V = \left\{ {{v_{\text{1}}},{v_2}, \cdots ,{v_n}} \right\} $为G的节点集；集合 $ E = \left\{ {{e_{\text{1}}},{e_2}, \cdots , {e_n}} \right\} $为G的边集. 任意一个边 $ {e_i} = \left\{ {{v_s},{v_i}} \right\} $都对应集合上的一个二元关系.

为了实现程序化计算，采用拓扑图的邻接矩阵来描述其结构. 若有向拓扑图 $ G(V,E) $包含的变电站节点数为n，则用 $ n \times n $的邻接矩阵A表示节点之间的直接连接情况.

(2) $ {\boldsymbol{A}} = [{a_{{\boldsymbol{ij}}}}]. $

式中：若 $a_{ij} = 1$，表示供电设施 $ {v_i} $与供电设施 $ {v_j} $存在相邻边； $a_{ij} = 0$，表示供电设施 $ {v_i} $与供电设施 $ {v_j} $不存在相邻边.

QMC模拟方法与常规Monte Carlo （MC）模拟方法相似，但是其理论基础不同. QMC模拟的基本思路是采用具有较低离散型的低偏差序列（low discrepancy sequence, LDS）替代MC模拟方法所对应的伪随机数序列^[13-15]. 在先前电力系统抗震可靠性研究中已经证明QMC模拟方法的高效性^[16]，因此其理论基础不再详述. 对于QMC模拟方法，LDS的构造是其核心内容. LDS的特性为构造过程简便，且序列均匀性较好，通常可以采用Sobol或Halton低偏差序列. 本研究采用Sobol序列作为QMC模拟的原始抽样序列，并基于线性矩阵加扰技术将Sobol低偏差序列进行随机化^[17]，保证其具有低偏差与随机性的双重特征. 随机过程分以下2步.

1）以一个加扰动的 ${{\boldsymbol{v'}_j}}$代替二进制向量 $ {{\boldsymbol{v}}_j} $，有

(3) $ {{{\boldsymbol{v}}}_{j}^{\prime (t)}}^{}={{{v}}}_{j,{1}}^{(t)} {{\boldsymbol{I}}}_{{1}}^{(t)}+{{{v}}}_{j,{2}}^{(t)} {{\boldsymbol{I}}}_{{2}}^{(t)}+\cdots . $

式中：t为Sobol序列的维度变量， $ {\boldsymbol{I}}_k^{(t)} $为以2为基的二进制向量.

2）对每个新的向量增加一个二进制相互独立的随机向量 $ {{\boldsymbol{e}}^{(t)}} $：

(4) $ {\boldsymbol{r}}_n^{(t)}{\text{ = }}{{\boldsymbol{e}}^{(t)}} \oplus {g_0}(n){\boldsymbol{v'}}_{\text{1}}^{(t)} \oplus {g_{\text{1}}}(n){\boldsymbol{v'}}_{\text{2}}^{(t)} \oplus \cdots \oplus {g_{n - 1}}(n){\boldsymbol{v'}}_n^{(t)}. $

式中： ${g_{n - 1}}(n){{\boldsymbol{v}}_n}$为格莱码的二进制表示形式，以n－1作为变量.

引入连通性矩阵的概念，连通性矩阵M指用矩阵形式来描述有向图中各节点之间经过一定数量边相连的状况， $m_{ij} = 1$与终点 $ {v_j} $之间存在通路. 与 $a_{ij}$不同之处在于 $m_{ij}$表示点 $ {v_i} $与点 $ {v_j} $之间存在通路的状态，通路可能由多条边组成. 在电力系统抗震连通可靠性的分析过程中，须求得假定的发电源点到每一个变电站节点的连通性关系，传统方法基于邻接矩阵的阶乘运算可以实现. 例如基于图论的矩阵阶乘法和模糊数学法，对邻接矩阵或模糊矩阵分别阶乘后，进行布尔加运算求解连通性，通常计算量随着矩阵阶数的增长而快速增加. Warshall算法是一种简单而广泛使用的算法，用于计算边缘加权定向图中的所有节点之间的最短路径. 其优点在于通过对邻接矩阵A的元素直接进行若干次布尔运算后即可求得节点间的最短路径. Warshall算法虽然原理简单、容易实现，但是变电站系统网络拓扑的邻接矩阵大多数为稀疏矩阵，在每次遍历过程中会产生大量的非0元素. 此外，在功能性分析中并不关心2点间的最短路径. 因此，本研究提出适用于变电站系统抗震可靠性分析的拟Warshall算法，其基本思路如下.

1）对于 $ G=\left\{V，E\right\} $有向拓扑图， $ V=\{{v}_{1}，{v}_{2}，\cdots ， {v}_{n}\} $为G的节点集， $ E = \left\{ {{e_{\text{1}}},{e_2}, \cdots ,{e_n}} \right\} $为G的边集，节点数为n，其邻接矩阵A表示为

(5) $ {a}_{ij}=\left\{\begin{array}{l}1\text{，}{v}_{i}邻接{v}_{j}\text{；}\\ 0\text{，}{v}_{i}不邻接{v}_{j} 或i=j.\end{array}\right. $

2）置新矩阵 $ {\boldsymbol{M}} = {\boldsymbol{A}} $；

3）定义所需输出的矩阵元素坐标，置其所在行号为s；

4）遍历s行的每一列j，若判定 $m_{sj}=1$，则对符合条件的j列， $ k=1，2，\cdots ，n $，有表达式如下：

(6) $ m_{sk}=m_{sk} \vee m_{jk}. $

5）重复检索，若存在未被选择的j，使得 $m_{sj}= 1$，则挑选出符合条件的j，对 $ k=1，2， \cdots ，n $执行式（6），否则停止，输出矩阵 ${\boldsymbol{M}}$.

本案例研究选择了一个在中国分布广泛的典型6进线10出线220/110/10 kV变电站^[18]. 其中，变电站的10 kV部分不易受到地震影响，因此假定该部分在地震后仍可正常工作. 220 kV部分中有6条进线，110 kV部分中有10条出线. 每条线路独立工作并连接至管母线，管母线用于收集和分配电能. 变电站系统各设备布局如图1所示. 图中，SF为变电站架构，用于支撑进出线路，DS、CT、CB、IS、PT、LA、BUS、TR分别对应为隔离开关、电流互感器、断路器、绝缘子、电压互感器、避雷器、管母线和变压器. 本研究的目的是验证所提出的可靠性分析流程的有效性，供电设施的初始易损性参数均来自震害统计，对震害统计处理后绘制的易损性曲线如图2所示. 图中，PGA为地震峰动值加速度，P为超越概率.

变电站功能F_s表示110 kV部分中10条出线的可靠性. 为了与功能单元的易损性相一致，文献[19]定义P（F_s）为地震中不超过N条出线可靠的概率. F_s的值为0~9/10，增量为1/10. 例如F_s=1/2表示10条线路中正常运转的线路为5条，F_s=9/10表示10条线路中正常运转的线路为9条.

变电站系统存在冗余性，其在正常运行条件下不会满负荷运行（通常低于75%）. 在这种情况下，完全运行的2台变压器足以维持供电. 同时假定1条进线的负荷可以支撑2条出线. 因此，变电站的功能可以表达为

(7) $ {F_{\rm{s}}}{\text{ = min}}\;\{{N_{{\text{in}}}}{\text{/5, }}{N_{{\text{Tra}}}}{\text{/2, }}{N_{{\text{out}}}}{\text{/10}}\}{\text{.}} $

式中： ${N_{{\text{in}}}}$为正常进线的数量， $ {N_{{\text{Tra}}}} $为变压器运转的数量， ${N_{{\text{out}}}}$为正常出线的数量.

为了展示各功能单元构成，将变电站划分为4部分：220 kV部分、母联部分、110 kV部分和变压器部分，如图3~6所示. 在立面图中红线为化简后的输电导线，实际上分为A、B、C三相传输. 这些设备彼此独立，但架构由不同线路的设备共享. 此外，母联部分用于组合2条不同的管母线，如图4所示给出了220 kV母联部分的示意图，110 kV母联与220 kV母联相似，只是设备间距不同. 考虑到该变电站系统中的大量设备和结构，共计664个组件. 首先采用功能单元概念进行化简，将各类结构设备依据不同功能划分为进线单元、变压器单元、出线单元、母联单元、母线单元、连接单元. 等价功能单元的破坏概率为各组成设施的联合破坏概率：

(8) $ P = 1 - \prod\limits_{i{\text{ = 1}}}^m {(1 - {P_{{c_i}}})} . $

式中： ${P_{{c_i}}}$表示组成设施 $ {c_i} $在设定地震下的破坏概率，m为组成设施的个数.

进线单元（见图7（a））由1个PT、1个DS、1个CT、1个CB和2个SF组成. SF1-1支撑进线单元1，而SF1-2支撑进线单元2~6. 为了简单起见，功能单元中统一标记为SF1，出线单元（见图7（b））构成与进线单元类似，不再详述.

连接1单元（见图8（a））与管母线相连，电能经过进线单元，通过各自的连接单元到达管母线. 且通常在远离设备侧的220 kV的DS和BUS I之间安装IS，以保持导线与地面的安全距离. 110 kV出线的设备间距较小，因此远离设备侧的110 kV的DS与BUS间未安装IS（连接2单元）. 此外，母联单元（见图8（c））用于连接2个不同的母线单元，并由2个DS、1个CT、1个CB、1个IS和2个SF组成.

220 kV母线单元（见图9（a））由SP、IS、BUS、PT和LA组成. BUS安装在ISs的顶部，与连接单元相连用于收集电力. PT和LA分别用于监控电位和保护BUS. 220 kV母线单元共计由12个SP和IS支撑，110 kV母线单元（见图9（b））类似于220 kV母线单元，不同之处在于110 kV母线单元共享16个SP和IS.

变压器单元（见图10）由3个SF、2个CB、2个CT、2个DS、2个LA和TR本身组成. 最易破坏的构件是安装在TR上部的瓷套管. 变压器连接到220 kV线路、110 kV线路和10 kV线路，因此存在3种类型的套管：高压套管（HB）、中压套管（MB）和低压套管（LB）. TB为变压器的油箱.

逻辑关系图可以考虑组件之间的相互关系，并给出变电站系统所有可能序列的完整描述. 将变电站系统按功能单元作出逻辑关系图，如图11所示. 图中的系统包含了69个节点与若干连接线路. 为了研究变电站系统地震作用下的工作性能，必须先建立其网络基本模型.

文献[19]和[20]在对变电站供电设施震害调研中发现，连接母线具有较高的可靠性，不易发生功能性破坏，因此本研究假定所有线路完全可靠，将点权模型进一步转换为边权网络模型，进而将多输入多输出网络以单输入多输出拓扑模型表示，网络复杂性有了较大程度的降低. 同时考虑到控制室的正常运转是整个系统工作的前提，添加边“①-②”代表控制室，如图12所示.

采用邻接矩阵法建立系统网络模型，以QMC法随机模拟变电站系统各功能单元在每一次模拟中的地震破坏状态，采用改进的拟Warshall算法对变电站系统进行抗震可靠性分析，基本步骤如下.

1）确定EDP范围e（本研究以PGA作为工程需求参数）、功能区间f及总模拟次数N.

2）求解在第e次循环下各功能单元对应的破坏概率P_c,ij，并生成初始邻接矩阵A.

3）以QMC法对变电站系统中所有功能单元产生低偏差数列r_ij，并与功能单元的破坏概率P_c,ij相比较，修正本次模拟的邻接矩阵A，其中，

(9) $ {a}_{ij}=\left\{\begin{array}{*{20}{c}}1,& {r}_{ij}>{P}_{{\rm{c}},ij};\\ 0,& {r}_{ij}\leqslant {P}_{{\rm{c}},ij}.\end{array}\right. $

4）定义进线单元、变压器单元及出线单元所对应的矩阵元素坐标，并以拟Warshall算法求解得到变电站系统的功能性矩阵M. 对所有的j列，若m_sj=1，则源点（控制室）s与功能单元j之间连通，若m_sj=0，则源点s与功能单元j之间不连通.

5）统计M矩阵中进线单元、变压器单元及出线单元对应元素m_sj=1的数量，并以式（7）判别统计，将统计结果累计到50×10阶破坏状态矩阵T中，矩阵T中的元素t_ef依次代表变电站系统在不同PGA下不同功能状态出现的次数. 在本案例研究中，PGA范围为0~0.5g，间隔为0.01g，因此e=50；变电站共划分为10个功能状态，取f=10. 矩阵T中元素表达式为

(10) $ \left.\begin{gathered} {t_{1{\text{1}}}}^{\left( {n+1} \right)}{\text{ = }}{m_{1{\text{1}}}}^{(n)}{\text+}1,\begin{array}{*{20}{c}} {}&{{F_{\rm{s}}} = 0/10}; \end{array} \\ {t_{2{\text{2}}}}^{\left( {n+1} \right)}{\text{ = }}{m_{2{\text{2}}}}^{(n)}{\text+}1,\begin{array}{*{20}{c}} {}&{{F_{\rm{s}}} = 1/10}; \end{array} \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{}&{}&{} \end{array}}&{}& \vdots &{} \end{array} \\ {t_{ef}}^{\left( {n+1} \right)}{\text{ = }}{m_{{\text{e}}f}}^{(n)}{\text+}1,\begin{array}{*{20}{c}} {}&{{F_{\rm{s}}} = \left( {f - 1} \right)/10.} \end{array} \\ \end{gathered}\right\} $

6）每个e值都重复步骤1）~5）N次，统计变电站系统各功能状态出现频率. 为了求得精确值，本研究取N＝20000. 频率表达式为

(11) $ {p_{ef}} = {t_{ef}}/{{{\text{2}}0\;000}}. $

7）将各功能状态出现频率逐级累加得到变电站各功能对应的超越概率P(F_s).

基于上述相关理论，变电站系统可靠性分析流程如图13所示.

变电站系统不同功能的易损性曲线如图14所示. 可以看出，当F_s＜8/10时，相邻F_s的结果高度相似（例如，当F_s＝6/10、7/10时，抗震能力中值分别为0.183 232、0.182 501，当F_s＝0、1/10时，抗震能力中值皆为0.295 230）. F_s=5/10和F_s=4/10的易损性曲线有较大差别，体现了变电站系统的关联性和冗余性. 理论上讲，1个进线单元可以支撑2个出线单元，因此当1个进线单元破坏应导致变电站系统功能下降2/10. 但是，进线单元中存在20%的冗余，因此在第2个进线单元发生破坏时系统功能才会受到影响. 同样，变压器单元中存在50%的冗余，除非3个变压器单元中的2个破坏，否则系统功能不会受到影响.

如表1所示为变电站系统不同功能的中值能力m_c和标准偏差s_d. 可以看出，当PGA中值低于0.15g时，大多数的功能单元依然可以正常运行，变电站系统的功能处于较高值（不小于8/10）. 此时，变电站系统是冗余的，其功能取决于出线单元，任何一个出线单元的破坏将使变电站系统的功能下降1/10. 随着PGA的增加，进线单元的冗余度降为0，一个进线单元的破坏将导致变电站功能下降2/10.

为了更清楚地说明这一点，如图15所示为变电站系统随PGA变化的不同功能的频率f. 当PGA＜0.15g时，变电站仍保持完全运行状态，而当PGA＞0.23g时，变电站丧失大半功能. 可以观察到，F_s=6/10并没有随着F_s=8/10，F_s=4/10，F_s=2/10呈现逐级下降规律，F_s=6/10下的峰值概率远小于F_s=4/10，表明变压器单元的易损性过大导致其冗余性并不能有效保证变电站的抗震可靠性. 因此，这组曲线清楚地表达了组件冗余对系统功能的影响. 结果表明，系统设计中考虑的冗余度不足以抵御0.23g的最大可信地震. 同时，当PGA较小时，变电站功能取决于出线单元及110 kV连接单元，当PGA较大时，变电站的功能主要取决于进线单元、220 kV连接单元与主变压器单元. 这是由于当PGA较大时，即便进线单元与变压器单元存在一定的冗余度，但进线单元与变压器单元的易损性远大于出线单元和110 kV连接单元.

通过逐个改变相关设备的易损性参数（能力中值提高20%），得到了相应的变电站系统易损性曲线，为了便于分析，本研究仅给出变电站系统抗震可靠性曲线的中值能力和对数标准差在不同设备类型下的变化趋势，如图16~19所示. 图中，P_m、P_s分别为中值能力和对数标准差的变化率.

由图16~17可知，110 kV设备和IS对变电站系统易损性曲线的中值能力和对数标准差影响不大（总体小于1.5%），因为这些设备具有相对较高的能力中值. 可以发现，除了110 kV的LA、IS，110 kV的CB、DS、CT、PT在F_s=9/10时对变电站的可靠性呈现整体上升，这是因为进线单元、220 kV连接单元和变压器单元在PGA较小时才是冗余的，而110 kV连接单元与出线单元包含1个CB、1个PT、1个CT和3个DS，因此110 kV的DS的变化率最大. 提高110 kV的DS的抗震性能可以显著提高变电站系统F_s=9/10的抗震可靠性. 此时，变电站功能状态由110 kV连接单元与出线单元决定，这与3.1节中的结果可以相互验证.

由图18可知，220 kV设备中CT对变电站系统易损性曲线的中值能力影响最大，随着F_s的减小在4%~7%变化. 随后是CB，中值能力随F_s的减小在3%~6%变化. 其余设备对中值能力的影响由大到小分别为PT、DS、TR、LA. 当F_s=9/10时220 kV设备的变化率显著减小，也验证了3.1节中的分析，即在系统抗震可靠性较高时，出线单元决定了系统的震后功能.

由图19可知，CB、CT、DS、PT对系统易损性曲线的对数标准差也有较大影响，但影响方式不同，即CB、CT减小对数标准差，并与F_s正相关，DS、PT增大对数标准差，与F_s成负相关. LA与PT、TR由于其较高的能力中值及较少的参与度，对变电站系统能力中值（<1.5%）的影响非常有限. 表明220 kV的LA与PT、TR单体设备具有较高的抗震冗余性. 因此，在进线单元、220 kV连接单元及变压器单元中，220 kV的CT、CB、DS被认为是变电站中最关键的功能设备，提高他们的抗震性能可以显著提高整个变电站系统的抗震可靠性.

（1）提出一种变电站系统抗震可靠性评估流程. 该流程利用邻接矩阵对系统拓扑和功能进行分析，引入“功能单元”概念，综合考虑各供电设施的易损性，大大降低了系统网络基本模型的复杂度，使其在拟蒙特卡罗仿真中具有可计算性. 针对变电站系统功能性计算，基于Warshall算法提出适用于各类变电站系统可靠性分析的功能性算法，并应用此算法对典型220 kV变电站系统进行抗震可靠性分析.

（2）通过关键设备分析确定220 kV的CT、CB、DS及110 kV的DS对系统功能影响最大，是变电站中最关键的设备. 本研究提出的方法为程序化评估变电站系统抗震可靠性提供了相应的思路，并为变电站系统以及电力系统抗震韧性评估奠定了研究基础.

（3）可靠性评估结果的准确性取决于输入的设施易损性参数，为了更精准地获得评估结果，后续计划开展更为详细的变电设施地震易损性研究.