矿井通风阻变型故障复合特征无监督机器学习模型

刘 剑，尹昌胜，黄 德，刘 丽，邓立军

(1.辽宁工程技术大学 安全科学与工程学院，辽宁 葫芦岛 125105; 2.辽宁工程技术大学 矿山热动力灾害与防治教育部重点实验室，辽宁 葫芦岛 125105)

矿井通风系统在生产中无法避免风门开关或破损、巷道冒落或变形、底臌或片帮等现象，这种现象会引起系统等效风阻发生持久或永久性变化，称为矿井通风系统阻变型故障[1]。故障发生地点所对应矿井通风网络分支所在的位置称为故障位置，对应故障位置发生等效风阻变化的大小称为故障量。矿井通风系统阻变型故障诊断的目的为准确识别阻变型故障位置及其故障量。

在工程应用中，针对不同领域有许多不同故障诊断方法，如支持向量机，贝叶斯，遗传算法，神经网络，模糊逻辑学，深度学习等[2-11]。目前，在水管网、 管网、电网及天然气管网等方面有较多的相关研究。有相关学者利用贝叶斯分类器对水管网的泄漏位置进行了诊断研究[10];文献[12]针对城市给水管网故障诊断，提出了实时故障诊断方法，利用人工神经网络技术诊断故障位置、故障程度和故障影响范围;文献[13]针对电网中的信息畸变、保护断路器不正确动作等异常信息，提出了通过训练最大熵隐马尔科夫模型对电网故障数据所隐藏的异常模式进行挖掘，从而实现电网的故障诊断;文献[14]利用随机矩阵谱理论，建立成品油管道管网模型，提出集中式管道管网故障诊断模糊决策方法和基于图论模型的管道管网故障诊断方法，用于成品油管道管网的故障诊断;文献[15]利用支持向量机对水管网局部堵塞频率进行了故障诊断研究;文献[16]利用变分模式分解和模糊相关分类方法对天然气管网的小泄漏进行诊断了研究;文献[17]利用贝叶斯网络模型针对天然气管网失效概率进行了研究;文献[18]利用粒子群优化支持向量机参数，基于特征熵对管道泄漏进行了定位研究;文献[19]提出了一种基于信度规则库模型的无线传感器网络故障诊断方法，并利用协方差矩阵自适应进化策略(CMA-ES)对无线传感器网络故障模型的初始参数进行优化。目前，对于矿井通风系统故障诊断方面的研究包括:凌六一等[20]针对矿井通风机机械故障，利用收集的故障样本，提出了结合BP神经网络和D-S证据理论对各证据相结合的矿井通风机机械故障诊断的方法;黄雷[21]针对矿井局部通风故障诊断，收集了矿井历史局部故障实例，基于Web技术开发了矿井局部通风故障诊断专家系统;龚晓燕等[22]收集了矿井局部通风系统故障样本，并利用遗传算法和神经网络将故障进行分类，对局部通风系统的故障类别和位置进行了诊断研究;赵丹等[23-24]利用专家系统理论技术及矿井通风监控系统建立了矿井通风智能诊断专家信息库，用以解决通风系统故障诊断存在主观判断和误判断的问题，并确定出矿井通风系统流体网络中的分支风量与分支风阻之间的数学关系，从而通过监测到的异常风速值间接求出故障分支的位置;笔者等[1]定义了矿井通风系统阻变型故障诊断，以各矿井通风流体网络每条分支的风量作为特征，以矿井通风仿真系统[25-26]仿真模拟的阻变型故障实例作为样本，利用支持向量机分类方法对阻变型故障进行了分类和回归研究;胡青伟[27]构建了矿井通风系统异常诊断指标体系，并实现对矿井通风系统的异常状态诊断，建立了矿井通风系统异常影响范围库，利用矿井通风在线监测数据和BP神经网络方法实现矿井通风系统异常原因识别的技术。

