CN104268126B - 一种矿井通风系统风道参数的多态自动识别方法 - Google Patents

Abstract

本发明提供一种矿井通风系统风道参数的多态自动识别方法，步骤为：1)按照风阻变化情况，将风道划分为三类；2)Ⅲ类风道构成的单连通区域变为一个虚拟节点，按风流的稳定程度找到最佳生成树，最佳生成树的连支风道布置风量传感器，进而获得所有风道分别在各通风状态下的风量值；结合主要风机工况的在线监测数据，得到各种状态下的降阶网，进而求出各种状态下所有风道风量；3)建立风道风阻识别基本回路矩阵，利用该矩阵的结构特点和Ⅰ、Ⅱ类风道风阻的动态变化特征及最小二乘计算原理，计算Ⅰ类风道固有风阻和Ⅱ类风道动态变化风阻，利用风机工况数据自动完成风机特性曲线拟合。为矿井通风系统仿真模拟和实时调风控风提供准确的基础数据。

Description

一种矿井通风系统风道参数的多态自动识别方法

技术领域

本发明属于矿井通风与安全技术领域，具体涉及一种矿井通风系统风道参数的多态自动识别方法。

背景技术

矿井通风是保障矿井安全的最主要技术手段之一。在矿井生产过程中，必须源源不断地将地面新鲜空气输送到井下各作业地点，以供给人员呼吸，并稀释和排除井下各种有毒、有害的气体和矿尘，从而创造良好的矿内工作环境，保障井下作业人员的身体健康和劳动安全。

矿井通风的主要任务是根据井下各个地点的温度、湿度、有害气体和矿尘浓度实时保证供风质量，满足正常时期和灾变时期各用风地点按时按需供风。

但是，一般大中矿井的通风系统均是由几百条，甚至上千条风道组成的非线性网络流体网络；因此，无论是通风系统模拟分析还是调风控风计算，不仅需要借助计算机软硬件技术，更重要的是，需要获取比较精确的计算基础参数，例如，包括各风道的风阻、自然风压(或火风压或热力风压)和各风机的工作特性曲线；否则，无论是人工计算还是计算机解算，均不可能得到正确的计算结果。

人工阻力测定方法虽然可以获得计算基础数据，然而，随着地面大气压和地温的变化、巷道的变形、掘进面和回采面的推进、通风设施的状态改变和各种车辆和设备的扰动，有许多风道的风阻和自然风压均是随时变动的，因此，由人工阻力测定工作获得的计算基础数据必然滞后矿井的实际情况。此种工作状况，不可能满足按时按需供风，更不可能满足火灾时期的调风和控风需要。

现有的解决方法一般是在风阻变化较大的风道中安装风速传感器，并在解算时将其设为固定风量风道，解算后，根据回路阻力平衡定律和风道阻力定律先后求出其阻力和风阻值。事实上，该种方法存在下列问题：(1)固定风量风道中的实测风量值不一定是由该风道的风阻决定的；(2)只用固定风量风道的实测风量和风机的特性曲线或风机的工况不一定能得到合理的或真实的固定风量风道的风阻值，特别是当固定风量风道的实测风量出现过大或异常时，用该方法解算出来的固定风量风道的风阻值可能为负数，这显然不符合实际。另外，也有人提出了测风求阻，但该方法也存在以下问题：(1)没有利用大部分风道风阻不变的特点(一次测定长久使用)，使得每次计算均需要大量的测风数据；(2)没有利用部分风道风阻的动态变化特征，需要人工调节部风风道风阻，既费时，又影响正常生产；(3)没有形成科学的计算方法，计算过程中很容易出现病态情况，无法得到正确的结果。也就是说，利用上述方法不能有效估计或自动跟踪矿井通风系统各风道风阻的真实值。

发明内容

针对现有技术存在的缺陷，本发明提供一种矿井通风系统风道参数的多态自动识别方法，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种矿井通风系统风道参数的多态自动识别方法，包括以下步骤：

步骤一，对待识别的矿井通风系统的节点由1开始进行顺序编号，设所述矿井通风系统共有J个节点和N条风道，则得到通风系统网络图为G＝[V，E]，其中V为节点集，且|V|＝J；E为风道集，且|E|＝N；

步骤二，按照风阻的变化情况，将所述矿井通风系统的风道划分为三类，分别为：Ⅰ类风道、Ⅱ类风道和Ⅲ类风道；其中，Ⅰ类风道指风阻未知且不随通风状态变化的风道；Ⅱ类风道指风阻未知且随通风状态变化的风道；Ⅲ类风道指在任何通风状态下风阻均已知且不随通风状态变化的风道；

步骤三，选取一种编号方式，对所述通风系统网络图的N条风道由1开始进行顺序编号，由此得到的风道分别记为：L₁、L₂、…、L_N；则可得到如下的M行N列的基本回路矩阵：

其中，M＝N-J+1，代表所述通风系统网络图的基本回路数量；且m1+m2+m3＝M，m1+m2+m3+n1+n2＝N；

并且，第1～m1条风道，即：L₁、L₂、…、L_m1为Ⅰ类风道；

第m1+1～m1+m2条风道，即：L_m1+1、L_m1+2、…、L_m1+m2为Ⅱ类风道；

第m1+m2+1～m1+m2+m3条风道，即：L_m1+m2+1、L_m1+m2+2、…、L_m1+m2+m3风道为Ⅲ类风道；

