基于超声波全断面测风的矿井风网实时解算方法

宋涛，王建文，吴奉亮，张国群，陈菲，冯雄，李龙清

（1.陕煤集团 神木柠条塔矿业有限公司，陕西 神木 719300；2.西安科技大学 安全科学与工程学院，陕西 西安 710054；3.西安奥文智联信息科技有限公司，陕西 西安 710018）

0 引言

目前我国矿井通风正加速向智能化方向发展[1-2]，矿井通风网络解算软件在此进程中发挥着重要作用。目前具有良好可视化效果与解算性能[3]的通风网络解算软件已广泛用于帮助工程技术人员理解矿井通风系统的复杂网络，如MVSS[4]、VentGIS[5]、VentSim[6]、Ventgraph[7]等。这些软件多是以某一时刻的矿井通风系统状态为假设条件进行计算，主要用于回答“假如在井下采取某种措施后，通风系统将会发生什么变化”等一系列的问题，计算不涉及实时风量，是一种静态解算。但井下风流受风门开闭、行人行车、通风机运行状态的影响，无时无刻不在变化，特别是在巷道冒顶、风门损坏等突发异常情况下，许多巷道的风量、风向都会发生显著变化。这些变化是不能通过以上软件来实时解算的，还需要用风速传感器获取以上变化数据。当前风速传感器对风网的全面感知实时性较差：一方面是由于风速传感器的精度与稳定性易受矿尘的影响，常表现为监测稳定性差；另一方面是因为风速传感器还存在定点监测、安装数量少的问题，不能全面监测风网。随着对矿井智能通风研究的深入，矿井风速监测稳定性差、风网感知不全面的问题已得到部分学者的重视。针对监测稳定性差的问题，刘鹏等[8]改进了压差风速传感器取压装置结构，张巍等[9]研究了风速传感器的数据降噪方法，王恩等[10]实现了测风装置自动在多点移动测取全断面风速，李秉芮等[11]优化了传感器在风网中的布设位置。以上研究虽然在一定程度上提高了风速测量的准确性，但均未解决矿尘对风速传感器精度与稳定性的不利影响。蔡峰等[12]研究了煤尘数量、粒径与超声波波长对超声波声速和衰减系数的影响，但并未探讨研究结果如何应用到超声波测速中，矿尘对风速传感器稳定性影响带来的监测问题仍然存在。针对风网监测不全面的问题，李伟等[13]、谈国文[14]建立了与安全监控系统风速传感器连接的实时解算系统，但未涉及监测值误差及异常风量的处理。鉴此，本文提出了基于超声波全断面测风的矿井风网实时解算方法，并在陕煤集团神木柠条塔矿业有限公司（以下简称柠条塔煤矿）进行了试验，验证了该方法的可靠性。该方法为提升矿井通风智能化水平提供了技术支撑。

1 超声波全断面精准测风

超声波在风流中的速度为其在静止空气中传播速度与空气流速之和。因此，在一定距离内，超声波顺风与逆风传播所需要的时间不同，利用这个时间差可求得风流的速度。超声波测风原理如图1所示。

图1 超声波测风原理Fig.1 Wind measurement principle of ultrasonic

一个风速监测点包括A、B 两个超声波收发点，A、B 横跨巷道断面，两点连线与风流方向夹角为θ。A 点的超声波发送器1、接收器2 向B 点发送超声波信号并接收来自B 点的超声波信号；同理在B 点的超声波发送器2、接收器1 向A 点发送超声波信号并接收来自A 点的超声波信号。

沿A 点到B 点的顺风方向，超声波传播速度为

式中：vc为超声波在静止空气中的传播速度，m/s；vAB为沿AB 方向的风速，m/s。

超声波顺风从A 点传播到B 点所需时间为

式中L为超声波传播距离，m。

同理，沿B 点到A 点的逆风方向，超声波传播速度为

超声波逆风从B 点传播到A 点所需时间为

由式（2）、式（4）可得AB 方向的风速值为

则巷道轴向风速为

基于以上原理的测风装置的实物安装如图2 所示。含有超声波接收器、发送器的2 个装置通过支架固定在图1 所示的巷道两帮A、B 点，此时L、θ是固定值，只要测得tAB、tBA就可求得风速。

图2 超声波全断面测风装置安装实物Fig.2 Installation material object of ultrasonic full-section wind speed measuring device

超声波全断面测风的精度与分辨率主要取决于测风装置CPU 对时间的控制精度，即可识别的最小时长。若取L=6 m，θ=45°，vc=340 m/s 等典型值进行分析，可得到风速v=0.1 m/s 时超声波顺风和逆风的传播时间：