目前在矿井通风系统故障诊断方面的研究主要集中在机械故障、矿井通风系统稳定性及异常诊断等方面，而对矿井通风系统阻变型故障诊断方面的研究较少，且这些故障诊断方法中难以避免故障训练样本的收集。因此，在矿井通风系统阻变型故障诊断方法中存在几个亟待解决的问题:① 需要建立专家系统或收集故障样本，在诊断过程中离不开样本的参与;② 故障位置诊断和故障量诊断无法同时进行，在进行故障诊断之前需要分别建立故障位置诊断和故障量诊断模型;③ 难以实现矿井通风系统阻变型故障实时在线诊断。基于此，笔者基于CMA-ES进化算法[28]，提出一种无监督学习的故障诊断方法，建立分类与回归一体化的数学模型，实现矿井通风系统阻变型无监督故障诊断。此方法在诊断过程中无需故障样本参与训练，避免因样本数据噪声产生的误差或出现欠拟合及过拟合的现象，同时可达到矿井通风系统阻变型故障实时诊断的目的。

1 故障诊断模型及求解方法

1.1 问题描述

在矿井通风系统阻变型故障诊断对应的通风系统流体网络中，将故障位置对应的分支称为故障分支;将为故障诊断提供数据的分支称为观测点。设V为矿井通风系统流体网络中m个节点的集合，E为矿井通风系统流体网络中n条分支的集合，则矿井通风系统流体网络为G=(V,E)，其中V={v1,v2,v3,…,vm}，E={e1,e2,e3,…,en}。设矿井通风系统中未发生阻变型故障的通风系统流体网络观测点的原始风阻集合对应的向量为R0;发生阻变型故障后的观测点风量集合对应的向量为Q′;风压集合对应的向量为H′;矿井通风系统流体网络中故障分支编号为ex;故障量大小为Δrx，N·s2/m8，其中x∈(1,2,3,…,n)，表示故障位置对应分支编号索引。矿井通风系统阻变型故障诊断的目标为：在已知正常情况下，各观测点的风阻R0=(r01,r02,r03,…,r0n)及故障时所监测到的观测点风量Q′=(q′1,q′2,q′3,…,q′n)和风压H′=(h′1,h′2,h′3,…,h′n)，求出x及Δrx。其中，风阻单位为N·s2/m8，风量单位为m2/s，风压的单位为Pa。

对于任何正常运行或出现阻变型故障的矿井通风系统，均满足:

(1)

(2)

由式(1)可知，任何分支的等效风阻发生变化会引起整个系统各分支对应的风量及风压发生变化。反之，当监控到系统中分支风量和风压发生变化时，系统对应的流体网络中至少有一条分支的等效风阻发生变化，且常见为一条分支的等效风阻发生变化[1]。

在笛卡尔坐标系中，假设p=(p1,p2,…,pz)，q=(q1,q2,…,qz)为欧几里得空间中的两点，则点p和点q之间的欧式距离表达式为

d(p,q)=d(q,p)=

(3)

根据欧式距离的定义，将矿井通风系统阻变型故障后的风量(风压)向量与预设等效风阻所求的风量(风压)向量作为欧几里得空间中的两点，则矿井通风系统阻变型故障诊断问题可转化为求故障前后风量(风压)向量的欧式距离问题。若所求得的欧式距离大，则说明预设的等效风阻与实际发生的故障位置和故障量不符合，反之若求得的欧式距离越小，则说明所预设的故障位置和故障量越接近实际阻变型故障发生的情况。

1.2 无监督学习故障诊断数学模型

矿井通风系统阻变型故障诊断实际为求解欧式距离问题，而距离小方可得到准确的故障位置和故障量，因此，矿井通风系统阻变型故障诊断模型为求解故障前后故障分支风量(风压)欧式距离最小的优化问题。流体网络中观测点的风量(风压)与等效风阻之间的关系符合平衡定律:

(q′1,q′2,…,q′x,…,q′nq,h′1,…,h′x,…;h′nh)=f(Rx)

(4)

Rx=(r01,r02,…,r0x+Δrx,…,r0n)

(5)

式中，nq∈[0,n]为监测风量的观测点数，nh∈[0,n]为监测风压的观测点数;r0i为编号i的分支原始风阻，N·s2/m8，i=1,2,3,…,n;q′x为故障后分支ex的风量，m2/s;h′x为故障后分支ex的风压，Pa。