第m1+m2+m3+1～m1+m2+m3+n1条风道，即：L_m1+m2+m3+1、L_m1+m2+m3++2、…、L_m1+m2+m3+n1为Ⅰ类风道；

第m1+m2+m3+n1+1～N条风道，即：L_{m1+m2+m3+n1+1}、L_{m1+m2+m3+n1+2}、…、L_N为Ⅲ类风道；

步骤四，通过步骤三，可确定m1、m2、m3、n1和n2的值；然后，确定满足下列关系式的通风系统的通风状态数量P；其中，通风状态用G^k表示，且k＝1,2,…,P；

此处，符号[]为取整数符号；

步骤五，将通风系统网络图G的Ⅲ类风道构成的单连通区域变为一个虚拟节点，从而构成一个新网络Gn，在Gn中按风流的稳定程度找到一个最佳生成树，所述最佳生成树的连支风道为较稳定风道，然后，在所述最佳生成树的连支风道中，根据在保证风量平衡的条件下，要求的传感器最少的原则，选取T个连支风道L₁、L₂、…、L_T，使所述最佳生成树的连支风道L₁、L₂、…、L_T为风量传感器的布置风道,即：在所得到的连支风道中分别安装风量传感器；其中，T≤N；L₁、L₂、…、L_T分别为第1条风道、第2条风道…第T条风道；

步骤六，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值：

(1)使用步骤五所安装的风量传感器，分别实际测量在通风状态G¹、G²…G^P下，连支风道L₁、L₂、…、L_T的风道风量值；由此得到L₁、L₂、…、L_T分别在通风状态G¹、G²…G^P下的风道风量值；

(2)以L₁、L₂、…、L_T的风道风量为已知值，根据风量平衡定律，计算通风系统网络图G中的各个悬挂风道分别在通风状态G¹、G²…G^P下的风道风量值；

(3)去除通风系统网络图G中的L₁、L₂、…、L_T风道以及悬挂风道，得到进出风量均已知的降阶网，由风道L₁、、L₂、…、L_T的获得方法可知，所述降阶网中每一个风道均为Ⅲ类风道，风阻为已知值，因此，采用固定风量风道法求出降阶网中各风道的风道风量值；

由此，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值；

步骤七，在通风状态G¹、G²…G^P下，将所述矿井通风系统的状态参数表示为： $T^k=(R^k,Q_f^k,H_f^k,Q^k,H^k,{H_z}^k),$ (k＝1,2,…,P)；

其中：R^k、Q^k、H^k和H_z ^k分别为风道风阻向量、风机风量向量、风机风压向量、风道风量向量、风道风压向量和风网的自然风压向量；

对于N条风道，假设每一条风道的流态指数σ为已知值，且σ≥¹，且各条风道的流态指数可不同；假设每一条风道的自然风压以及风机风压在各种通风状态下的值均为已知值；

则：通过以下步骤求解Ⅰ类风道L₁、L₂、…、L_m1以及L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1的风阻值，其中，由于M＝m1+m2+m3，则L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1即为风道L_M+1、L_M+2、…、L_M+n1；即求解R_j(j＝1,..,m1；M+1,...,M+n1)；

以及，还求解Ⅱ类风道分别在通风状态G¹、G²…G^P下的风阻值，即：求解 $R_j^k(j=m1+1,...,m1+m2;k=\mathrm{1,2},...,P):$

(1)建立下列基本回路定律方程组：

$Q_i^k\left|Q_i^k\right|R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=\mathrm{1,2},...,m1;k=\mathrm{1,2},..,P---\left(\mathrm{EQ}10\right)$

$Q_i^k\left|Q_i^k\right|R_i^k+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=m1+1,...,m1+m2;k=\mathrm{1,2},...,P---\left(\mathrm{EQ}20\right)$

其中：

$b_i^k=-\underset{j=N-n2+1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j+\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}(H_{\mathrm{zj}}^k+H_{\mathrm{fj}}^k),i=\mathrm{1,2},...,M-m3,k=\mathrm{1,2},...,P$

(2)对于每个通风状态G^k，方程EQ10均包含了m1个方程，因此，可将方程EQ10分开并改写成：

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^1{R}_j=c_i^1,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-1)$

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^2{R}_j=c_i^2,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-2)$

……

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^P{R}_j=c_i^P,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-P)$

其中：

$B_{\mathrm{ij}}^1=\frac{B_{\mathrm{ij}}{Q}_j^1{\left|Q_j^1\right|}^{\mathrm{\σ}-1}}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},B_{\mathrm{ij}}^2=\frac{B_{\mathrm{ij}}{Q}_j^2{\left|Q_j^2\right|}^{\mathrm{\σ}-1}}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,B_{\mathrm{ij}}^P=\frac{B_{\mathrm{ij}}{Q}_j^P{\left|Q_j^P\right|}^{\mathrm{\σ}-1}}{Q_P^P{\left|Q_P^P\right|}^{\mathrm{\σ}-1}}$

$c_i^1=\frac{b_i^1}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},c_i^2=\frac{b_i^2}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,c_i^P=\frac{b_i^P}{Q_P^P{\left|Q_P^P\right|}^{\mathrm{\σ}-1}}$