超声波顺风和逆风传播时间差为7.3 μs。常用的C51 单片机CPU 可识别的最小时间粒度可达到1 μs，故超声波时差法测风精度可以达到0.1 m/s，完全能满足矿井最低风速0.15 m/s 的要求；若取v=0.13 m/s，即风速产生0.03 m/s 的变化，此时tAB=17 642.3 μs、tBA=17 651.8 μs，在2 个方向超声波的传播时间改变量均为1.1 μs，都超过了仪器所能捕获的时间1 μs，因此基于本文方法的测风装置的分辨率可以达到0.03 m/s。

采用超声波全断面测风具有不受声波速度、温湿度和气压等参数影响的优点。另外，与传统风速传感器相比，巷道全断面测风的过风口大，不存在粉尘堵塞测风道的问题；风道长度变长，增加了超声波传播时间，降低了捕捉超声波传播时间差的难度。

2 矿井风网实时解算方法

2.1 风流在井巷中流动遵循的物理模型

风流在井巷中流动总是遵循节点流量平衡与回路风压平衡定律。对于含有N条分支、M个节点的矿井通风网络，节点流量平衡与回路风压平衡定律可分别描述如下：

式中：Fq（Q）为节点风量代数和，m3/s，Q为风量真值的列向量，Q=[q1q2…qN]T；bij为分支j（j=1,2,…,N）与节点i（i=1,2,…,M-1）的连接关系，当风流由分支j流入节点i时bij=-1，风流从节点i流出时bij=1，分支j与节点i不相连时bij=0；qj为分支j的准确风量值，m3/s；Fh（R）为回路风压代数和，Pa，R为分支风阻列向量，R=[r1r2…rN]T；ckj为分支j与回路k（k=1,2,…,P，P=N-M+1，为独立回路或余树枝的个数）的关系，分支j在回路k中且与回路同向时ckj=1，与回路反向时ckj=-1，不在回路k中时ckj=0；rj为分支j的风阻，N·s2/m8；Hj为分支j中的通风机风压（若无此项则为0），Pa。

2.2 风网实时解算中的固定风量法

2.2.1 实时解算原理

风网的实时解算是相对于静态解算而言的。静态解算是将通风机性能曲线、分支风阻作为已知数，风量作为未知数，求解式（9）、式（10）的计算过程，解算方法主要有牛顿法、斯拷特-恒斯雷法[15]等，静态解算主要用于矿井通风设计、通风系统预测。与静态解算相比，实时风网解算的已知条件发生了变化，它是在不断采集主要通风机风量、风压实时工况和部分井巷实时风量的条件下解算风网，主要功能是由部分监测数据推演全风网实时状况。牛顿法、斯拷特-恒斯雷法在求解风网时都需要先找到风网的P条余树枝来分别形成一个独立回路，并以余树枝的风量Qy=[qy1qy2…qyP]T为未知数建立一组非线性方程组，Qy中元素qyk的下标增加的前缀y 表示k分支是余树支。根据通风网络理论，风网中任一分支j的风量可通过Qy来计算，即

只要求得余树枝的风量，就可以推算出全风网的风量。以简化的矿井通风网络（图3）为例，分支7，8，4 即是该风网的一组余树枝，只要已知或先求得这3 条分支的风量，就可用式（11）计算其他分支的风量。

图3 简化的矿井通风网络Fig.3 Simplified mine ventilation network

用固定风量法进行实时风网解算就是将通过风速传感器获得实时风量的分支（采用加边法找回路时，将这些分支排在最后）尽可能选作余树枝，这样不仅将监测风量引入式（11）来推演其他分支的风量，也减少了使用牛顿法求Qy时的未知数个数。这种在Qy含有已知风量（固定风量）的情况下进行风网解算的方法，称为固定风量法。理想情况下，当有P个测风装置正好布置在可以同时作为余树枝的分支上时，可以直接使用式（11）求得全风网风量。由于实践中风速传感器多是布置在总回风巷和采掘工作面等重要位置，在数量上也很难覆盖余树的所有分支。所以，在不理想条件下，实时风网解算时也需要先求解Qy，再推演全风网风量。关于求解Qy，特别是保证迭代收敛的方法在许多文献中均有详细介绍[15]，此处不再赘述。下面采用图3 算例进一步说明实时解算原理。

2.2.2 算例分析

图3 算例风网实时解算结果见表1。图3 中各分支的风阻见表1 中的R*，通风机分支8 的风阻为0，设分支1，6，7 有风量监测装置，某时刻监测到的通风机风量、风压分别为87.0 m3/s、2 298 Pa，分支监测风量列向量为表1 中的QM。通风机性能函数J=2 305.8-3.057 218 102q8+0.033 670 849q82-0.000 312 697q83，Q0为用J及R*完成的静态解算结果，相对于实时解算结果而言，Q0可以视作各分支的设计风量。Q1为基于固定风量法的实时解算结果。