设待求解的向量为x=(x,Δrx)，其中x为故障位置，属于分类问题;Δrx为故障量，N·s2/m8，属于回归问题，利用式(6)作为矿井通风系统阻变型故障诊断无监督学习模型的目标函数，求得的最优解既包括分类结果，也包括回归结果，实现矿井通风系统阻变型故障诊断的分类与回归一体化。

(6)

1.3 故障诊断模型的求解

矿井通风系统阻变型故障诊断为一个求最小欧式距离的优化问题，优化求解方法有很多，其中CMA-ES算法是一种最成功的无导数全局优化算法，被归为一种进化策略算法[29]。因此，笔者所提的阻变型故障诊断模型求解方法为CMA-ES，下面介绍求解的关键流程。

xi=mt+σtyi

(7)

其中，t为迭代次数；i=1,2,3,…,λ，yi=BDN(0,I);I为单位矩阵。

(2)求解目标函数值。利用式(6)和式(7)求解各个解对应的目标函数值f(xi:λ)，其中i=1,2,3,…,λ，并对其进行从小到大排序操作:

f(x1:λ)≤f(x2:λ)≤…≤f(xλ:λ)

(8)

(3)更新均值参数。设μ为被选择解的个数，ωi为重组的正权系数，则均值参数的更新如下:

(9)

其中，ω1≥…≥ωμ>0。

(4)更新协方差矩阵。构造进化路径公式为

(10)

(11)

(5)更新步长参数。根据式(12)和(13)对步长参数进行更新。

(12)

(13)

其中，s0=0;cσ为学习速率的累加步长控制参数;dσ为更新步长的阻尼参数。

(6)终止判断。判断是否达到设定的最大迭代代数，如果达到则终止，输出x，否则根据式(7)～(13)继续迭代计算，直至达到终止条件终止。

2 模拟仿真试验

为了验证所提的矿井通风系统阻变型故障诊断无监督学习模型的可行性和求解方法的有效性。本文使用文献[1]的矿井通风系统及其数据，其网络图如图1所示，利用网络解算迭代，模拟矿井通风系统阻变型故障，并使用CMA-ES对数学模型进行求解。

为了证明CMA-ES方法对矿井通风系统阻变型故障诊断的适用性和泛化性，阻变型故障模拟覆盖除风井分支外的整个矿井通风网络，即所有试验的故障分支的集和为矿井通风中除风井分支外的所有分支。在进行试验之前需要对每条故障分支的风阻上下限进行估算，以避免网络解算过程中出现无法找到风机工况点的现象。根据文献[1]所提供的分支编号及对应的风阻，通过对每条故障分支的风阻不断进行增值和降值操作，最终得到各故障分支的风阻变化上下限见表1，其中风井变化范围为[0,0]，x为故障分支编号索引，即故障位置，lx和ux分别为各故障分支的等效风阻变化下限和上限，N·s2/m8。

模拟仿真试验过程中，每次试验选择一条分支作为故障分支，并在故障分支对应的等效风阻变化范围内随机产生一个浮点数作为故障量。如当x=3，故障分支为e3，l3=-0.003 N·s2/m8，u3=3.0 N·s2/m8，在-0.003与3.000之间随机取浮点数0.349 6作为故障量，即:Δr3=0.349 6 N·s2/m8。利用r03=r03+Δr3=0.03+0.349 6=0.379 6 N·s2/m8替换原始风阻，并进行矿井通风网络解算可得到一组风量集合向量Q′和一组风压集合向量H′，模拟产生了故障位置和故障量对应的监控风量和监控风压。初始化n=100，及各故障分支对应的风阻变化上、下限，使用式(7)～(13)进行迭代求最优解，当达到终止条件时，将输出的x和Δr与模拟产生的故障位置和故障量进行比较，若输出的故障位置与模拟产生的故障位置一致，表明故障位置诊断结果正确，否则故障位置诊断失败;若输出的故障量与模拟产生的故障量误差较大，表明故障量诊断结果不可靠，反之，误差越小说明故障量诊断结果越可靠。