(3)当i＝1时，用方程(EQ10-2)减去方程(EQ10-1)；用方程(EQ10-3)减去方程(EQ10-1)…用方程(EQ10-P)减去方程(EQ10-1)；

当i＝2时，用方程(EQ10-2)减去方程(EQ10-1)；用方程(EQ10-3)减去方程(EQ10-1)…用方程(EQ10-P)减去方程(EQ10-1)；

依此类推，当i＝m1时，用方程(EQ10-2)减去方程(EQ10-1)；用方程(EQ10-3)减去方程(EQ10-1)…用方程(EQ10-P)减去方程(EQ10-1)；

从而消去R₁，R₂，…，R_m1，得到关于R_M+1，R_M+2，…,R_M+n1的m1×(P-1)个方程，即：

$\underset{j=1}{\overset{n1}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{R}_{M+j}=B_i,i=\mathrm{1,2},...,m1\×(P-1)---\left(\mathrm{EQM}10\right)$

其中，系数A_ij和B_i为方程合并得到的系数；

(4)由P的确定条件可知，方程组EQM10的方程个数不小于未知个数，因此，用最小二乘法解方程组EQM10，得到R_M+1，R_M+2，…,R_M+n1；

(5)利用方程EQ10-1，求出R₁，R₂，…，R_m1；

(6)利用方程EQ20，分别求出通风状态G¹、G²…G^P下对应的 $R_j^k(j=m1+1,...,m1+m2;k=\mathrm{1,2},...,P).$

优选的，还包括：

步骤八：利用每台主要风机监测的多状态工况数据(k＝1,2,…,P)，自动拟合得到风量—风压特性曲线：

H_f＝A_f|Q_f-B_f|(Q_f-B_f)+C_f；

其中，A_f、B_f和C_f均为风机特性曲线系数。

本发明的有益效果如下：

本发明提供的矿井通风系统风道参数的多态自动识别方法，不仅可以用较少的通风传感器代替人工阻力测定工作，而且可以实时准确的估计出矿井通风系统中各风道的风阻、自然风压下各风机的工作特性曲线等参数，并且，具有识别过程的计算复杂度低、识别速度快等优点，为矿井通风系统仿真模拟和实时调风控风提供准确的基础数据。

附图说明

图1为本发明提供的一种具体示意的矿井通风系统网络拓扑结构图；

图2为图1的通风状态1下得到的降阶网络拓扑结构图。

具体实施方式

以下结合附图对本发明进行详细说明：

本发明提供一种矿井通风系统风道参数的多态自动识别方法，包括以下步骤：

步骤一，对待识别的矿井通风系统的节点由1开始进行顺序编号，设所述矿井通风系统共有J个节点和N条风道，则得到通风系统网络图为G＝[V，E]，其中V为节点集，且|V|＝J；E为风道集，且|E|＝N；

步骤二，按照风阻的变化情况，将所述矿井通风系统的风道划分为三类，分别为：Ⅰ类风道、Ⅱ类风道和Ⅲ类风道；其中，Ⅰ类风道指风阻未知且不随通风状态变化的风道；Ⅱ类风道指风阻未知且随通风状态变化的风道；Ⅲ类风道指在任何通风状态下风阻均已知且不随通风状态变化的风道；

本步骤中，在对被识别的矿井通风系统的风道进行类别划分时，可通过不同时段和不同通风状态下的多次精细阻力测定和平差计算，按照风阻的变化情况将风道划分为Ⅰ、Ⅱ、Ⅲ类。

步骤三，选取一种编号方式，对所述通风系统网络图的N条风道由1开始进行顺序编号，由此得到的风道分别记为：L₁、L₂、…、L_N；则可得到如下的M行N列的基本回路矩阵：

其中，M＝N-J+1，代表所述通风系统网络图的基本回路数量；且m1+m2+m3＝M，m1+m2+m3+n1+n2＝N；

并且，第1～m1条风道，即：L₁、L₂、…、L_m1为Ⅰ类风道；

第m1+1～m1+m2条风道，即：L_m1+1、L_m1+2、…、L_m1+m2为Ⅱ类风道；

第m1+m2+1～m1+m2+m3条风道，即：L_m1+m2+1、L_m1+m2+2、…、L_m1+m2+m3风道为Ⅲ类风道；

第m1+m2+m3+1～m1+m2+m3+n1条风道，即：L_m1+m2+m3+1、L_m1+m2+m3++2、…、L_m1+m2+m3+n1为Ⅰ类风道；

第m1+m2+m3+n1+1～N条风道，即：L_{m1+m2+m3+n1+1}、L_{m1+m2+m3+n1+2}、…、L_N为Ⅲ类风道；

本步骤中，需要多次尝试不同的风道编号方式，直到使编号后的各风道满足上述的基本回路矩阵。

对于任意一个基本回路矩阵B，其通风状态一定满足下列平衡定律：

风量平衡定律

$Q_j^k=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_i^k,j=\mathrm{1,2},...,N;k=\mathrm{1,2},...,P---\left(1\right)$

基本回路定律：

$\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}(H_j^k-H_{\mathrm{zj}}^k-H_{\mathrm{fj}}^k)=0,i=\mathrm{1,2},...,M;k=\mathrm{1,2},...,P---\left(2\right)$

基本阻力定律：

$H_j^k=R_j^k{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1},j=\mathrm{1,2},...,N;k=\mathrm{1,2},...,P---\left(3\right)$