表1 图3 算例风网实时解算结果Table 1 Real-time calculation results of the example ventilation network of figure 3

解算Q1时采用加边法找回路，分支按2，3，4，5，1，6，7，8 排序，分支1，6，7，8 后置是为了尽可能将它们选为余树枝。最终分支4，7，8 选为余树枝，形成的3 组独立回路分别为（4，-2，3）、（7，3，-2，-5）和（8，1，2，5，6），回路中负号表示对应分支与余树枝风流方向相反。可以看出：①Q1符合节点流量平衡方程，只要不断地用分支7，8 的监测值进行计算就可以产生动态的Q1。② 监测值还未得到充分利用，尽管分支1，6 有监测风量，但未能被选为余树枝，这些冗余监测分支在QM与Q0中的值还有偏差，这既可能是由测风装置引起的，也可能是由分支风阻的准确度引起的。严格地讲，矿井风网是一个非稳态系统，分支风量无时不在波动，尽管风速传感器采用了全断面精准测风技术，监测风量仍然不会严格满足节点流量平衡定律，因此有必要对实时解算风量进行修正。③R*，Q1未能使以上3 个独立回路达到风压平衡，这是因为井下行人、行车等都会引发风阻波动，因此风阻也需要进行实时修正。

2.3 风量和风阻的实时修正

2.3.1 实时解算结果的最优化修正模型

为实时监测风量，对有冗余监测风量的分支用QM中的元素替换Q1中的元素，替换后的向量记为Q*=[q1*q2*…qN*]T，对于图3 算例需要进行替换的分支是1 和6，替换后的Q*见表2，Q*不符合节点流量平衡定律，需要进行修正。记修正量的列向量为ΔQ，ΔQ=Q-Q*=[Δq1Δq2… ΔqN]T。用fq（ΔQ）表示风量修正量的加权平方和，即

表2 图3 算例风网实时解算修正结果Table 2 Real-time corrected calculation results of the example ventilation network of figure 3

式中wj为第j条分支待修正数据的准确度或可信度，初值越可信，wj值越大。

引入权对角矩阵Wq，则修正量加权平方和的矩阵运算式为fq（ΔQ）=ΔQTWqΔQ。将Q=ΔQ+Q*代入式（9），将Fq（Q）=0 改写为Fq（ΔQ）=0，这样在加权最小二乘法意义下，风量修正的数学模型为

求解式（13）可得到风量的修正量ΔQ。同理可求解到向量R*，R，风阻修正量列向量ΔR=R-R*=[Δr1Δr2… ΔrN]T，风阻修正量的加权平方和fr（ΔR）=ΔRTWrΔR，Wr为风阻修正时的权对角矩阵。以回路风压平衡方程为约束条件，对分支风阻进行修正的模型为

求解式（14）可得到风阻的修正量ΔR。经过风量、风阻修正后的Q，R才是完全符合节点流量平衡与回路风压平衡方程的实时解算结果。

2.3.2 修正模型求解

式（13）和式（14）在数学形式上相同，因此，可将Q，R统一用向量X来表示，则两式统一表示为

将f（ΔX）的极值转变为多元函数φ的极值问题，故需用式（16）对ΔX，Z求一阶偏导数，并令偏导数为零，整理得到

将式（17）代入IΔX+Y=0，求得Z=-（IWX-1IT）-1Y，将Z代回式（17）可得ΔX。由于求风阻的修正量要用到风量，所以先修正风量，再修正风阻，修正时权值可分别取风量、风阻的倒数。

2.3.3 算例与异常值分析

对于图3 算例，在2.2.2 节计算结果的基础之上，经实时修正，得到表2 中的Q和R。经验证，Q，R完全符合式（9）、式（10）的约束，实时解算的风量既与监测值高度吻合，又严格遵循回路风压平衡与节点流量平衡定律。