图1 试验矿井通风系统流体网络

表1 故障分支等效风阻变化范围

续 表

针对不同的nq和nh大小，分别进行多次模拟试验，以分析风量特征、风压特征及风量-风压复合特征下的矿井通风系统阻变型无监督学习故障诊断性能及故障诊断模型的工程应用适用性。因此，针对所有的分支监控数据，设计4组模拟试验:①nq=100，nh=0;②nq=0，nh=100;③nq=100，nh=100;④ 随机给定nq和nh，并随机选择风量和风压监测的位置。每组试验过程中，分别以每条分支作为故障分支，模拟产生故障量，进行最优求解，并以故障位置诊断准确度和故障量误差大小作为故障诊断的评价指标。

3 结果分析

根据表1给定的故障分支编号及其对应的等效风阻变化下限和上限，除风井对应的分支外，每次不重复的以其中一条分支作为故障分支进行一次故障模拟，模拟得到的故障分支和故障量见表1，其中故障量Δrx的单位为N·s2/m8。所有故障分支的模拟故障量均在设定的等效风阻变化上下限范围内，符合试验设定条件。

根据表1提供的故障位置和故障量，进行矿井通风网络解算得到发生阻变型故障情况下的风量(风压)，利用所提出的阻变型故障诊断模型进行故障诊断，得到的故障位置和故障量诊断结果如图2所示。图中左侧为故障位置诊断结果与故障位置模拟结果的对比，右侧为故障量诊断结果与故障量模拟结果对比。其中在故障位置结果图中，若“模拟位置”与“预测位置”重合在一个坐标点上，说明故障位置的诊断结果与模拟结果一致，否则认为故障位置诊断失败，这种情况下的故障量诊断也被认定为失败;剔除故障量诊断失败的结果，得到图中右侧的故障量诊断结果，横坐标为模拟的故障量，纵坐标为预测的故障量，黑色的圆点离红色虚线的垂直距离越小，说明故障量诊断的结果越好，故障量诊断误差越小。如图3所示为故障量成功诊断情况下的误差分析，其中横坐标为试验序号，左侧纵坐标对应诊断和模拟故障量的相对误差，右侧纵坐标对应诊断和模拟故障量的绝对误差。

根据图2，3可得到，以风量为故障诊断的特征时(nq=100，nh=0)，故障位置和故障量诊断性能均不如以风压作为故障诊断特征(nq=0，nh=100)进行故障诊断得出的诊断性能好，而以风量-风压同时作为故障诊断特征(nq=100，nh=100)时，故障诊断的性能最佳，在图中表现为“模拟位置”和“预测位置”坐标重合点数更多，故障量诊断误差更小;相对误差呈现为上下波动少，与发生故障的故障量大小无关，故障诊断的绝对误差较其他更为平稳;相对误差整体比分别以风量或风压单独作为故障诊断特征时小，且相对误差在±5%内的试验次数最多。

通过对各故障分支发生故障进行统计分析得出表2，其中EMA表示平均绝对误差，EMS表示均方误差。① 当风量作为故障诊断的特征，并以所有分支作为观测点时，即nq=100，nh=0。故障位置诊断的准确率为86.71%，故障量诊断的平均绝对误差为0.432 3 N·s2/m8，均方误差为4.439 9 N2·s4/m16，R2为0.977 5。② 当风压作为故障诊断的特征，并以所有分支作为观测点时，即nq=0，nh=100。故障位置诊断的准确率为94.74%，故障量诊断的平均绝对误差为0.224 7 N·s2/m8，均方误差为1.052 7 N2·s4/m16，R2为0.995 1。③ 当风量-风压同时作为故障诊断的特征，并以所有分支作为观测点时，即nq=100，nh=100。故障位置诊断的准确率为96.55%，故障量诊断的平均绝对误差为0.118 4 N·s2/m8，均方误差为0.158 7 N2·s4/m16，R2为0.999 1。

图2 故障位置和故障量诊断结果

图3 故障量诊断误差