σ≥1为流态指数，且每一个风道的流态指数可以不一样。

步骤四，通过步骤三，可确定m1、m2、m3、n1和n2的值；然后，确定满足下列关系式的通风系统的通风状态数量P；其中，通风状态用G^k表示，且k＝1,2,…,P；

此处，符号[]为取整数符号；

步骤五，风量传感器最佳布置方案：

将通风系统网络图G的Ⅲ类风道构成的单连通区域变为一个虚拟节点，从而构成一个新网络Gn，在Gn中按风流的稳定程度找到一个最佳生成树，所述最佳生成树的连支风道为较稳定风道，然后，在所述最佳生成树的连支风道中，根据在保证风量平衡的条件下，要求的传感器最少的原则，选取T个连支风道L₁、L₂、…、L_T，使所述最佳生成树的连支风道L₁、L₂、…、L_T为风量传感器的布置风道,即：在所得到的连支风道中分别安装风量传感器；其中，T≤N；L₁、L₂、…、L_T分别为第1条风道、第2条风道…第T条风道；

步骤六，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值：

(1)使用步骤五所安装的风量传感器，分别实际测量在通风状态G¹、G²…G^P下，连支风道L₁、L₂、…、L_T的风道风量值；由此得到L₁、L₂、…、L_T分别在通风状态G¹、G²…G^P下的风道风量值；

(2)以L₁、L₂、…、L_T的风道风量为已知值，根据风量平衡定律，计算通风系统网络图G中的各个悬挂风道分别在通风状态G¹、G²…G^P下的风道风量值；

(3)去除通风系统网络图G中的L₁、L₂、…、L_T风道以及悬挂风道，得到进出风量均已知的降阶网，由风道L_1、、L₂、…、L_T的获得方法可知，所述降阶网中每一个风道均为Ⅲ类风道，风阻为已知值，因此，采用固定风量风道法求出降阶网中各风道的风道风量值；

由此，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值；

步骤七，在通风状态G¹、G²…G^P下，将所述矿井通风系统的状态参数表示为： $T^k=(R^k,Q_f^k,H_f^k,Q^k,H^k{H_z}^k),$ (k＝1,2,…,P)；

其中：R^k、Q^k、H^k和H_z ^k分别为风道风阻向量、风机风量向量、风机风压向量、风道风量向量、风道风压向量和风网的自然风压向量；

对于N条风道，假设每一条风道的流态指数σ为已知值，且σ≥1，且各条风道的流态指数可不同；假设每一条风道的自然风压以及风机风压在各种通风状态下的值均为已知值；

则：通过以下步骤求解Ⅰ类风道L₁、L₂、…、L_m1以及L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1的风阻值，其中，由于M＝m1+m2+m3，则L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1即为风道L_M+1、L_M+2、…、L_M+n1；即求解R_j(j＝1,..,m1；M+1,...,M+n1)；

以及，还求解Ⅱ类风道分别在通风状态G¹、G²…G^P下的风阻值，即：求解 $R_j^k(j=m1+1,...,m1+m2;k=\mathrm{1,2},...,P):$

(1)建立下列基本回路定律方程组：

$Q_i^k\left|Q_i^k\right|R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=\mathrm{1,2},...,m1;k=\mathrm{1,2},..,P---\left(\mathrm{EQ}10\right)$

$Q_i^k\left|Q_i^k\right|R_i^k+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=m1+1,...,m1+m2;k=\mathrm{1,2},...,P---\left(\mathrm{EQ}20\right)$

其中：

$b_i^k=-\underset{j=N-n2+1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j+\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}(H_{\mathrm{zj}}^k+H_{\mathrm{fj}}^k),i=\mathrm{1,2},...,M-m3,k=\mathrm{1,2},...,P$

(2)对于每个通风状态G^k，方程EQ10均包含了m1个方程，因此，可将方程EQ10分开并改写成：

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^1{R}_j=c_i^1,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-1)$

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^2{R}_j=c_i^2,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-2)$

……

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^P{R}_j=c_i^P,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-P)$

其中：

$B_{\mathrm{ij}}^1=\frac{B_{\mathrm{ij}}{Q}_j^1{\left|Q_j^1\right|}^{\mathrm{\σ}-1}}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},B_{\mathrm{ij}}^2=\frac{B_{\mathrm{ij}}{Q}_j^2{\left|Q_j^2\right|}^{\mathrm{\σ}-1}}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,B_{\mathrm{ij}}^P=\frac{B_{\mathrm{ij}}{Q}_j^P{\left|Q_j^P\right|}^{\mathrm{\σ}-1}}{Q_P^P{\left|Q_P^P\right|}^{\mathrm{\σ}-1}}$

$c_i^1=\frac{b_i^1}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},c_i^2=\frac{b_i^2}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,c_i^P=\frac{b_i^P}{Q_P^P{\left|Q_P^P\right|}^{\mathrm{\σ}-1}}$

(3)当i＝1时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

当i＝2时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

依此类推，当i＝m1时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

从而消去R₁，R₂，…，R_m1，得到关于R_M+1，R_M+2，…,R_M+n1的m1×(P-1)个方程，即：

