刘兆伦1,2张春兰2武尤2王海羽2刘彬1,2
1.燕山大学河北省特种光纤与光纤传感实验室，秦皇岛，0660042.燕山大学信息科学与工程学院，秦皇岛，066004
关键词：贝叶斯增量学习；结构学习；篦冷机；水泥熟料换热；二次风温；故障诊断
针对上述现状，本文将专家知识与数据学习相结合得到原结构和原参数，构造WTUN(whethertoupdatenetwork)函数用于判断新数据集到达时原结构是否改变，若是，则定义Affect函数用于衡量数据对结构中各节点的平均影响程度，利用爬山算子对影响较大节点的马尔可夫毯结构进行修正得到候选结构，利用改进的评分函数选择评分最大的结构作为最优结构，无论结构是否改变均利用EM算法更新参数，得到一种贝叶斯增量学习算法(incrementallearningalgorithm,ILA)。将该算法应用到篦冷机工艺参数建模中，利用模型分析参数间的影响，对原模型进行增量维护，最后对二次风温进行故障诊断，找出导致故障发生的原因并采取措施，可提升热回收效率从而达到节约能源的目的。
数据学习要求获得大量数据样本，且所得结构不准确，专家知识学习主观性较强，易产生较大误差，将两种方法相结合会相对准确且能提高效率[10]，本文结合数据学习与专家知识获得原贝叶斯结构G，利用最大似然估计(maximumlikelihoodestimation,MLE)算法学习得到原参数θijk。
(1)
(2)
其中，P(D/G)表示似然概率，即给定结构情况下关于数据的后验概率；qi是变量Xi父节点取值组合的数目；ri是变量Xi取值的数目；Nijk是Xi取k值、Xi的父节点取j值时的样本数目，也称为充分统计因子。由式(2)可以判断数据与结构的拟合程度，其实质是对每个样本在结构中有相互影响关系的变量条件概率的求和[11]，依据其思想，构建WTUN函数：
WTUN(D|G)=
(3)
依据马尔可夫毯原则[12]对式(3)进行化简。在BN中节点xi的马尔可夫毯结构包括它的父节点、子节点及子节点的父节点，当节点xi的马尔可夫毯给定时，在该网络中xi与其他变量条件独立，因此，网络中的条件概率可化简为
P(xi|x1,x2,…，xi-1,xi+1,…,xn)=
P(xi|MB(xi))
(4)
依据乘法公式，可将式(4)变换为
P(xi|MB(xi))=P(xi,MB(xi))/P(MB(xi))
(5)
由于式(5)中P(MB(xi))不包含xi，即无论xi取何值，P(MB(xi))的取值相同，视为常量，根据随机变量形式链式规则，将式(5)化简为
(6)
根据马尔可夫毯原则将式(6)中的分成3部分，由于其他节点的条件概率项与节点xi无关，可视为常量，故式(6)可化简为
P(xi|MB(xi))=bP(xi|π(xi))·
(7)
式中，P(xi|π(xi))为xi的条件概率；为xi子节点的条件概率。
将式(7)代入式(6)，再代入式(4)可得
(8)
式中，π(xi)为节点xi的父节点；xj为贝叶斯结构中节点xi的子节点；π(xj)为节点xi的子节点的父节点。
由式(8)可看出，WTUN函数是对结构中有相互影响关系变量的条件概率的求和，与BIC函数计算结果相比，只是某些概率项的系数有所增大，故可用于判断数据与原结构的拟合程度。对于给定的Do、Dn，以及正数阈值β∈[0,1]，若
(9)
则需更新结构，否则只更新参数。
若式(9)不成立，即原结构不需要调整，只需对BN进行参数学习，根据EM算法思想[13]，将原参数变为先验参数，将以下结果作为新的参数：
(10)
式中，为先验参数；为增量数据的充分统计因子。
若式(9)成立，则需要进行结构调整，在新的结构上更新参数。构造Affect函数，计算当前数据对各节点的平均影响程度，给定阈值δ，存储并记录Affect值大于阈值δ的节点集S，通过爬山算子S中各节点的马尔可夫范围进行修正，得到候选结构集Gq，利用改进评分函数对Gq进行评分，选取评分最高的结构作为最优结构Ggood，并利用式(10)更新参数，完成BN的增量维护。
已知P(Dn|G)可以体现新数据与原结构的拟合程度，新数据对G的所有影响由P(Dn|G)的变化相对于参数θijk的变化组成[14]，则当前数据对原结构中各节点的平均影响表示为
(11)
式中，A为新数据对原结构每个节点的平均影响；B为原数据对原结构每个节点的平均影响。
对A进行化简，有
(12)
其中，dn为新数据集样本个数；P为结构对应的条件概率集，又已知
则式(12)可化简为
(13)
同理可得B，则Affect函数表达式为
(14)
θijk=P(Xi=k|π(Xi=j))
式中，为原数据集的充分统计因子；k为该节点条件概率表中概率的个数；α为新数据集所占比率。
利用式(14)计算每个节点的Affect值，给定阈值δ，存储并记录Affect值大于阈值δ的节点集S，利用增加弧、反转弧、删除弧三种操作集合op={op(1),op(2),…,op(k)}修正S中每个节点的马尔可夫毯结构，得到候选结构集：
Gq=op(j)(MB(Si))
(15)
式中，op(j)为3种操作之一；MB(Si)为节点集S中某个节点的马尔可夫毯结构。
对Gq中的结构进行评分：
(16)
式中，Nn、No分别为新旧数据的数量；λ为调整新旧数据学习速度及趋势的因素，如果新数据与原结构不匹配，则λ增大，否则λ减小；g为结构集Gq中的个体；pen(g)为结构g的惩罚函数，防止数据与结构的过度拟合。
选取Gq中评分最大的结构作为最优增量结构Ggood，即
(17)
得到Ggood后，利用式(10)更新参数，完成BN的增量学习过程。
ILA算法的流程图见图1。对ILA算法分步进行简要描述，算法步骤如下：
图1ILA算法流程图Fig.1TheflowchartofILAalgorithm
(1)输入初始数据集、批量式算法与专家知识结合获得原贝叶斯结构G，利用MLE算法学习得到原参数；
(2)检测到新数据，利用WTUN函数判断G是否改变，若是则转步骤(3)，否则利用EM算法更新参数，转步骤(6)；
(3)利用Affect函数找到数据对节点影响大的节点集S，并得到每个节点的马尔可夫毯结构的集合Gp；
(4)在原结构中利用op(j)对步骤(3)得到的Gp中的结构进行修改，得到候选结构集Gq；
(5)使用Score函数对步骤(4)中的每个候选结构进行评分，以评分最大的结构作为增量结构Ggood，利用EM算法更新参数；
(6)输出增量学习后的结构和参数。
有新数据的WTUN函数值越大，说明贝叶斯网络对新数据的拟合程度越好，随着新数据的增加，与新旧结构的WTUN值在不断变化，若与新结构的相适性越来越好，则可将WTUN函数用于判断数据与结构拟合程度的可行性。
利用经典Asia、Car及Alarm网络进行仿真验证，以Asia网络为例，其标准网络如图2所示。利用Netica软件的smulatecases功能生成8000组数据作为原数据集Do，修改概率分布表，在不产生环路的条件下，删除和增加一定数量的边后得到的网络如图3所示，用同样的方法再产生8000组数据作为新数据集Dn，部分数据结果见表1。
图2原Asia网络Fig.2TheoriginalAsianetwork
图3新Asia网络Fig.3ThenewAsianetwork
表1Asia网络部分数据举例
Tab.1SomedatasamplesinAsia
IDnumTuberculosisXRayTbOrCaCancerDyspneaBronchitisVisitAsiaSmoking1absentnormalfalseabsentabsentabsentno_visitnon_smoker2absentabnormalfalsepresentpresentpresentno_visitsmoker3presentabnormaltrueabsentpresentabsentno_visitsmoker4absentabnormaltruepresentpresentabsentno_visitsmoker5absentnormalfalseabsentabsentabsentvisitnon_smoker………………………
对图2所示的网络，分别加入0～100%的新数据，随着新数据的增加，分别计算含不同比例新数据下的WTUN(Dn|G)及WTUN(Dn|Gn)，对Car、Alarm与Asia进行相同操作，WTUN值的变化情况分别如图4和图5所示。
图4WTUN值随原结构变化情况Fig.4TheWTUNvaluevarieswiththeoriginalstructure
图5WTUN值随新结构变化情况Fig.5TheWTUNvaluevarieswiththenewstructure
由图4可知，随新数据的增加，与原结构的WTUN值越来越小，即数据与原结构的适应性越来越差。由图5可知，与新结构的WTUN值越来越大，即数据与新结构的拟合程度越来越好，而且变化比较平稳，从而验证了WTUN函数的有效性。
利用Asia网络验证Affect函数的正确性。同样利用图2所示网络作为原网络，生成2000组原数据，利用图3所示网络作为新网络，生成8000组新数据，即考虑10000组数据中含20%原数据和80%新数据的情况，说明α的值为0.8。利用式(14)计算数据对原网络各个节点的Affect值，并给定阈值δ为1，以节点Bronchitis为例说明式(14)的具体计算过程，该节点的原参数与新参数的学习结果分别见表2和表3。
当Bronchitis=present，Smoking=smoker时，节点Bronchitis的Affect值为
表2节点Bronchitis的原条件概率
Tab.2TheoriginalconditionalprobabilityofnodeBronchitis
Smoking=smokerSmoking=nonsmokerBronchitis=present0.60.4Bronchitis=absent0.30.7
表3节点Bronchitis的新条件概率
Tab.3ThenewconditionalprobabilityofnodeBronchitis
Smoking=smokerSmoking=nonsmokerBronchitis=present0.57240.4276Bronchitis=absent0.33870.6613
其中，497和156分别为Do与Dn中满足Bronchitis=present、Smoking=smoker的样本数量。同理可计算另外3种情况下的数值，这4个值求和的结果即节点Bronchitis的Affect值。因此，可得到每个节点的Affect值，结果见表4。
表4Asia网络中数据对各节点的影响情况
Tab.4TheinfluenceofdataoneachnodeinAsianetwork
节点及简称Affect值与阈值相比是否要修改与原结构对比是否判断正确VisitToAsia(v)1.862是否Smoking(S)39.467是是Tuberculosis(T)2.474是是LungCancer(L)2.5888是是Bronchitis(B)12.396是是TuberculosisorBronchitis(O)0.032否否XRay(X)0.24否是Dyspnea(D)185.7682是是
由表4结果可知，Affect函数判断数据对结构影响度的准确率为75%。为了排除数据的偶然性，在上述条件下进行20次试验取平均值，最终得到准确率为72.5%，由此可以验证Affect函数的正确性。
在Affect函数正确性验证实验中，由表4可知，当加入8000组新数据时，根据Affect值大小判断需要修正的节点分别为D、S、B、v、L、T，对这些节点马尔可夫毯结构用爬山算子进行修改，利用Score函数进行评分，得到评分最大的最优网络结构，如图6所示。
图6ILA算法得到的最优结构Fig.6Optimalstructureobtainedbyincrementalalgorithm
(18)
(19)
(20)
其中，*代表与ILA对比的算法，Log-less值的计算公式为代表数据集，PB(d)代表学习到贝叶斯网络的预测概率，Dtest为采样的数据集。
Tab.5ComparisonofthequalityoflearningresultsandtimeconsumptionbetweenILAandHCalgorithm
Tab.6ComparisonofthequalityoflearningresultsandtimeconsumptionbetweenILAandIHCalgorithm
根据Alarm网络生成的训练数据，分别读入k为500、1000个事例进行一次增量学习，分别比较ILA算法与HC算法、增量遗传算法(incrementalgeneticalgorithm，IGA)[7]的存储量，如图7所示。
图7三种算法的存储量对比Fig.7Thecomparisonofthestoragecapacityofthreealgorithms
由图7可以看出，批量HC算法的存储量随事例数的增大呈线性增长，这是因为批量HC算法每次要将新旧数据合在一起重新学习，IGA需要的存储量先上升达到一个顶峰后下降的过程，虽然后期存储量趋于常数，但前期存储量较大，而本文所提算法通过Affect函数判断出新数据对原结构哪些节点影响大，只对这些节点进行了局部修改，极其有效地提高了效率，减少了存储空间，存储开销优于HC与IGA算法。综上所述，ILA算法是有效的。
篦冷机是新型干法水泥生产线的关键设备，担负着出窑熟料的冷却、输送和热回收任务[15]，其工作状态直接影响到水泥厂的效益及水泥用户的利益，但其工艺参数较多且影响错综复杂，且故障具有很大不确定性，人工排查方式效率很低。由于BN强大的推理能力及方便的决策机制，故将ILA算法应用到篦冷机熟料换热系统，对诊断模型进行维护。首先根据篦冷机原理确定BN的变量，其次结合专家知识与数据学习建立篦冷机原始模型，然后利用ILA算法的增量机制更新初始模型，最后利用变量消元算法进行故障诊断。
篦冷机结构如图8所示，根据其工艺原理选取变量，其范围与状态分类的对应关系见表7。利用MATLAB对这8个参数的数据进行离散化处理，然后通过专家知识与SHC算法[16]结合的方法建立原始模型，利用EM算法对所建立的结构和原数据进行参数学习，得到原始模型，如图9所示。
图8篦冷机结构示意图Fig.8Schematicdiagramofgratecooler
表7变量范围与变量状态分类的对应关系
Tab.7Rangeofvariablesandcorrespondenceofvariablestateclassification
篦冷机速度v(次/min)小于12Low(2)12～14Norm(1)大于14High(3)一段篦下压力Gp(Pa)小于5400Low(2)5400～6000Norm(1)大于6000High(3)平衡风机转速R(r/min)小于980Low(2)980～1090Norm(1)大于1090High(3)生料下料量M(t/h)小于320Low(2)320～340Norm(1)大于340High(3)二次风温度St(℃)小于870Low(2)870～1070Norm(1)大于1070High(3)三次风温度Tt(℃)小于890Low(2)890～950Norm(1)大于950High(3)窑头负压Ph(Pa)小于-80Low(2)-80～-25Norm(1)大于-25High(3)窑主机电流I(A)小于430Low(2)430～530Norm(1)大于530High(3)
图9篦冷机原始模型Fig.9Originalmodelofgratecooler
在原始模型的基础上，检测到新数据时，利用ILA算法，加入1000组新数据和加入10000组新数据时的学习结果分别如图10、图11所示，实现了诊断模型的更新。
图10加入1000组新数据时增量学习结果Fig.10Incrementallearningresultswhenadding1000newsetsofdata
图11加入10000组新数据时增量学习结果Fig.11Incrementallearningresultswhenadding10000newsetsofdata
由图10和图11可以看出，1000组新数据不足以引起结构的改变，只是参数进行了更新，而加入10000组新数据则引起了结构的变化，同时参数也发生改变。根据初始模型，利用ILA、IHC算法进行增量学习，利用SHC算法进行批量学习得到的诊断模型如图12～图14所示。由图12～图14可知：利用SHC得到的诊断模型偏差较大，如忽略了喂料量对二次风温的影响，错误添加了二室风温对二室篦下压力的影响。
图12利用ILA算法得到的诊断模型Fig.12DiagnosismodelbasedonILAalgorithm
图13利用IHC算法得到的诊断模型Fig.13DiagnosismodelbasedonIHCalgorithm
图14利用SHC算法得到的诊断模型Fig.14DiagnosismodelbasedonSHCalgorithm
图16三种算法得到篦冷机结构的存储量Fig.16Comparisonofthethreealgorithmsforthestoragecapacityofthestructureofthegratecooler
二次风温是篦冷机的重要参数，由篦下高压冷却风机与高温熟料进行热交换得到，二次风在窑头负压作用下，进入回转窑供窑内燃烧，二次风温偏低会使热回收效率变差，造成燃料煤的大量浪费，同时产生大量废气污染环境，提高二次风温可、节省燃煤，达到节约能源的目的。如监测到二次风温忽然偏低，由图12可知导致二次风温忽然偏低的原因为V、R、M及Gp，以10000组数据为例，利用变量消元法计算后验概率：
α(0.06480.11740.8178)
α(0.35240.21910.4285)
=
α0.24160.65450.1039
α0.32480.21590.4393
通过比较后验概率α可知，二次风温忽然降低是由篦速偏高导致的，应及时减小篦速，最大程度上减少燃煤的浪费与氮氧化物等污染废气的排放。从数据集中筛选出满足二次风温偏小征兆的数据1000、3000、5000、8000组进行测试性诊断实验，记录每一组数据对应的诊断结果，部分验证结果见表8。
表8ILA算法的部分验证结果
Tab.8PartialverificationresultsofILAalgorithm
序号二次风温篦速诊断结果1St=2R=3正确2St=2R=3正确3St=2R=1错误4St=2R=3正确
由诊断结果可得到各算法在每一组测试集的正确诊断故障数：
H=Dture/Dtest
(21)
其中，H为正确诊断率；Dture为诊断正确的数据数；Dtest为测试集大小。最终得到3种算法对比结果，如图17所示。由图17可知：在每组测试数据中，利用ILA对二次风温进行故障诊断的准确率要高于SHC算法以及IHC算法，能够基本满足实际生产中二次风温故障诊断的需求。
图17三种算法诊断准确率对比Fig.17Comparisonofdiagnosticaccuracyofthreealgorithms
参考文献：
[1]杨昌昊，竺长安，胡小建.基于贝叶斯网的复杂系统故障诊断方法[J].中国机械工程,2009,20(22):2726-2732.
YANGChanghao,ZHUChang’an,HUXiaojian.MethodsonComplexSystemFaultDiagnosisBasedonBayesianNetwork[J].ChinaMechanicalEngineering,2009,20(22):2726-2732.
[2]SAMETS,MINA,GRANGERE,etal.IncrementalLearningofPrivacy-preservingBayesianNetworks[J].Appl.SoftComput.,2013,13(8):3657-3667.
[3]刘彬,刘永记,刘浩然,等.基于贝叶斯改进结构算法的回转窑故障诊断模型研究[J].中国机械工程,2017,28(18):2143-2151.
LIUBin,LIUYongji,LIUHaoran,etal.ResearchonFaultDiagnosisModelofRotaryKilnBasedonImprovedBayesianStructureAlgorithm[J].ChinaMechanicalEngineering,2017,28(18):2143-2151.
[4]ZHUY,LIMINH,LUJ.BayesianNetworks-basedApproachforPowerSystemsFaultDiagnosis[J].IEEETransactionsonPowerDelivery,2006,21(2):634-639.
[5]田凤占,黄丽,于剑,等.包含隐变量的贝叶斯网络增量学习方法[J].电子学报,2005,33(11):1925-1928.
TIANFengzhan,HUANGLi,YUJian,etal.IncrementalLearningMethodforBayesianNetworkswithHiddenVariables[J].JournalofElectronics,2005,33(11):1925-1928.
[6]ALCOBEJR.AnIncrementalAlgorithmforTree-shapedBayesianNetworkLearning[C]//EureopeanConferenceonArtificialIntelligence.Lyon,2002:350-354.
[7]王飞,刘大有,王淞昕.基于遗传算法的Bayesian网结构增量学习的研究[J].计算机研究与发展,2005,42(9):1461-1466.
WANGFei,LIUDayou,WANGSongxin.ResearchonIncrementalLearningofBayesianNetsBasedonGeneticAlgorithm[J].ComputerResearchandDevelopment,2005,42(9):1461-1466.
[8]ZHUY,LIUD,CHENG,etal.MathematicalModelingforActiveandDynamicDiagnosisofCropDiseasesBasedonBayesianNetworksandIncrementalLearning[J].Mathematical&ComputerModelling,2013,58(3/4):514-523.
[9]YASINA,LERAYP.IncrementalBayesianNetworkStructureLearninginHighDimensionalDomains[C]//InternationalConferenceonModeling,SimulationandAppliedOptimization.Hammamet,2013:1-6.
[10]慕春棣,戴剑彬,叶俊.用于数据挖掘的贝叶斯网络[J].软件学报,2000,11(5):660-666.
MUChunli,DAIJianbin,YEJun.BayesianNetworksforDataMining[J].JournalofSoftware,2000,11(5):660-666.
[11]LISH,ZHANGJ,SUNBL,etal.AnIncrementalStructureLearningApproachforBayesianNetwork[C]//ControlandDecisionConference.Changsha,2014:4817-4822.
[12]MADDENMG.EvaluationofthePerformanceoftheMarkovBlanketBayesianClassifierAlgorithm[R].Galway:NationalUniversityofIreland,2002.
[13]王宝帅,杜兰,和华,等.基于复高斯模型的样本缺失窄带雷达信号重构算法[J].电子与信息学报,2015,37(5):1065-1070.
WANGBaoshuai,DULan,HEHua,etal.ReconstructionAlgorithmofNarrowbandRadarSignalswithMissingSamplesBasedonComplexGaussModel[J].JournalofElectronicsandInformation,2015,37(5):1065-1070.
[14]YUEK,FANGQ,WANGX,etal.AParallelandIncrementalApproachforData-intensiveLearningofBayesianNetworks[J].IEEETransactionsonCybernetics,2015,45(12):2890-2904.
[15]王美琪,刘彬,闻岩,等.篦冷机内水泥熟料温度的软测量[J].机械工程学报,2016,52(6):159-165.
WANGMeiqi,LIUBin,WENYan,etal.SoftMeasurementofCementClinkerTemperatureinGrateCooler[J].JournalofMechanicalEngineering,2016,52(6):159-165.
[16]刘浩然,吕晓贺,李轩,等.基于Bayesian改进算法的回转窑故障诊断模型研究[J].仪器仪表学报,2015,36(7):1554-1561.
LIUHaoran,LYUXiaohe,LIXuan,etal.AneffectiveInferenceAlgorithmforFaultDiagnosis[J].ControlandDecision,2015,30(11):2033-2040.
LIUZhaolun1,2ZHANGChunlan2WUYou2WANGHaiyu2LIUBin1,2
1.HebeiProvinceKeyLaboratoryofSpecialOpticalFiberandOpticalFiberSensing,YanshanUniversity，Qinhuangdao,Hebei,0660042.InformationScienceandEngineeringCollege，YanshanUniversity,Qinhuangdao,Hebei,066004
Keywords:Bayesianincrementallearning;structurelearning;gratecooler;clinkerheattransfer;secondaryairtemperature;faultdiagnosis
中图分类号：TH165.3
DOI：10.3969/j.issn.1004-132X.2019.10.005
开放科学(资源服务)标识码(OSID)：
收稿日期：2017-11-28
基金项目：国家自然科学基金资助项目(51641609)
(编辑陈勇)
作者简介：刘兆伦，男，1978年生，副教授。研究方向为智能算法及故障诊断。发表论文10余篇。E-mail:liuzhaolun@ysu.edu.cn。