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

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

摘要:为解决矿井通风系统阻力测定过程中部分巷道受限于现场测试条件而导致的风阻测不准、甚至无法测定的问题，研究了基于节点压能的风量反演风阻算法，提出使用两组风量数据及部分节点压能数据反演风阻。针对算法的病态性质，对算法病态的原因进行了分析，得出反演风阻方程系数矩阵稀疏性和通风系统拓扑结构对算法的影响是导致算法病态的两个主要原因。基于此，提出基于角联子网的风量反演风阻病态改良算法与基于贪心策略的风压传感器优化布置算法，给出了算法步骤及程序框图，能够解决风量反演风阻算法的病态问题，实现布点合理、数量较少的风压传感器的优化布置。将上述算法应用于23条分支、16个节点的通风网络，验证了所有传感器布置方案下的风量反演风阻算法，其中65.2%的方案下的算法都是病态的。笔者应用本文提出的算法给出了角联子网划分方案及9个布点的风压传感器优化布置方案，在此基础上进行了风量反演风阻计算。结果表明:算法在实现了风压传感器优化布置的基础上，准确地反演了风阻，误差控制在10-4以内。算法在现场应用时结合矿井监测监控系统，可以进一步实现实时的通风系统基础参数(风量、风压、风阻)的采集、校验、计算、存储与反馈，是矿井智能通风的基础性研究。

关键词:风量反演风阻;病态算法;角联;贪心策略;传感器优化布置;通风参数智能采集

作者简介:李雨成(1978—)，男，内蒙古赤峰人，教授，博士。E-mail:lyclntu@163.com

Improved algorithm of air quantity calculating resistance based on diagonal subnetwork

(1.College of Safety and Emergency Management Engineering,Taiyuan University of Technology,Taiyuan 030024,China; 2.College of Safety Science and Engineering,Liaoning Technical University,Fuxin 123000,China)

Abstract:The wind resistance of some roadways are not accurate or even impossible to be measured in the process of measuring the resistance of mine ventilation system,which is caused by underground condition.Therefore,the calculation of resistance algorithm based on the airflow of node pressure is studied.It is proposed that using two sets of air volume data and part of node pressure data to inverse the resistance.According to the ill-posed characteristics of the algorithm,it is concluded that the sparsity of resistance inversion equation coefficient matrix and ventilation system topological structure are the two major reasons for the ill-condition of the algorithm.Based on this,the improved algorithm of airflow calculating resistance based on diagonal subnetwork and optimal layout algorithm of wind pressure sensors with algorithm steps and program chart are proposed,which solve the ill-posed problem and achieve the optimal scheme with the least amount and rational layout of wind pressure sensors.The above algorithm is applied to a ventilation network with 23 branches and 16 nodes,the schemes with all combination of wind pressure sensors layout are checked,65.2% of which are ill-posed.Based on the algorithm proposed in this paper,the division scheme of the angular network subnet and the optimal arrangement scheme of the wind pressure sensor of 9 nodes are presented.The results show that the algorithm accurately calculates the wind resistance based on the optimal arrangement of the wind pressure sensor,and the error is controlled within 10-4.Combined with the mine monitoring system,the algorithm can further realize the collection,verification,calculation,storage and feedback of the basic parameters (air volume,air pressure and wind resistance) of the ventilation system in real time,which is the basic research of mine intelligent ventilation.

Key words:air quantity calculating resistance;ill algorithm;diagonal;greedy strategy;sensor optimized layout;intelligent collection of ventilation parameters

李雨成，李俊桥，邓存宝,等.基于角联子网的风量反演风阻病态改良算法[J].煤炭学报,2019,44(4):1147-1153.doi:10.13225/j.cnki.jccs.2018.0741

LI Yucheng，LI Junqiao，DENG Cunbao,et al.Improved algorithm of air quantity calculating resistance based on diagonal subnetwork[J].Journal of China Coal Society,2019,44(4):1147-1153.doi:10.13225/j.cnki.jccs.2018.0741

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

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

1 风量反演风阻算法原理

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

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

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

式中，

和

为节点i和j第1(2)组巷道静压能;Kij定义为“风比平方系数”，无量纲，

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

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

2 算法病态性质分析

2.1 算法病态判别方法

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

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

2.2 算法病态原因分析

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

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

为了不影响井下实际生产，调节构筑物的位置和数量会受到限制，导致获取两组风量的手段有局限性。对于不存在调节构筑物的并联巷道，其风量是大致成相同比例变化的，其风比平方系数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:计算通风网络中所有的角联结构，优先选择分支个数少的角联子网作为待求子网，选择的角联子网D1,D2,…,Dn应尽可能覆盖通风网络除进、回风井外的所有分支(进、回风井由于其拓扑性质，不会包含在角联结构中)，即

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

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

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

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

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

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

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

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

步骤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′集合中删除。

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 整体直接求解

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

表2 节点压能
Table 2 Node pressure

因此理论上需要6个已知压能节点即可满足求解条件。调用Matlab软件中的nchoosek函数计算所有传感器布置组合方案个数为

分别对8 008种传感器布置下的测风求阻算法进行测试，发现仅有2 789组计算结果非病态，即仅有34.8%的计算结果满足精度要求。而且，矿井通风系统是一个动态系统，当系统拓扑结构改变或调节构筑物的方式改变时，满足精度要求的传感器布置能否保证精度仍未可知，这种结果显然是不能接受的。