$\underset{j=1}{\overset{n1}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{R}_{M+j}=B_i,i=\mathrm{1,2},...,m1\×(P-1)---\left(\mathrm{EQM}10\right)$

其中，系数A_ij和B_i为方程合并得到的系数；

(4)由P的确定条件可知，方程组EQM10的方程个数不小于未知个数，因此，用最小二乘法解方程组EQM10，得到R_M+1，R_M+2，…,R_M+n1；

(5)利用方程EQ10-1，求出R₁，R₂，…，R_m1；

(6)利用方程EQ20，分别求出通风状态G¹、G²…G^P下对应的 $R_i^k(i=m1+1,...,m1+m2;k=\mathrm{1,2},...,P).$

还包括：

步骤八：利用每台主要风机监测的多状态工况数据(k＝1,2,…,P)，自动拟合得到风量—风压特性曲线：

H_f＝A_f|Q_f-B_f|(Q_f-B_f)+C_f；

其中，A_f、B_f和C_f均为风机特性曲线系数。

拟合得到的风量—风压特性曲线供各种分析、优化使用。

本发明提供的矿井通风系统风道参数的多态自动识别方法，对任意矿井的通风系统，只需要监测风机工况、部分风道的风量，利用本发明的方法，均可以自动识别出矿井通风系统风道参数和风机特性，利用该识别结果不仅可以有效地解决矿井通风系统故障诊断和预警等问题，更重要的是，为实现通风系统自动化奠定了基础。

下面以图1为例，介绍应用本发明提供的矿井通风系统风道参数和风机特性的多态自动识别方法的具体实施例：

步骤一，待识别的矿井通风系统共有7个节点，11条风道，因此，J＝7，N＝11；由此得到通风系统网络图为G＝[V，E]；

步骤二，通过不同时段和通风状态下的多次精细阻力测定和平差计算，按照风阻的变化情况将11条风道分为三类，即：Ⅰ类风道指风阻未知且不随通风状态变化的风道；Ⅱ类风道指风阻未知且随通风状态变化的风道；Ⅲ类风道指在任何通风状态下风阻均已知且不随通风状态变化的风道；

步骤三，通过不断尝试风道编号方式，最终得到图1所示的矿井通风系统网络拓扑结构图，风网拓扑关系和风道分类见表1：

表1风网拓扑关系和风道分类

由图1和表1可知，本例共包括11条风道，风道1为主风机所在风道；风道3和风道4为回采工作面及其顺槽所在风道，风阻随工作面推进而变化，且未知待求；风道5、9、10和11为Ⅲ类风道，风阻固定且已知；风道2、6、7和8为Ⅰ类风道，风阻不随通风状态变化，但未知待求。

按照表1的风道类型，通过调整风道编号或通过列交换，可建立如下满足本发明要求的基本回路矩阵；

其中m1＝2，m2＝2，m3＝1，n1＝3，n2＝3，从而有：

$I_1=\left[\begin{array}{cc}1& \\ & 1\end{array}\right],I_2=\left[\begin{array}{cc}1& \\ & 1\end{array}\right],I_3=\left[1\right],B_{11}=\left[\begin{array}{ccc}0& 1& 1\\ 1& -1& -1\end{array}\right],B_{12}=\left[\begin{array}{ccc}0& 1& 0\\ 1& -1& 0\end{array}\right]$

$B_{21}=\left[\begin{array}{ccc}1& -1& -1\\ 0& 0& -1\end{array}\right],B_{22}=\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 1\end{array}\right],B_{32}=\left[\begin{array}{ccc}1& -1& -1\end{array}\right]$

对于该基本回路矩阵，与图1的对应关系为：基本回路矩阵的行数代表图1中的基本回路个数，并且，每一行中，非0元素所对应的风道构成图1中的一个基本回路，并且，该行非0元素用1和-1表示，1所对应的风道号的方向与该行对应的基本回路方向相同，-1所对应的风道号的方向与该行对应的基本回路方向相反。另外，在该基本回路矩阵中，I₁、I₂和I₃代表单元矩阵，I₁、I₂和I₃的阶数分别为m1、m2和m3；B₁₁的列数为n1，B₁₂的列数为n2。

具体的，对于本例的基本回路矩阵B，共有5行，代表图1中共有5个基本回路；对于每一行，由非0元素所对应的风道号即组成一个基本回路，由此得到的5个基本回路分别为：由风道1-7-8-10构成的第1个基本回路；由风道2-6-7-8-9-10构成的第2个基本回路；由风道3-6-7-8-11构成的第3个基本回路；由风道4-8-11构成的第4个基本回路；由风道5-9-10-11构成的第5个基本回路。另外，基本回路矩阵中，每一行除0外，包括1或-1，1所对应的风道号的方向与该行组成的基本回路的方向相同，-1所对应的风道号的方向与该行组成的基本回路的方向相反。例如，对于第三行所对应的由风道3-6-7-8-11构成的第3个基本回路，则：风道6、3和11的风道方向与第3个基本回路的方向相同，因此，在B中的第三行中，风道号6、3和11对应的元素为1；而风道8和11的风道方向与第3个基本回路的方向相反，因此，在B中的第三行中，风道号8和11对应的元素为-1。

步骤四，通过下式确定通风系统的通风状态数量P：

本例中，由于m1＝2，n1＝3，因此，因此，即，P≥3；本例中，选取P＝3；即利用状态1、状态2和状态3的实测数据进行计算。

步骤五，风量传感器最佳布置方案的确定：