根据实时解算结果可对2 种故障进行报警。一种是辨识传感器故障或受到瞬时干扰（如车、人经过A、B 测点之间）的无效测值，这时监测风量值之间往往互相矛盾，必然导致计算结果中风量实时修正值过大（如ΔQ在Q中的占比超过30%）。另一种是阻变型故障引发风量异常，此时传感器无故障，但实时解算值Q与设计值Q0之间有较大偏差，这时必然会出现计算结果中风阻值修正量过大的情况，如表2 中解算风量与设计值偏差最大的是分支1，6，最可能的原因就是分支1，6 风阻增加，导致风阻的修正量最大，如果这种情况只是短时发生，致因可能是行人、行车，如果长期存在，则可能是这些分支风阻已发生永久性变化（如巷道变形、断面缩小或有堆积物）；如果实时解算风量与设计风量绝对改变量过大，即出现异常风量（往往是由于巷道冒顶或通风构筑物损毁这种突发的阻变型故障导致），这时必然会对风阻产生一个较大的修正量，因此可从风阻改变量最大的分支中优先定位阻变型故障的位置。

3 柠条塔煤矿精准测风与风网实时解算试验

柠条塔煤矿采用分区对角式通风系统，矿井井田范围大，通风系统复杂，共有4 进2 回6 条井筒，风网共有分支1 319 条，节点945 个。

3.1 精准测风

本次试验在井下27 个测风站安装全断面超声波测风装置，形成智能测风站，对风速进行连续测定。风速测值随时更新，并动态生成1，5，10 min 内的平均值，软件界面如图4 所示，测风装置的分辨率为0.03 m/s。矿井采煤工作面N1204 回风巷的测风装置在2021 年6 月的历史测值统计如图5 所示，在试验点能观察到对0.1 m/s 低风速的响应。5-7 月每旬的人工测风与同时段超声波全断面测风监测值的对比结果如图6 所示。从图5、图6 可看出，在6 月，采煤工作面分风量比较固定，风速在1.4 m/s 上下波动；5-7 月3 个月中，受网络配风不均及回采时采煤工作面通风距离变短、风阻减少、风量增大的影响，采煤工作面配风量呈上升趋势，此时人工测风与超声波测风结果仍是吻合的，证明了超声波全断面测风装置的可靠性。超声波测风装置一经安装到位，理论上不需要调校，除非A、B 基座发生变位。各测点的装置已稳定运行超过6 个月，监测结果准确、可靠。

图4 风量实时测值显示界面Fig.4 Display interface of real-time values of air volume

图5 风速监测值统计Fig.5 Statistical chart of wind speed monitoring values

图6 人工测风与超声波测风对比Fig.6 Comparison between manual wind speed measurement values and ultrasonic wind speed measurement values

3.2 全风网实时解算

27 个测风站主要分布在进回风井、盘区大巷、采煤工作面进回风巷等位置。软件采用云计算架构实现，通过设置Web 服务器、监测主机和计算服务器将软件服务、实时数据采集、风网实时解算集成为一体；前端采用Html 5、WebGL 等技术实现多种信息的三维展示，运行效果如图7 所示。为消除瞬时异常测值对解算结果的影响，测点风量采用1 min内的平均值参与实时解算，解算平均迭代次数约为105，在配置2.4 GHz 主频、10 核CPU 和32 GB 内存的服务器上仅需0.9 s 即可完成1 次解算，解算结果在宏观上可随时间不断更新。尽管本次试验只设置了少量监测点，但与静态网络解算相比，实时解算可形象、准确地展示复杂风网风量变化的实况。随着矿井智能化建设的深入，矿井将建设更多风速精准监测点，使系统的解算能力进一步提升。

图7 基于WebGL 的风网实时解算前端显示界面Fig.7 The front-end display interface of real-time calculation results of mine ventilation network based on WebGl technology

4 结论

（1）利用超声波在两点间顺风、逆风传播的时间差实现巷道全断面测风，风速测定结果与声速无关，具有不受声波速度、温湿度和气压等参数影响的优点，精度高，稳定性好。全断面测风避免了传统风速传感器的风道易受矿尘堵塞的难题，测风装置的分辨率可达0.03 m/s，试验点可观察到对0.1 m/s 风速的响应和超过6 个月不用调较的稳定运行期。

（2）全风网实时解算方法通过不断采集主要通风机风量、风压实时工况和部分井巷实时风量解算风网，具有风量、风阻双实时解算能力。使用固定风量法将监测风量融入通风网络中，解算得到全风网实时风量，采用最优化方法修正风量与风阻解算结果，解决了监测冗余分支引起的节点风量不平衡与分支风阻波动引起的回路风压不平衡问题。通过算例验证了该方法实时解算结果与监测值高度吻合，同时又严格遵循回路风压平衡与节点流量平衡的约束。

（3）在2.4 GHz 主频CPU 的服务器上对柠条塔煤矿含1 319 条分支、945 个节点的风网进行实时解算，1 次解算仅用时0.9 s，解算迭代收敛次数约为105，解算用时少，速度快，结果稳定可靠，可全面、形象地展示复杂矿井通风系统的真实情况。