4.2 划分角联子网求解

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

不考虑节点在通风系统中的实际情况，仅根据节点在角联子网中出现的次数，设定节点权重。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

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

表4 反演风阻结果
Table 4 Calculation of resistanceN·s2/m8

5 结论

(1)系数矩阵的稀疏性及系统拓扑结构是导致风量反演风阻算法病态的主要原因。实例中，整体直接求解风阻的8 008种传感器布置方案中，仅有34.8%的计算结果非病态。而对于划分子网求解，8个子网中仅有2个是病态的。这也说明算法病态原因分析是合理的。

(2)针对基于节点压能的风量反演风阻算法病态特征，提出了基于角联子网的风量反演风阻病态改良算法，并给出了计算实例。与整体直接求解风阻相比，虽然增加了风压传感器布置数目，但却有效解决了病态问题，误差控制在10-4以内，准确地实现了风阻反演。

(3)提出了基于贪心策略的风压传感器优化布置算法，能够制定满足风量反演风阻计算要求的传感器最优布置方案，为风量反演风阻病态改良算法与监测监控系统的结合提供了经济、有效的方案。

[2] 赵丹,刘剑,潘竞涛,等.矿井通风系统风速故障源诊断技术在大明矿的应用研究[J].安全与环境学报,2012,12(3):204-207.

ZHAO Dan,LIU Jian,PAN Jingtao,et al.Application study of air velocity fault source diagnosis technology for ventilation system in Daming Mine[J].Journal of Safety and Environment,2012,12(3):204-207.

[3] 邓立军.矿井通风阻力系数反演研究[D].阜新:辽宁工程技术大学,2014.

DENG Lijun.Study on mine ventilation resistance coefficient inversion[D].Fuxin:Liaoning Technical University,2014.

[4] 司俊鸿.矿井通风系统风流参数动态监测及风量调节优化[D].徐州:中国矿业大学,2012.

SI Junhong.Dynamic monitoring of airflow parameters and air quantity regulation optimization for mine ventilation system[D].Xuzhou:China University of Mining and Technology,2012.

[5] 李雪冰,刘剑,宋莹,等.井巷断面内单点风速与平均风速转换机制[J].安全与环境学报,2018,18(1):123-128.

LI Xuebing,LIU Jian,SONG Ying,et al.On the conversion between the mean airflow velocity and that of the individual point in the underground mine tunnels[J].Journal of Safety and Environment,2018,18(1):123-128.

[6] 王国法,王虹,任怀伟,等.智慧煤矿2025情景目标和发展路径[J].煤炭学报,2018,43(2):295-305.

WANG Guofa,WANG Hong,REN Huaiwei,et al.2025 scenarios and development path of intelligent coal mine[J].Journal of China Coal Society,2018,43(2):295-305.

[8] 刘泽功.利用通风系统调风和阻力测定求算复杂通风网路分支风阻[J].煤矿安全,1991(1):1-7.

[9] 周丽红,吕军,刘小军.测风求阻法的原理及实现[J].西安科技学院学报,2004,24(2):148-150.

ZHOU Lihong,LÜ Jun,LIU Xiaojun.The principle and implementation of calculating resistance through surveying airflow quantity[J].Journal of Xi’an University of Science & Technology,2004,24(2):148-150.

[10] 司俊鸿,陈开岩.基于Tikhonov正则化的矿井通风网络测风求阻法[J].煤炭学报,2012,37(6):994-998.

SI Junhong,CHEN Kaiyan.Measuring airflow & evaluating resistance model of the mine ventilation network based on Tikhonov regularization[J].Journal of China Coal Society,2012,37(6):994-998.

[11] 李雨成,刘天奇,周洋,等.基于风量反演风阻的节点压能解析方法[J].煤炭学报,2015,40(5):1076-1080.

LI Yucheng,LIU Tianqi,ZHOU Yang,et al.Study of node pressure energy analytical method based on inversion from air volume to wind resistance[J].Journal of China Coal Society,2015,40(5):1076-1080.

[13] 吴杰,李明峰,余腾.测量数据处理中病态矩阵和正则化方法[J].大地测量与地球动力学,2010,30(4):102-105.

WU Jie,LI Mingfeng,YU Teng.Ill matrix and regularization method in surveying data processing[J].Journal of Geodesy and Geodynamics,2010,30(4):102-105.

[14] 徐文,陈义,游为.奇异值分解法在病态问题中的应用[J].测绘通报,2016(1):62-63.

XU Wen,CHEN Yi,YOU Wei.Application of SVD method in ill-posed problem[J].Bulletin of Surveying and Mapping,2016(1):62-63.

[15] 刘剑,贾进章,郑丹.基于无向图的角联结构研究[J].煤炭学报,2003,28(6):613-616.

LIU Jian,JIA Jinzhang,ZHENG Dan.Study on diagonal structures in a non-directed graph[J].Journal of China Coal Society,2003,28(6):613-616.

[16] 贾进章.矿井通风系统可靠性、稳定性、安全性理论[M].北京:科学出版社,2016.

[17] 刘剑,贾进章,郑丹.流体网络理论[M].北京:煤炭工业出版社,2002.