由表1可知，风道5、9、10和11为Ⅲ类风道，因此，风道5、9、10和11构成一个新网络Gn；在Gn中按风流的稳定程度找到一个最佳生成树，该最佳生成树共有5个连支风道，分别为风道1、2、3、4和8；然后，根据在保证风量平衡的条件下，要求的传感器最少的原则，将连支风道8去掉，选取风道1、2、3和4为风量传感器的布置风道；又由于本例中假设风道1为风机风道，由于风道1的风量值可测，因此，不需要额外安装风量传感器，只需要在风道2、3和4中安装风量传感器即可。

步骤六，由于步骤4确定的通风状态数为3，因此，本步骤中，获得所有风道分别在通风状态G¹、G²和G³下的风道风量值，具体确定方法为：

(1)使用步骤五所安装的风量传感器，分别实际测量在通风状态G¹、G²和G³下，连支风道L₁、L₂、L₃、L₄的风道风量值；由此得到风道1、风道2、风道3和风道4分别在通风状态G¹、G²和G³下的风道风量值，结果见表2：

表2状态监测数据表

(2)利用表2的已知数据，结合图1的拓扑结构，根据风量平衡定律，对于每个通风状态，依次计算出风道6、7、8的风量，得到表3所示数据：：

表3由风量平衡定律推算的风道风量

风道号	6	7	8
				状态1的计算风量	14.9672598	14.17108698	4.835471644
状态2的计算风量	14.91823946	14.27796812	4.548753366
				状态3的计算风量	14.7997211	14.45566635	4.200266725

(3)根据表2和表3可知，风道1、2、3、4、6、7、8的风量均已知，在图1中去掉风道1、2、3、4、6、7、8，得到所有进出风量均已知的降阶网，例如状态1对应的降阶网见图2。同样可获得状态2和状态3对应的降阶网。

降阶网中每一个风道均为Ⅲ类风道，风阻为已知值，因此，采用固定风量风道法求解如图2的降阶网，即可获得风道5、9、10、11的风量，结果见表4。

表4求解降阶网获得的风道风量

至此所有风道的风量均已求出，并且满足风量平衡定律。

步骤七，由于假设流态指数σ＝2，而且各条风道的自然风压和风机风压均已知，所以，只剩下求解Ⅰ类风道和Ⅱ类风道的风阻值，即求解R_j(j＝1,2,6,7,8)和步骤为：

(1)建立下列基本回路定律方程组：

