基于角联子网的风量反演风阻病态改良算法

2019-05-08李雨成李俊桥邓存宝刘蓉蒸

煤炭学报 2019年4期

李雨成，李俊桥，邓存宝，刘蓉蒸

(1.太原理工大学安全与应急管理工程学院，山西太原 030024; 2.辽宁工程技术大学安全科学与工程学院，辽宁阜新 123000)

摩擦风阻是通风系统重要参数，是网络解算的基础，由于井下大部分故障均可归结于风阻改变，其变化也能够反映某些巷道是否发生故障[1-2]。目前，风阻主要通过人工测定的方式获取，测定过程不仅费时费力，也存在诸多弊端。一方面，限于测试条件，某些巷道风阻测不准，甚至无法测定[3-4]。例如:井下巷道风流存在湍流脉动现象[5]，读数时不可避免地产生观测误差;测定人员技术不熟练，也会影响测试准确度;又如，回风井为立井或预留回风小井的矿井，由于风硐无检查门或安全问题人员无法测定。另一方面，井下采掘是一个动态变化的过程，风门开启、矿车运行、提升设备升降、巷道冒落等都会引起风阻的变化，难以做到风阻的实时更新。相比于人工测试，使用监测监控数据(风量、风压)对风阻进行反演，不仅能够将技术人员从繁重的阻力测定任务中解放出来，还能实现风阻的实时更新，更能提高测试的安全性与准确度，因此是矿井风阻参数智能采集技术的发展方向[6-7]。

近年来，许多学者对测风求阻算法进行了研究。1991年，刘泽功[8]首次提出利用通风系统调风和阻力测定计算通风风阻，并给出了测算步骤与计算实例;2004年，周丽红等[9]提出使用风量数据能够计算出等同于独立回路个数的风阻，降低了风阻参数获取的难度;2012年，司俊鸿等[10]提出使用Tikhonov正则化法对测风求阻病态模型进行修正，可以减轻模型的病态程度;2015年，李雨成等[11]提出已知两组风量及部分节点压能数据，能够反演未知巷道风阻。至此，前人对测风求阻算法的研究已经十分成熟，但应用于大规模通风网络时，病态问题仍然没有得到彻底的解决，不能满足精度要求，导致无法实际应用。基于此，笔者引入文献[11]中的风量反演风阻算法，针对该算法的病态性质，提出将大规模通风网络划分为若干角联子网分别求解风阻，解决了算法的病态问题。同时给出了风压传感器优化布置算法，结合井下监测监控系统，能够实时反演风阻，不仅使矿井风阻参数智能采集成为了可能，更是未来实现风阻异常分支智能诊断的基础理论。

1 风量反演风阻算法原理

已知巷道风阻求解风量具有唯一解，而已知巷道风量求解风阻则属于多解问题，为了使风量反演风阻具有唯一解，可以限定部分节点压能。基于节点压能的风量反演风阻算法是已知两组巷道风量及部分节点压能，通过构建平衡方程组求解未知节点压能，根据通风阻力定律就能实现反演风阻。

对于通风系统任一巷道，能够得到平衡方程

hi-hj+Gij=rijQijQij

(1)

式中，hi(hj)为巷道节点静压能，Pa;Gij为巷道节点间位压差，Pa;rij为巷道风阻，N·s2/m8;Qij为巷道风量，m3/s。

通过调节某些构筑物的方式，获取两组风量。对于风阻不发生改变的巷道，将风量调节前后的平衡方程作比，消去不变量rij，得到风量反演风阻计算式(2)，解出未知节点压能，就能根据通风阻力定律实现风阻反演。

(2)

理论上，对于m条分支，n个节点的矿井通风系统，设已知压能节点个数为y，风阻发生改变巷道(一般为获取两组风量时调节的构筑物巷道)个数为z。为使方程定解，方程个数应不少于未知数个数，即

(3)

已知压能节点个数y满足式(3)时，风量反演风阻算法能够得到唯一解。但在算法的实际应用中，笔者发现该算法具有明显的病态特征，即使已知节点压能个数满足条件，也得不到正确的结果，计算风阻与实际风阻常常相去甚远。

2 算法病态性质分析

2.1 算法病态判别方法

