基于贝叶斯网络的风力发电机故障诊断方法
2021-05-07唐小峰文永康
胡 宇,唐小峰,文永康,邹 建
(成都天奥测控技术有限公司，成都 611731)
0 引言
贝叶斯网络方法是由贝叶斯公式发展而来，贝叶斯公式是用于解决“逆概”的问题一种方法，即通过“结果”计算“原因”的概率[1]。贝叶斯方法作为概率统计学的重要分支，广泛应用于机器学习、数据挖掘、数据融合、医疗诊断、邮件过滤和产品客户需求分析等场景。故障诊断的过程通常是利用测试“结果”推理出故障“原因”，这和处理“逆概”问题的概念基本相同，因此，贝叶斯公式和贝叶斯网络也经常应用于故障诊断领域。刘钦文利用贝叶斯网络处理传统多信号模型中故障传递和故障检测中的不确定性问题[2]，文献[3]使用朴素贝叶斯算法分类提高了滚动轴承故障诊断的准确率，文献[4]提出基于重要度分级的概念优化传统贝叶斯网络模型有利于故障特性的提取，文献[5]提出基于贝叶斯网络的民机起落架系统故障诊断方法提升了诊断效率和精度。
基于知识的故障诊断方法是目前最为流行的故障诊断方法，这类方法主要是通过故障征兆和系统模型等方式来描述故障现象与故障原因之间的联系，系统模型可分为定量模型和定性模型，复杂系统的精确定性模型往往难以获得[6]。常用的故障树和相关性矩阵都是基于定性模型的故障诊断方法，故障分辨率低[7]，处理复杂系统多重并发故障的能力较差。贝叶斯网络具有强大的不确定性推理能力，可以利用历史数据建立较为精确的故障概率模型，因此在复杂系统故障诊断中具有一定的优势[8]。风能作为一种无污染、可再生的绿色能源，近年来得到了快速发展[9]，我国建立了大量的风力发电场并投入使用。经过长期运转，风力发电机维护成本越来越高，对状态监测和故障诊断的要求也越来越高，传统手段捉襟见肘。风力发电系统是集机电、机械、电气和自动控制技术为一体的复杂系统，故障原因不确定度较高，可采用贝叶斯网络方法进行故障诊断。
1 贝叶斯网络故障诊断模型
1.1 网络结构
贝叶斯网络结构用于定性的描述各网络节点之间的关联，在贝叶斯网络模型中采用“有向无环图-DAG”来实现。贝叶斯网络DAG由一系列代表节点的圆圈和代表关联的箭头组成，“有向”是指节点之间的连接是有方向的，“无环”是指从任一节点出发的连接不会再回到该节点，如图1所示。
图1 有向无环图
贝叶斯网络中的节点通常代表某个事件或随机变量的某种状态，在故障诊断模型中可以是故障原因、故障单元、故障模块、故障状态、故障现象、和测试结果等等。节点间的连线代表了因果关系，连线起端的节点为“因”，称为“父节点”，连线终端箭头指向的节点为“果”，称为“子节点”。图1中节点“A”、“B”连接并指向“C”，表示事件“A”、“B”有概率导致事件“C”发生，或者表示变量“A”、“B”的不同状态组合会使“C”进入某种状态。没有“父节点”的节点为根节点，也称为一级节点，一级节点的子节点为二级节点，以此类推n级节点的子节点为n+1级节点，当某个节点拥有多个父节点时，n取其中最大的那个值。例如：图1中节点“A”、“B”为一级节点，节点“C”为二级节点，节点“D”为三级节点，节点“E”既是二级节点“C”的子节点，也是三级节点“D”的父节点，取最大值“3”再加“1”，因此节点“E”为四级节点。
未知的贝叶斯网络结构难以通过样本数据学习的方法获得，一般采用工程经验进行构建[6]。以某照明电路为例，电路包括电池、3只保险管和指示灯，保险管“FUSE2”和“FUSE3”并联后与保险管“FUSE1”串联，当电路中所有器件正常时指示灯会点亮，如图2所示。
图2 照明电路原理图
假设三只保险的连接关系未知，但经验表明其故障状态会影响指示灯的亮/灭状态。以变量F1、F2、F3分别代表保险管“FUSE1”、“FUSE2”、“FUSE3”的故障状态，变量T代表指示灯“LED”的燃亮状态。F1、F2、F3的故障会导致T熄灭，即存在因果关系，“故障”为因，“熄灭”为果，电路的网络结构如图3所示。
图3 照明电路贝叶斯网络结构图
1.2 条件概率分布CPD
贝叶斯网络结构定性的描述了节点之间的关联，为了更加准确的描述节点的关联程度，还需要用到条件概率参数。条件概率是指事件B已经发生的条件下，事件A发生的概率，记为P(A|B)[1]。贝叶斯网络中的条件概率是指父节点发生或处于某种状态下，子节点发生或进入某种状态的概率。当子节点拥有多个父节点时，需要关注父节点全部组合状态下的子节点进入各个状态的概率，这些概率值的集合即为条件概率分布。图3所示网络结构中， F1、F2、F3是代表“FUSE1”、“FUSE2”、“FUSE3”的故障状态的变量，每个变量有两个取值，“1”为“故障”，“2”为“正常”。T是代表指示灯“灭”、“亮”的变量，共两个取值，“1”为“灭”，“2”为“亮”。当“FUSE1”、“FUSE2”、“FUSE3”均“故障”时，指示灯“灭”，即F1=1、F2=1、F3=1时T=1的概率为1，根据历史数据，照明电路的条件概率分布如表1所示。
表1 照明电路CPD
贝叶斯网络中的根节点状态存在不确定度，但它没有父节点，因此不能利用条件概率来描述根节点所处状态的概率，在贝叶斯统计推断中通过先验概率来表达根节点的状态置信度。照明电路贝叶斯网络的根节点是反映保险管故障状态的变量，保险管故障的概率取决于其可靠性，假设在使用一段时间过后，保险管“FUSE1”故障的概率为10%，无故障的概率为90%，则先验概率记为：P(F1=1)=0.1，P(F1=2)=0.9。设三只保险管型号和生产批次一致，可靠性指标完全相同，则其先验概率如表2所示。
表2 保险管先验概率
1.3 推理算法
贝叶斯推理的基本公式包括：
P(A|B)=P(AB)/P(B)[1]
(1)
P(A)=∑iP(A|Bi)/P(Bi)[1]
(2)
P(Bi|A)=P(A|Bi)/P(Bi)/∑jP(A|Bj)/P(Bj)[1]
(3)
公式(1)为条件概率公式，表示事件B条件下发生事件A的概率为事件A和事件B同时发生的概率除以事件B发生的概率。公式(2)为全概率公式，表示发生事件A的概率为事件B各种状态下发生事件A的概率乘以事件B该状态下概率的乘积之和。公式(3)为贝叶斯公式，由公式(1)、(2)推导得出，用于已知条件概率P(A|B)和先验概率P(B)时，计算后验概率P(B|A)。
贝叶斯网络模型中最重要的公式为联合概率分布公式，用于表述网络结构中所有节点事件同时发生的概率，可由公式(1)推导得出。图1所示的网络结构的联合概率表达为：
PJ=P(A)*P(B)*P(C|A,B)*P(D|C)*P(E|C,D)
(4)
此求解方法为消元法，是一种精确求解的方法，即通过求和某个节点各种状态的先验概率与条件概率的乘积，消除联合概率表达式中非待求解的节点。
例如：图3所示照明电路的联合概率分布为：PJ=P(F1,F2,F3,T) =P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)；
指示灯“LED”不亮时保险管“FUSE1”故障的概率表述为：P(F1=1|T=1)，求解过程为：
P(F1=1|T=1) =P(F1=1,T=1)/P(T=1) =
∑F1=1,T=1,F2,F3PJ/∑T=1,F1,F2,F3PJ
(5)
∑F1=1,T=1,F2,F3PJ=
∑F1=1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)
(6)
∑F1,T=1,F2,F3PJ=
∑F1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)
(7)
将{F1,T,F2,F3}依次取值{1,1,1,1}、{1,1,1,2}、{1,1,2,1}、{1,1,2,2}代入式(6),得出：
∑F1=1,T=1,F2,F3PJ=[P(F1=1)*P(F2=1)*P(F3=1)*P(T=1|F1=1,F2=1,F3=1)]+[P(F1=1)*P(F2=1)*P(F3=2)*P(T=1|F1=1,F2=1,F3=2)] +[P(F1=1)*P(F2=2)*P(F3=1)*P(T=1|F1=1,F2=2,F3=1)] + [P(F1=1)*P(F2=2)*P(F3=2)*P(T=1|F1=1,F2=2,F3=2)]=0.1*0.1*0.1*1 + 0.1*0.1*0.9*1 + 0.1*0.9*0.1*1 + 0.1*0.9*0.9*1 =0.001+0.009+0.009+0.081 = 0.1。
将{F1,T,F2,F3}依次取值{1,1,1,1}、{1,1,1,2}、{1,1,2,1}、{1,1,2,2}、{2,1,1,1}、{2,1,1,2}、{2,1,2,1}、{2,1,2,2}、代入式(7),得出：
∑F1,T=1,F2,F3PJ=∑F1,T=1,F2,F3P(F1)*P(F2)*P(F3)*P(T|F1,F2,F3)= [P(F1=1)*P(F2=1)*P(F3=1)*P(T=1|F1=1,F2=1,F3=1)] +[P(F1=1)*P(F2=1)*P(F3=2)*P(T=1|F1=1,F2=1,F3=2)] +[P(F1=1)*P(F2=2)*P(F3=1)*P(T=1|F1=1,F2=2,F3=1)] + [P(F1=1)*P(F2=2)*P(F3=2)*P(T=1|F1=1,F2=2,F3=2)] + [P(F1=2)*P(F2=1)*P(F3=1)*P(T=1|F1=2,F2=1,F3=1)] + [P(F1=2)*P(F2=1)*P(F3=2)*P(T=1|F1=2,F2=1,F3=2)] + [P(F1=2)*P(F2=2)*P(F3=1)*P(T=1|F1=2,F2=2,F3=1)] +[P(F1=2)*P(F2=2)*P(F3=2)*P(T=1|F1=2,F2=2,F3=2)] = 0.1*0.1*0.1*1 + 0.1*0.1*0.9*1 + 0.1*0.9*0.1*1 + 0.1*0.9*0.9*1 + 0.9*0.1*0.1*1 + 0.9*0.1*0.9*0 + 0.9*0.9*0.1*0 + 0.9*0.9*0.9*0 = 0.001+0.009+0.009+0.081+0.009+0+0+0 = 0.109。
将上述结果代入式(5)，得出：
P(F1=1|T=1)=0.1/0.109=91.74%，同理可求出P(F2=1|T=1)=17.43%、P(F3=1|T=1)= 17.43%，即当指标灯不亮的时候保险管“FUSE1”故障的概率为91.74%，“FUSE2”或“FUSE3”故障的概率为17.43%。
2 风力发电机贝叶斯网络MATLAB建模
2.1 网络结构建模
风力发电机故障高发组件、部件包括：冷却风扇、集电环组件、编码器组件、接地碳刷组件、发电机本体、发电机轴承、测温电阻等，冷却风扇由电机、电容、电缆、机械结构等零部件组成。集电环组件由主碳刷、集电环和其他配件等组成；编码器组件主要包括编码器、电缆、机械结构和附件等零部件；接地碳刷组件由接地碳刷和机械结构组成；发电机本体包括接线端子、绕组、冷却系统和线缆等零部件；发电机轴承包括轴承、线缆等零部件；测温电阻包括绕组测温电阻、轴承测温电阻和冷却系统测温电阻等。发电机运行过程中可实时监视电机转速、绕组温度、轴承温度和冷却系统温度等参数，用于评估系统是否正常运行。风力发电机系统贝叶斯网络结构如图4所示。
图4 风力发电机系统贝叶斯网络结构图
风力发电机系网络节点较多，人工处理较繁琐，因此采用Matlab FullBNT-1.0.7贝叶斯网络工具箱进行建模。根据FullBNT的要求将网络中的节点从小到大进行编号，且父节点的编号必须小于子节点编号，从1开始共50个节点编号如图4所示。风力发电机系统网络结构Matlab建模部分代码如下：
N = 50;
dag = zeros(N,N);
N1=1;N2=2;N3=3;N4=4;N5=5;N6=6;N7=7;N8=8;N9=9;N10=10;
N11=11;N12=12;N13=13;N14=14;N15=15;N16=16;N17=17;N18=18;N19=19;N20=20;
N21=21;N22=22;N23=23;N24=24;N25=25;N26=26;N27=27;N28=28;N29=29;N30=30;
N31=31;N32=32;N33=33;N34=34;N35=35;N36=36;N37=37;N38=38;N39=39;N40=40;
N41=41;N42=42;N43=43;N44=44;N45=45;N46=46;N47=47;N48=48;N49=49;N50=50;
dag(N1,N36)=1
dag(N2,N36)=1
dag(N3,N36)=1
dag(N4,N36)=1
dag(N5,N36)=1
………………………………
g=bnt.dag;
draw_graph(g);
其中：代码dag(Na,Nb)表示节点a指向节点b，代码运行结果如图5所示。
图5 风力发电机系统贝叶斯网络结构图
2.2 概率分布建模
在进行概率分布建模前，需要对网络中各节点的状态取值类型和数量进行定义。风力发电机系统贝叶斯网络各节点均为离散变量，节点44：“绕组温度”、节点45：“轴承温度”、节点46：“系统冷却温度”均有4个状态，分别取值“1”代表“超量程”、“2”代表“正常”、“3”代表“过热”、“4”代表“其他故障”。节点44、45、46之外的节点均只有2个状态，取值“1”代表“故障”、“2”代表“正常”。风力发电机系统贝叶斯网络中共有35个一级节点，8个二级节点，6个三级节点，1个四级节点。根据历史故障数据统计，部分1级节点先验概率如表3所示。
表3 1级节点先验概率
二级及以上节点的条件概率的数量m，取决于其父节点的数量n以及节点状态的数量k，设kn+1为该节点自身状态数量：
(8)
例如节点N41有节点N22、节点N23两个父节点，节点N41、N22、N23均只有两个状态，因此N41共有8个条件状态参数；节点50有8个父节点，5个节点有两种状态，3个节点为4种状态，其条件概率参数达4 096个。为了将条件概率分布参数导入FullBNT工具箱，需按一定的规律对各节点的条件概率参数进行排序，建立每个节点的条件概率表。条件概率表用于描述当前节点及相关节点处于不同状态组合条件下的概率，表的列从左至右以节点编号从小到大的顺序排列，最后一列为条件概率值；表的行从上之下以节点变量状态值从小到大的顺序排列。节点N41“接地碳刷组件”的条件概率表如表4所示，其中第一行代表节点22、节点23、节点41依次取值{1,1,1}时的概率为0.999 99。
表4为完成其余节点的条件概率表，再将这些概率参数填入名为FLFDJCPT.xls的EXCEL表格中，便于Matlab导入，部分数据如表5所示。
表4 节点N41条件概率表
风力发电机系统条件概率分布建模Matlab部分代码如下：
discrete_nodes = 1:N;
node_sizes=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 2 2 2 2];
bnt = mk_bnet(dag,node_sizes,'names',{'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10','N11','N12','N13','N14','N15','N16','N17','N18','N19','N20','N21','N22','N23','N24','N25','N26','N27','N28','N29','N30','N31','N32','N33','N34','N35','N36','N37','N38','N39','N40','N41','N42','N43','N44','N45','N46','N47','N48','N49','N50'},'discrete',discrete_nodes);
bnt.CPD{N1} = tabular_CPD(bnt,N1,[0.08258 0.91742]);
bnt.CPD{N2} = tabular_CPD(bnt,N2,[0.00965 0.99035]);
………………………………
bnt.CPD{N35} = tabular_CPD(bnt,N35,[0.00021 0.99979]);
表5 风力发电机条件概率分布
bnt.CPD{N36} = tabular_CPD(bnt,N36,[xlsread('FLFDJCPT.xls','sheet1','A2:A65')]);
bnt.CPD{N37} = tabular_CPD(bnt,N37,[xlsread('FLFDJCPT.xls','sheet1','B2:B33')]);
bnt.CPD{N38} = tabular_CPD(bnt,N38,[xlsread('FLFDJCPT.xls','sheet1','C2:C33')]);
bnt.CPD{N39} = tabular_CPD(bnt,N39,[xlsread('FLFDJCPT.xls','sheet1','D2:D17')]);
bnt.CPD{N40} = tabular_CPD(bnt,N40,[xlsread('FLFDJCPT.xls','sheet1','E2:E65')]);
bnt.CPD{N41} = tabular_CPD(bnt,N41,[xlsread('FLFDJCPT.xls','sheet1','F2:F9')]);
bnt.CPD{N42} = tabular_CPD(bnt,N42,[xlsread('FLFDJCPT.xls','sheet1','G2:G65')]);
bnt.CPD{N43} = tabular_CPD(bnt,N43,[xlsread('FLFDJCPT.xls','sheet1','H2:H9')]);
bnt.CPD{N44} = tabular_CPD(bnt,N44,[xlsread('FLFDJCPT.xls','sheet1','I2:I65')]);
bnt.CPD{N45} = tabular_CPD(bnt,N45,[xlsread('FLFDJCPT.xls','sheet1','J2:J65')]);
bnt.CPD{N46} = tabular_CPD(bnt,N46,[xlsread('FLFDJCPT.xls','sheet1','K2:K65')]);
bnt.CPD{N47} = tabular_CPD(bnt,N47,[xlsread('FLFDJCPT.xls','sheet1','L2:L17')]);
bnt.CPD{N48} = tabular_CPD(bnt,N48,[xlsread('FLFDJCPT.xls','sheet1','M2:M5')]);
bnt.CPD{N49} = tabular_CPD(bnt,N49,[xlsread('FLFDJCPT.xls','sheet1','N2:N17')]);
bnt.CPD{N50} = tabular_CPD(bnt,N50,[xlsread('FLFDJCPT.xls','sheet1','O2:O4097')]);
代码“discrete_nodes = 1:N;”表示所有节点均为离散变量；代码“node_sizes =[x x x x x]”描述各网络节点的状态取值数量，其中节点44、45、46状态数量为4，其余节点均为2；代码“bnt.CPD{Nx} = xxx”用于导入各节点的概率参数，其中先验概率直接键入，条件概率参数较多，通过读取“FLFDJCPT.xls”文件导入。
3 实验结果与分析
3.1 故障诊断结果
在完成风力发电机贝叶斯网络结构和条件概率分布建模后，即可使用FullBNT工具箱进行推理。例如需要得出系统故障时各零部件故障的概率，则可将节点50置为1，其他节点设为隐含状态，调用“联合树推理引擎”分别求解各节点取值为“1”时的概率，Matlab代码如下：
engine = jtree_inf_engine(bnt);
evidence =cell(1,N);
evidence{N50} = 1;
[engine, loglik] = enter_evidence(engine,evidence);
fori=1:50
marg = marginal_nodes(engine,i);
NBP(i)=marg.T(1)
end
bar(NBP)
代码执行结果如图6所示。
图6 风力发电机系统零部件故障概率分布
从图6可以看出，当系统故障时，除去节点50本身，节点36“冷却风扇”故障概率最高，节点47“集电环组件”和节点1“冷却风扇电机”故障概率次之。节点36“冷却风扇”是故障概率最高的组件，节点1“冷却风扇电机”是故障概率最高的零部件。
假设节点44“绕组温度”超量程，节点46“冷却系统温度”超量程，求其他节点故障概率，则将代码“evidence{N50} = 1;”改为“evidence{N44} = 1; evidence{N46} = 1;”，执行结果如图7所示。
图7 绕组及冷却系统温度超量程故障概率分布
从图7可以看出，当节点44 “绕组温度”和节点46“冷却系统温度”同时出现超量程故障时，除去节点44和节点46本身，故障概率最高的节点依次为节点26、节点50和节点42。节点26“线缆”是节点42“发电机本体”组成零部件，节点50“风力发电机系统”为系统本身，节点26故障会导致节点42组件故障，进而导致节点50系统故障。节点26“发电机本体-线缆”是故障概率最高的零部件，故障概率为99.55%，远高于其他零部件的故障概率。
3.2 与D矩阵诊断结果的对比
3.2.1D矩阵诊断示例1
D矩阵故障诊断方法是典型的定性模型方法，其核心是建立故障原因与故障现象的依赖关系[10]，再通过故障现象推导故障原因。D矩阵依赖模型与贝叶斯网络结构类似，同样利用带有箭头的连线来表达故障原因与故障现象之间的关联。若某故障原因会导致某故障现象出现，则在它们之间通过箭头进行连接。D矩阵故障诊断通常将依赖模型转化到一个布尔矩阵中，利用“1”、“0”来代表“依赖”或“不依赖”关系。风力发电机系统节点50与其父节点的依赖矩阵如表6所示。
表6 节点50依赖矩阵
表6中第一列为故障原因的节点编号，第一行为故障现象的节点编号，单元格中的“1”表示对应行的首列的节点故障会导致对应列的首行的节点出现故障。如节点41所在行中对应列节点41和节点50的单元格为“1”，则表示节点41故障时，节点41自身和节点50会出现故障。D矩阵诊断推理基本规则为：当依赖矩阵中某一故障现象存在，则具有依赖关系的故障原因可能存在；反之，当依赖矩阵某一故障现象不存在时，则具有依赖关系的故障原因肯定不存在。如表6所示第一行的最后一列节点50出现故障时，则该列单元格为“1”的对应行的节点41、节点42、节点44、节点45、节点46、节点47、节点48、节点49、节点50故障现象均可能存在。同理，通过节点41、节点42、节点44、节点45、节点46、节点47、节点48、节点49的依赖矩阵，可以推导出他们的父节点也可能故障。因此，节点50“风力发电机系统”故障时，采用D矩阵诊断方法推理得到的结果是：节点1到节点49均可能存在故障，即系统中所有零部件都可能故障。
3.2.2D矩阵诊断示例2
对于节点44 “绕组温度”和节点46“冷却系统温度”超量程故障，D矩阵诊断过程如下。首先建立节点44、节点46及其父节点的依赖矩阵，如表7所示。
根据表7所示依赖关系，当节点44、节点46故障时，可直接推导出节点26、27、33、35、36可能故障，再通过节点36故障可以推导出节点1、2、3、4、5可能故障。因此，节点44 “绕组温度”和节点46“冷却系统温度”超量程故障，D矩阵诊断结果为：节点1、2、3、4、5、26、27、33、35、36均有可能故障。
表7 节点44、46依赖矩阵
表8 诊断结果对比表
3.2.3 诊断结果对比
通过上述两个示例可以看出，D矩阵诊断通常只能得出可能的故障原因的集合，而贝叶斯网络诊断不仅可以给出可能的故障原因，还能给出其置信度。贝叶斯网络故障诊断结果与D矩阵故障诊断结果对比如表8所示。
5 结束语
贝叶斯网络基于概率统计学理论具有严谨的理论基础，是处理不确定关系问题的首选模型，适用于复杂系统的故障诊断及维修决策。通过对风力发电机系统的建模和模拟故障诊断结果可以看出，相比于传统相关性建模故障诊断方法具有更高的故障分辨率，其输出的故障原因置信度可作为复杂系统维护保障的有力参考。贝叶斯网络建模较D矩阵建模复杂，需要大量的统计数据支撑，数据获取和处理难度较大，实际应用时还要考虑模型修正、概率参数随时间的动态变化等复杂因素。目前贝叶斯网络在装备故障诊断领域的实际应用案例比较少，相关方法、模型和算法有待优化、完善，还需投入大量的研究。