$Q_i^k\left|Q_i^k\right|R_i+\underset{j=6}{\overset{8}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k\left|Q_j^k\right|R_j=b_i^k,i=\mathrm{1,2};k=\mathrm{1,2,3}---\left(\mathrm{EQ}11\right)$

$Q_i^k\left|Q_i^k\right|R_i^k+\underset{j=6}{\overset{8}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k\left|Q_j^k\right|R_j=b_i^k,i=\mathrm{3,4};k=\mathrm{1,2,3}---\left(\mathrm{EQ}12\right)$

其中：

$b_i^k=-\underset{j=9}{\overset{11}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k\left|Q_j^k\right|R_j+\underset{j=1}{\overset{11}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}(H_{\mathrm{zj}}^k+H_{\mathrm{fj}}^k),i=\mathrm{1,2,3,4},k=\mathrm{1,2,3}$

对于每个状态k，将已知数据代入(EQ11)，并进行适当整理就得到关于R₁、R₂、R₆、R₇、R₈的方程组如下：

R₁+0.236524705532877 R₇+0.0275389808880287 R₈＝2.75530627321999  (11)

R₂+4.47035018910809 R₆-4.00740536342693 R₇-0.466589143258263 R₈＝1.62868510648224  (12)

R₁+0.239155256261787 R₇+0.0242734716339818 R₈＝2.74003992727891  (21)

R₂+4.6946391478839 R₆-4.30031118840001 R₇-0.4364674365955 R₈＝1.68228259263307  (22)

R₁+0.244154381866162 R₇+0.0206130393120945 R₈＝2.72594350878233  (31)

R₂+4.84816816764002 R₆-4.6253738896757 R₇-0.390502980500613 R₈＝1.68426240106135  (32)

分别由(21)-(11)，(22)-(12)，(31)-(11)，(32)-(12)得到下列最小二乘问题：

0.00263055072891 R₇-0.0032655092540469 R₈＝-0.0152663459410802

0.22428895877581 R₆-0.29290582497308 R₇₊0.030121706662763 R₈＝0.0535974861508299

0.00762967633328501 R₇-0.0069259415759342 R₈＝-0.0293627644376602

0.37781797853193 R₆-0.61796852624877 R₇₊0.07608616275765 R₈＝0.05557729457911

解上述最小二乘问题得到R₆，R₇，R₈；见表5：

表5解最小二乘问题得到风阻值

风道号	6	7	8
				风阻值	1.37300037511283	1.47100029283764	5.86000049466361

将表5的值代入EQ11，得到R₁，R₂，见表6。

表6风阻值

风道号	1	2
			风阻值	2.24599992049141	2.24599992049141

将表5的值代入EQ12，得到风道3和4在各种状态下的风阻值。见表7。

表7风阻值

风道号	3	4
			状态1的风阻值	1.17699977514089	0.980000132709415

状态2的风阻值	0.879999788519993	0.660000108128249
			状态3的风阻值	0.649999809481371	0.34000008297715

由表5、6、7和表1可看出，利用该方法，按照风道分类已求出了所有风道的所有状态的风阻值。

还包括：

步骤八：利用监测的三种状态风机工况点(29.13834684，2658.89511)、(29.19620758，2654.810192)和(29.25538745，2650.539216)，拟合得到风机的风量—风压特性曲线为：

H_f＝A_f|Q_f+B_f|(Q_f+B_f)+C_f＝-13.41665|Q_f-26.536244|(Q_f-26.536244)+2749.73666

其中，A_f、B_f和C_f均为风机特性曲线系数。

由上述曲线形式可知，其优点是单调下降的，不仅在合理工作区与风机实际工况相吻合，而且在任何网络风阻的情况下，均不会出现多解的情况。

综上所述，本发明提供的矿井通风系统风道参数的多态自动识别方法，不仅可以用较少的通风传感器代替人工阻力测定工作，而且可以实时准确的估计出矿井通风系统中各风道的风阻、自然风压下各风机的工作特性曲线等参数，并且，具有识别过程的计算复杂度低、识别速度快等优点，为矿井通风系统仿真模拟和实时调风控风提供准确的基础数据。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (2)

1.一种矿井通风系统风道参数的多态自动识别方法，其特征在于，包括以下步骤：

步骤一，对待识别的矿井通风系统的节点由1开始进行顺序编号，设所述矿井通风系统共有J个节点和N条风道，则得到通风系统网络图为G＝[V，E]，其中V为节点集，且|V|＝J；E为风道集，且|E|＝N；

步骤二，按照风阻的变化情况，将所述矿井通风系统的风道划分为三类，分别为：Ⅰ类风道、Ⅱ类风道和Ⅲ类风道；其中，Ⅰ类风道指风阻未知且不随通风状态变化的风道；Ⅱ类风道指风阻未知且随通风状态变化的风道；Ⅲ类风道指在任何通风状态下风阻均已知且不随通风状态变化的风道；

步骤三，选取一种编号方式，对所述通风系统网络图的N条风道由1开始进行顺序编号，由此得到的风道分别记为：L₁、L₂、…、L_N；则可得到如下的M行N列的基本回路矩阵：

$B={\left[B_{\mathrm{ij}}\right]}_{M\×N}=\left[\begin{array}{ccccc}{I}_1& 0& 0& B_{11}& B_{12}\\ 0& I_2& 0& B_{21}& B_{22}\\ 0& 0& I_3& 0& B_{32}\end{array}\right];$

其中，M＝N-J+1，代表所述通风系统网络图的基本回路数量；且m1+m2+m3＝M，m1+m2+m3+n1+n2＝N；在该基本回路矩阵中，I₁、I₂和I₃代表单元矩阵，I₁、I₂和I₃的阶数分别为m1、m2和m3；B₁₁的列数为n1，B₁₂的列数为n2；

并且，第1～m1条风道，即：L₁、L₂、…、L_m1为Ⅰ类风道；

第m1+1～m1+m2条风道，即：L_m1+1、L_m1+2、…、L_m1+m2为Ⅱ类风道；

第m1+m2+1～m1+m2+m3条风道，即：L_m1+m2+1、L_m1+m2+2、…、L_m1+m2+m3风道为Ⅲ类风道；

第m1+m2+m3+1～m1+m2+m3+n1条风道，即：L_m1+m2+m3+1、L_m1+m2+m3++2、…、L_m1+m2+m3+n1为Ⅰ类风道；

第m1+m2+m3+n1+1～N条风道，即：L_{m1+m2+m3+n1+1}、L_{m1+m2+m3+n1+2}、…、L_N为Ⅲ类风道；

步骤四，通过步骤三，可确定m1、m2、m3、n1和n2的值；然后，确定满足下列关系式的通风系统的通风状态数量P；其中，通风状态用G^k表示，且k＝1,2,…,P；

此处，符号[]为取整数符号；

步骤五，将通风系统网络图G的Ⅲ类风道构成的单连通区域变为一个虚拟节点，从而构成一个新网络Gn，在Gn中按风流的稳定程度找到一个最佳生成树，所述最佳生成树的连支风道为较稳定风道，然后，在所述最佳生成树的连支风道中，根据在保证风量平衡的条件下，要求的传感器最少的原则，选取T个连支风道L₁、L₂、…、L_T，使所述最佳生成树的连支风道L₁、L₂、…、L_T为风量传感器的布置风道,即：在所得到的连支风道中分别安装风量传感器；其中，T≤N；L₁、L₂、…、L_T分别为第1条风道、第2条风道…第T条风道；

步骤六，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值：

(1)使用步骤五所安装的风量传感器，分别实际测量在通风状态G¹、G²…G^P下，连支风道L₁、L₂、…、L_T的风道风量值；由此得到L₁、L₂、…、L_T分别在通风状态G¹、G²…G^P下的风道风量值；

(2)以L₁、L₂、…、L_T的风道风量为已知值，根据风量平衡定律，计算通风系统网络图G中的各个悬挂风道分别在通风状态G¹、G²…G^P下的风道风量值；

(3)去除通风系统网络图G中的L₁、L₂、…、L_T风道以及悬挂风道，得到进出风量均已知的降阶网，由风道L_1、、L₂、…、L_T的获得方法可知，所述降阶网中每一个风道均为Ⅲ类风道，风阻为已知值，因此，采用固定风量风道法求出降阶网中各风道的风道风量值；

由此，获得所有风道分别在通风状态G¹、G²…G^P下的风道风量值；

步骤七，在通风状态G¹、G²…G^P下，将所述矿井通风系统的状态参数表示为： $T^k=(R^k,{Q_f}^k,{H_f}^k,Q^k,H^k,{H_z}^k),$ 其中，k＝1,2,…,P；

其中：和H_z ^k分别为风道风阻向量、风机风量向量、风机风压向量、风道风量向量、风道风压向量和风网的自然风压向量；

对于N条风道，假设每一条风道的流态指数σ为已知值，且σ≥1，且各条风道的流态指数可不同；假设每一条风道的自然风压以及风机风压在各种通风状态下的值均为已知值；

则：通过以下步骤求解Ⅰ类风道L₁、L₂、…、L_m1以及L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1的风阻值，其中，由于M＝m1+m2+m3，则L_m1+m2+m3+1、L_m1+m2+m3+2、…、L_m1+m2+m3+n1即为风道L_M+1、L_M+2、…、L_M+n1；即求解R_j(j＝1,..,m1；M+1,...,M+n1)；

以及，还求解Ⅱ类风道分别在通风状态G¹、G²…G^P下的风阻值，即：求解 $R_j^k(j=m1+1,...,m1+m2;k=\mathrm{1,2},...,P):$

(1)建立下列基本回路定律方程组：

$Q_j^k\left|Q_i^k\right|R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=\mathrm{1,2},...,m1;k=\mathrm{1,2},..,P---\left(\mathrm{EQ}10\right)$

$Q_i^k\left|Q_i^k\right|R_i^k+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j=b_i^k,i=m1+1,...,m1+m2;k=\mathrm{1,2},..,P---\left(\mathrm{EQ}20\right)$

其中：

$b_i^k=-\underset{j=N-n2+1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}{Q}_j^k{\left|Q_j^k\right|}^{\mathrm{\σ}-1}{R}_j+\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}(H_{\mathrm{zj}}^k+H_{\mathrm{fj}}^k),i=\mathrm{1,2},...,M-m3,k=\mathrm{1,2},...,P$

(2)对于每个通风状态G^k，方程EQ10均包含了m1个方程，因此，可将方程EQ10分开并改写成：

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^2{R}_j=c_i^1,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-1)$

$R_i+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^2{R}_j=c_i^2,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-2)$