对于线性方程组Ax=b，若矩阵A或向量b有小的扰动，引起解向量x很大的误差，即解向量x对A与b高度敏感，称方程组Ax=b为病态方程组[12]。算法实现过程中，数据扰动是不可避免的，因此病态问题的数值解是极其不稳定的。病态的界限比较模糊，条件数是衡量方程组病态程度的指标，条件数越大，病态越严重，解的误差也就越大。为了判断风量反演风阻算法的误差是否能够接受，需要设置条件数阈值，小于该阈值时，认为计算结果符合精度。精度要求高时，阈值需要设置的小一些，精度要求不高则可以设置的大一些。本文设置条件数的数量级小于103时，计算结果是非病态的。条件数计算方法见式(4)。

Cond(A)=‖A‖·‖A-1‖

(4)

式中，Cond(A)为矩阵A的条件数，无量纲;‖A‖(‖A-1‖)为矩阵A(A-1)的范数，无量纲。

2.2 算法病态原因分析

想要改良风量反演风阻算法病态问题，需要了解算法病态的原因。笔者经过大量的尝试与分析，总结出两点主要原因:

(1)系数矩阵稀疏性对算法病态的影响。

对于大规模矿井通风系统，分支通常多达数百条。在构建反演风阻方程时，任意风阻不发生改变的、包含未知压能节点的分支都能作为方程组的一个等式;未知节点压能构成了方程组的未知数，而同一条分支只能关联两个节点，即最多只能有4个未知节点压能。这说明无论系统规模多大，系数矩阵每一行至多只能有4个系数是非零的。这导致随着系统规模的增大，系数矩阵的稀疏性也将随之增大，而高维矩阵中元素的不确定性往往会导致矩阵的病态。可以从降低系统规模的角度，降低矩阵的稀疏性，进而改善算法的病态程度。

(2)通风系统拓扑结构对算法病态的影响。

为了不影响井下实际生产，调节构筑物的位置和数量会受到限制，导致获取两组风量的手段有局限性。对于不存在调节构筑物的并联巷道，其风量是大致成相同比例变化的，其风比平方系数Kij非常接近。这会使方程组的行向量过于线性相关，系数矩阵接近奇异矩阵，最终导致测风求阻算法的病态。由于系统整体拓扑结构已经固定，为了方便求解风阻而调整系统拓扑结构显然是不现实的。

3 风量反演风阻病态改良及风压传感器优化布置算法

目前，求解病态方程的手段有很多种，如正则化解法、奇异值分解法[13-14]等。但这些方法都是通过改善条件数的方法降低算法的病态程度，不能从根本上解决算法的病态问题，实际效果并不理想。应用于本文的测风求阻算法时，仍然难以满足求解精度要求。因此，笔者考虑从算法模型的建立上降低算法的病态程度。

3.1 基于角联子网的风量反演风阻病态改良算法

角联结构普遍存在于通风系统中，具有风流不稳定的特点[15-17]。图1为简单角联结构，记作D={e1(v1,v2),e2(v1,v3),e3(v2,v3),e4(v2,v4),e5(v3,v4)}。

图1 简单角联结构Fig.1 A simple diagonal structure

将大型矿井通风系统划分为若干小规模的、算法模型良态的角联子网，进而对各角联子网分别求解风阻。这样不仅降低了系统规模，系数矩阵的稀疏性大大降低，而且角联子网构建的方程组行向量不具有相关性，因此是解决算法病态问题的有效手段。算法流程如下。

步骤1:计算通风网络中所有的角联结构，优先选择分支个数少的角联子网作为待求子网，选择的角联子网D1,D2,…,Dn应尽可能覆盖通风网络除进、回风井外的所有分支(进、回风井由于其拓扑性质，不会包含在角联结构中)，即

(5)

式中，EDi为角联结构Di的分支集合;∪为计算各集合的并集;E为通风网络分支集合;Ef为进风井分支集合;Et为回风井分支集合。

步骤2:赋予各节点权重，根据节点权重计算风压传感器最优布置方案，在3.2节中将详细介绍传感器优化布置算法。布置风压传感器的节点作为已知压能节点参与后续风量反演风阻计算。