不同观测点比例下的故障位置诊断准确率如图4所示，其中横坐标表示在进行故障诊断时，选择的观测点数占所有分支总数的比例，纵坐标表示不同观测点比例下的故障位置诊断准确率。由故障诊断结果可知，即使观测点比例小于50%，故障位置诊断的准确率均超过80%，高于文献[1]中使用支持向量机及所有观测点进行故障位置诊断的准确率(78.11%)。当选择部分观测点进行故障位置诊断时，也可达到选择所有观测点进行故障位置诊断的准确率，且观测点比例大小与故障诊断的准确率无直接影响关系。整体故障位置诊断结果呈现出的趋势为:风量-风压混合特征下的故障诊断优于风压单一特征的故障诊断，而风量单一特征下的故障诊断性能相对不佳，但是其具有一定的稳定性，即使20%以下的观测点也可达到nq=100下的故障位置诊断准确率。

表2 试验结果统计

图4 不同观测点比例下的故障位置诊断结果

上述试验结果证明矿井通风系统阻变型故障位置诊断和故障量预测，风量-风压混合特征下的故障诊断试验结果均最优，全部分支的风量和风压监测下，故障位置准确率达到96.55%，故障量诊断的平均绝对误差为0.118 4 N·s2/m8，均方误差为0.158 7 N2·s4/m16，R2为0.999 1，而低于50%观测点比例下故障位置诊断准确率也可达到甚至超过96.55%。选择不同数量的风量分支和风压分支进行故障位置和故障量模拟和故障诊断试验，其结果见表3。当观测点比例超过50%，低于65%时，故障位置诊断的准确率超过97%，故障量诊断R2值超过0.99，平均绝对误差低于0.32 N·s2/m8，均方误差低于1.89 N2·s4/m16。结果表明，本文所提的方法可用于矿井通风系统阻变型故障诊断。

表3 不同风量和风压观测点数故障诊断结果

课题组在金川集团二矿进行实际应用试验。通过打开关闭状态下的风门和停机运行状态下的风机机站方式进行故障模拟，利用金川二矿矿井通风系统获取的矿井通风网络图，分支数为1 006条，节点数为870个，经过节点映射网络等技术后，矿井通风网络图分支数为178条，节点数为106个。利用观测点位置优化计算得到观测点数为总分支数的10%，但在复杂的井下环境中，其中3个观测点不具备人工测试数据的条件，共得到15观测点。本次共进行了5个阻变型故障工业试验，其中4个风门打开试验位置分别为:850水平材料车场风门、850水平通主斜坡道风门、1050水平措施道与回风道联络巷风门、1350水平充填回风道风门;1个风机停机试验位置为:850水平西翼风机站。5个故障点的故障诊断准确率均达100%。

综上述，本文所提出的无监督学习阻变型故障诊断模型可使用风量或风压单一特征进行矿井通风系统阻变型故障诊断，要获得更高的故障诊断准确率，需使用风量-风压复合特征进行故障诊断，即使选择部分观测点对应的巷道进行风量(风压)监控也具有较高的故障诊断准确率。

4 结 论

(1)使用CMA-ES对建立的分类与回归一体化的矿井通风系统阻变型故障诊断数学模型进行求解，可实现无需样本参与训练的无监督学习故障诊断，进而实现矿井通风系统阻变型故障的实时诊断。

(2)风量-风压复合特征比风量或风压单一特征下的矿井通风系统阻变型故障诊断可获得更高准确率的故障位置诊断和更低误差的故障量诊断。在实际应用中，结合风速和压差传感器同时配备使用可实现更为精准的阻变型故障诊断。

(3)在本文试验条件下，使用部分观测点条件下的风量-风压复合特征矿井通风阻变型故障诊断，故障位置诊断的准确率最高可达到99.29%，故障量诊断的R2值超过0.99，其故障诊断的性能优于风量或风压单一特征下的阻变型故障诊断。

(4)观测点比例大小与故障位置诊断准确率及故障量诊断误差无直接影响关系，但是与观测点的合理选配(传感器布置位置优化选定)有关，合理的观测点选配方案可实现在少量观测点数的情况下达到100%准确率的效果，此问题另文介绍。