0.   引言

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

1.   超声波全断面精准测风

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

图  1  超声波测风原理

Figure  1.  Wind measurement principle of ultrasonic

下载: 全尺寸图片 幻灯片

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

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

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

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

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

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

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

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

则巷道轴向风速为

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

图  2  超声波全断面测风装置安装实物

Figure  2.  Installation material object of ultrasonic full-section wind speed measuring device

下载: 全尺寸图片 幻灯片

超声波全断面测风的精度与分辨率主要取决于测风装置CPU对时间的控制精度，即可识别的最小时长。若取L=6 m，θ=45º，v_c=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的变化，此时t_AB=17 642.3 μs、t_BA=17 651.8 μs，在2个方向超声波的传播时间改变量均为1.1 μs，都超过了仪器所能捕获的时间1 μs，因此基于本文方法的测风装置的分辨率可以达到0.03 m/s。

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

2.   矿井风网实时解算方法

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

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

式中：F_q（Q）为节点风量代数和，m³/s，Q为风量真值的列向量，Q=[q₁ q₂ … q_N]^T；b_ij为分支j（j=1, 2, …, N）与节点i(i=1, 2, …, M−1)的连接关系，当风流由分支j流入节点i时b_ij=−1，风流从节点i流出时b_ij=1，分支j与节点i不相连时b_ij=0；q_j为分支j的准确风量值，m³/s；F_h（R）为回路风压代数和，Pa，R为分支风阻列向量，R=[r₁ r₂ … r_N]^T；c_kj为分支j与回路k(k=1, 2, …, P，P=N−M+1，为独立回路或余树枝的个数)的关系，分支j在回路k中且与回路同向时c_kj=1，与回路反向时c_kj=−1，不在回路k中时c_kj=0；r_j为分支j的风阻，N·s²/m⁸；H_j为分支j中的通风机风压（若无此项则为0），Pa。

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

2.2.1   实时解算原理

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

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

图  3  简化的矿井通风网络

Figure  3.  Simplified mine ventilation network

下载: 全尺寸图片 幻灯片

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

2.2.2   算例分析

图3算例风网实时解算结果见表1。图3中各分支的风阻见表1中的R^*，通风机分支8的风阻为0，设分支1，6，7有风量监测装置，某时刻监测到的通风机风量、风压分别为87.0 m³/s、2 298 Pa，分支监测风量列向量为表1中的Q_M。通风机性能函数J=2305.8−3.057218102q₈+0.033670849 q₈²−0.000312697 q₈³，Q₀为用J及R^* 完成的静态解算结果，相对于实时解算结果而言，Q₀可以视作各分支的设计风量。Q₁为基于固定风量法的实时解算结果。

表  1  图3算例风网实时解算结果

Table  1.  Real-time calculation results of the example ventilation network of figure 3

分支号	R* 元素/ （${\rm{N}} \cdot {{\rm{s}}^{2}} \cdot {{\rm{M}}^{-8}}$）	QM元素/ （${\rm{m} }^{3}\cdot {\rm{s} }^{-1}$）	Q0元素/ (${\rm{m}}^{3}\cdot {\rm{s}}^{-1}$)	Q1元素/ （${\rm{m}}^{3}\cdot {\rm{s}}^{-1}$)
1	0.0075	83.0	88.0	87.0
2	0.3800	−	47.5	47.0
3	0.5000	−	40.5	40.0
4	0.0500	−	27.8	27.3
5	0.2000	−	75.3	74.3
6	0.0057	85.0	88.0	87.0
7	7.0769	12.7	12.7	12.7
8	0	87.0	88.0	87.0

下载: 导出CSV 

| 显示表格

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

2.3   风量和风阻的实时修正

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

为实时监测风量，对有冗余监测风量的分支用Q_M中的元素替换Q₁中的元素，替换后的向量记为Q^*=[q₁^* q₂^* … q_N^*]^T，对于图3算例需要进行替换的分支是1和6，替换后的Q^*见表2，Q^*不符合节点流量平衡定律，需要进行修正。记修正量的列向量为ΔQ，ΔQ=Q−Q^*=[Δq₁ Δq₂ … Δq_N]^T。用f_q(ΔQ)表示风量修正量的加权平方和，即

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

表  2  图3算例风网实时解算修正结果

Table  2.  Real-time corrected calculation results of the example ventilation network of figure 3

分支号	Q*元素/ （m3·s−1）	Q元素/ （m3·s−1）	R元素/ （N·s2·m−8）	|Q−Q0|元素/ （m3·s−1）	|R−R*|/R 元素/%
1	83.0	85.80355	0.018582	2.20	59.6
2	47.0	46.1907	0.449975	1.31	15.6
3	40.0	39.61286	0.589756	0.89	15.2
4	27.3	27.15569	0.046953	0.64	6.5
5	74.3	73.34639	0.203945	1.95	1.9
6	85.0	85.80355	0.014123	2.20	59.6
7	12.7	12.45717	7.293318	0.24	3.0
8	87.0	85.80355	0	2.20	−

下载: 导出CSV 

| 显示表格

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

求解式（13）可得到风量的修正量ΔQ。同理可求解到向量R^*，R，风阻修正量列向量ΔR=R−R^*=[Δr₁ Δr₂ … Δr_N]^T，风阻修正量的加权平方和f_r(ΔR)=ΔR^TW_rΔR，W_r为风阻修正时的权对角矩阵。以回路风压平衡方程为约束条件，对分支风阻进行修正的模型为

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

2.3.2   修正模型求解

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

式中W_X为权对角矩阵。

设矩阵I的行数为m，当X表示Q时，矩阵I的元素a_ij=b_ij，m=M−1,向量D=0；当X表示R时，a_kj=c_kj(q_j)²，m=P，向量D 的元素为${d}_{k}={\displaystyle \sum _{j=1}^{N}{c}_{kj}（{r}_{j}{}^{*}（{q}_{j}{）}^{2}-{H}_{j})} $。

为求式（15），由条件极值原理引入乘子向量Z，建立拉格朗日函数：

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

将式（17）代入IΔX+Y=0，求得Z=−(IW_X⁻¹I^T)⁻¹Y，将Z代回式（17）可得ΔX。由于求风阻的修正量要用到风量，所以先修正风量，再修正风阻，修正时权值可分别取风量、风阻的倒数。

2.3.3   算例与异常值分析

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

根据实时解算结果可对2种故障进行报警。一种是辨识传感器故障或受到瞬时干扰（如车、人经过A、B测点之间）的无效测值，这时监测风量值之间往往互相矛盾，必然导致计算结果中风量实时修正值过大（如ΔQ在Q中的占比超过30%）。另一种是阻变型故障引发风量异常，此时传感器无故障，但实时解算值Q与设计值Q₀之间有较大偏差，这时必然会出现计算结果中风阻值修正量过大的情况，如表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  风量实时测值显示界面

Figure  4.  Display interface of real-time values of air volume

下载: 全尺寸图片 幻灯片

图  5  风速监测值统计

Figure  5.  Statistical chart of wind speed monitoring values

下载: 全尺寸图片 幻灯片

图  6  人工测风与超声波测风对比

Figure  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的风网实时解算前端显示界面

Figure  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，解算用时少，速度快，结果稳定可靠，可全面、形象地展示复杂矿井通风系统的真实情况。