步骤3:使用风量反演风阻算法计算各角联子网分支风阻。对于各子网，当参与计算的分支数等于未知压能节点数的2倍时，具有唯一解;分支数大于未知压能节点数的2倍时，具有最小二乘解;分支数小于未知压能节点数的2倍时，不具有唯一解。

步骤4:检查、筛选计算结果。① 舍去条件数≥103的角联子网的病态计算结果。② 检查方程系数矩阵任一列的非零元个数。若存在某一列非零元个数<2，一般存在于包含风阻改变分支的角联结构，由于该分支不参与运算，导致该列中非零元对应的节点压能是无解的，即使条件数符合要求，计算结果也是不可靠的。③ 一条分支存在于多个角联子网时，选择条件数较小的计算结果。

理论上，只要风量值和压能值是准确的，计算风阻可以达到很高的精度。但实际上传感器观测数据可能不准确，导致误差增大。由于有些分支存在于多个角联子网中，可以利用该特性进行检查，若某一传感器故障，会出现各角联子网计算结果不一致的现象。因此，该算法具有一定的自我纠错的能力。另外，为了尽可能减小误差，对传感器数据进行预处理也是必要的，如节点风量平衡修正、数据概率平均修正等[5]。

3.2 基于贪心策略的风压传感器优化布置算法

不合理的风压传感器布置方案既不经济、也无法满足3.1节中风量反演风阻病态改良算法的要求。因此，风压传感器的优化布置是值得研究的。为避免风量反演风阻算法病态，从控制系统规模的角度考虑，需要控制方程的系数矩阵的列数不超过4，即每一个角联子网中的未覆盖传感器节点不能超过两个。为了使风压传感器的布置数量最少，而且满足各角联子网的计算要求，笔者给出了基于贪心策略的风压传感器优化布置算法。首先设置节点权重为节点在各角联子网中出现的次数，节点权重越大，说明有更多的子网包含这个节点。权重可以根据矿井节点实际情况进行调整，能够使风压传感器布置更加合理。算法思路为:逐个选择节点，每一步选择节点时，优先选择权重最大的节点;权重相同时，优先选择能够使布点均匀分布于各角联子网的节点。图2为基于贪心策略的风压传感器优化布置算法程序框图。

图2 传感器最优布置算法Fig.2 Algorithm of sensor optimized layout

步骤1:初始化角联子网节点集合VD1,VD2,…,VDn;角联子网中未覆盖传感器节点最多个数N=2;待选节点集合V，备选节点集合V′，布点集合V″，节点权重Weight。

步骤2:在V中寻找最大权重节点集合va。

步骤3:判断va是否为空集。若不是空集，转步骤4;若是空集，转步骤5。

步骤4:计算va中使varVDj-V″-va[t]最小的节点，记作va[t]，将其加入布点集合V′，加入该节点后，可能存在布点个数满足要求的子网，将这些子网中未覆盖传感器的节点加入到备选节点集合V′。其中，var为计算方差，varVDj-V″-va[t] 为衡量布点是否分布均匀的指标，方差越小，说明布点越均匀。在节点权重相同的节点中优先选择使布点均匀分布于角联子网的节点。

步骤5:判断是否存在传感器个数未达到要求的子网，其节点集合记作VDk。若存在，转步骤6;若不存在，转步骤7。

步骤6:在V′与VDk的公共节点中寻找最大权重节点集合，若计算结果包括多个节点，优选使varVDj-V″-va[t] 最小的节点，将其加入布点集合V″，并从备选节点V′集合中删除。

步骤7:返回V″，程序结束。

4 实例分析

如图3所示的通风系统包括23条分支、16个节点、6个构筑物、1个进风井、1个回风井。通风方法为抽出式，通风机特性曲线方程为H=-0.052 3Q2+8.560 8Q+3 061.6。为便于计算，假设各节点标高相同，即不考虑位压能。调节分支e3，e19上的构筑物获得两组风量，分支e3风阻由30.058 640变为0.058 640，分支e19风阻由18.028 300变为0.028 300。分支实际风阻及调节构筑物前后两组风量见表1，节点压能见表2。限于条件，两组风量及压能均是根据实际风阻进行网络解算得到，虽然不是传感器数据，但同样能够对算法进行验证。