……

$R_1+\underset{j=M+1}{\overset{M+n1}{\mathrm{\Σ}}}{B}_{\mathrm{ij}}^P{R}_j=c_i^P,i=\mathrm{1,2},...,m1---(\mathrm{EQ}10-P)$

其中：

$B_{\mathrm{ij}}^1=\frac{B_{\mathrm{ij}}{Q}_j^1{\left|Q_j^1\right|}^{\mathrm{\σ}-1}}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},B_{\mathrm{ij}}^2=\frac{B_{\mathrm{ij}}{Q}_j^2{\left|Q_j^2\right|}^{\mathrm{\σ}-1}}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,B_{\mathrm{ij}}^P=\frac{B_{\mathrm{ij}}{Q}_j^P{\left|Q_j^p\right|}^{\mathrm{\σ}-1}}{Q_P^P{\left|Q_p^p\right|}^{\mathrm{\σ}-1}}$

$c_i^1=\frac{b_i^1}{Q_1^1{\left|Q_1^1\right|}^{\mathrm{\σ}-1}},c_i^2=\frac{b_i^2}{Q_2^2{\left|Q_2^2\right|}^{\mathrm{\σ}-1}},...,c_i^P=\frac{b_i^P}{Q_P^P{\left|Q_P^P\right|}^{\mathrm{\σ}-1}}$

(3)当i＝1时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

当i＝2时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

依此类推，当i＝m1时，用方程EQ10-2减去方程EQ10-1；用方程EQ10-3减去方程EQ10-1…用方程EQ10-P减去方程EQ10-1；

从而消去R₁，R₂，…，R_m1，得到关于R_M+1，R_M+2，…，R_M+n1的m1×(P-1)个方程，即：

$\underset{j=1}{\overset{n1}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{R}_{M+j}=B_i,i=\mathrm{1,2},...,m+\×(P-1)---\left(\mathrm{EQM}10\right)$

其中，系数A_ij和B_i为方程合并得到的系数；

(4)由P的确定条件可知，方程组EQM10的方程个数不小于未知个数，因此，用最小二乘法解方程组EQM10，得到R_M+1，R_M+2，…,R_M+n1；

(5)利用方程EQ10-1，求出R₁，R₂，…，R_m1；

(6)利用方程EQ20，分别求出通风状态G¹、G²…G^P下对应的 $R_j^k(j=m1+1,...,m1+m2;k=\mathrm{1,2},...,P).$

2.根据权利要求1所述的矿井通风系统风道参数和风机特性的多态自动识别方法，其特征在于，还包括：

步骤八：利用每台主要风机监测的多状态工况数据其中，k＝1,2,…,P，自动拟合得到风量—风压特性曲线：

H_f＝A_f|Q_f-B_f|(Q_f-B_f)+C_f；

其中，A_f、B_f和C_f均为风机特性曲线系数。