应用Matlab软件编写了计算程序calculResis.m，将两组风量及部分节点压能作为已知条件，分别使用整体直接计算风阻和划分子网计算风阻两种算法求解。通过对比计算风阻与实际风阻，验证直接求解风阻的病态特征及划分角联子网求解风阻的可行性。

4.1 整体直接求解

根据式(3)，计算得到满足求解条件的已知压能节点个数为

图3 通风网络Fig.3 Ventilation network graph

表1 分支实际风阻及风量Table 1 Real resistance and air quantity of branches

表2 节点压能Table 2 Node pressure

y≥0.5(2n-m+z)=0.5×(2×16-23+2)=5.5

4.2 划分角联子网求解

(1)角联子网划分。

计算通风网络中所有的角联结构，从中选择了8个分支数最少的角联结构作为待求子网，其中子网分支数最多不超过9条，覆盖了通风网络除进、回风井外的所有分支。选择的各角联子网如下:

D1={e15(v9,v11),e16(v9,v12),e18(v11,v13),

e19(v12,v13),e20(v12,v14),e21(v13,v15),e22(v14,v15)}

D2={e5(v3,v4),e6(v3,v8),e7(v4,v5),e9(v5,v8),

e10(v5,v9),e13(v8,v10),e15(v9,v11),e17(v10,v11)}

D3={e5(v3,v4),e6(v3,v8),e7(v4,v5),e9(v5,v8),

e10(v5,v9),e14(v8,v14),e16(v9,v12),e20(v12,v14)}

D4={e2(v2,v3),e3(v2,v6),e5(v3,v4),e6(v3,v8),

e7(v4,v5),e8(v4,v6),e9(v5,v8)}

D5={e2(v2,v3),e4(v2,v7),e5(v3,v4),e6(v3,v8),

e7(v4,v5),e8(v4,v6),e9(v5,v8),e11(v6,v7)}

D6={e9(v5,v8),e10(v5,v9),e13(v8,v10),e15(v9,v11),

e16(v9,v12),e17(v10,v11),e18(v11,v13),e19(v12,v13)}

D7={e9(v5,v8),e10(v5,v9),e14(v8,v14),e16(v9,v12),

e19(v12,v13),e20(v12,v14),e21(v13,v15),e22(v14,v15)}

D8={e3(v2,v6),e4(v4,v7),e7(v4,v5),e8(v6,v4),

e10(v5,v9),e11(v7,v6),e12(v7,v10),e15(v9,v11),

e17(v10,v11)}

(2)风压传感器优化布置方案。

计算各子网节点集合:

VD1={v9,v11,v12,v13,v14,v15}

VD2={v3,v4,v5,v8,v9,v10,v11}

VD3={v3,v4,v5,v8,v9,v12,v14}

VD4={v2,v3,v4,v5,v6,v8}

VD5={v2,v3,v4,v5,v6,v7,v8}

VD6={v5,v8,v9,v10,v11,v12,v13}

VD7={v5,v8,v9,v12,v13,v14,v15}

VD8={v2,v4,v5,v6,v7,v9,v10,v11}

不考虑节点在通风系统中的实际情况，仅根据节点在角联子网中出现的次数，设定节点权重。v2～v15节点权重分别为w2=3，w3=4，w4=5，w5=7，w6=3，w7=2，w8=6，w9=6，w10=3，w11=4，w12=4，w13=3，w14=3，w15=2。使用基于贪心策略的风压传感器布置算法计算传感器布置方案，选择节点v2，v4，v5，v7，v8，v9，v11，v12，v15为布置风压传感器节点，计算过程见表3，序号0代表初始化。传感器布置方案完成后，传感器节点的压能作为已知条件参与风阻计算。

表3 传感器布置方案计算过程Table 3 Computing process of sensor layout

(3)反演风阻。

对各角联子网进行风阻求解。计算完成后，检查计算结果的可靠性。D4，D7条件数太大，D6虽然条件数较小，但方程系数矩阵结构不合理，因此舍去这3组计算结果。最终整理计算结果见表4，可以直观地看到计算风阻与实际风阻非常接近，误差控制在10-4以内，可以忽略不计。计算结果不包括进、回风井分支，若要实现进、回风井分支反演，可以在进、回风井口节点v1，v16增加风压传感器，就能实现所有分支风阻的反演。