关于煤矿瓦斯和煤尘的监测与控制的数学模型

1、.高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛的题目是： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）：

2、参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：;1高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：关于煤矿瓦斯和煤尘的监测与控制的数学模型摘要：本文是对煤矿开采的实际情况进行研究。通过Matlab对所给数据进行处理，得出所需要数据，并导入Excel中。利用此数据分别求所要求的结果，最后根据实际情况约束要求得出

3、优化模型，求出最佳通风量，从而使实际问题得到解决。对于问题（1）根据所给数据对总回风巷一天中早班、中班、晚班取其平均根据公式算出相对瓦斯涌出量及绝对瓦斯涌出量，与“高瓦斯矿井”及“低瓦斯矿井”的判别条件进行比较，判断出该矿井类型为高瓦斯矿井。 对于问题（2），我们首先对瓦斯浓度进行相应计算，判断瓦斯爆炸情况，借助0-1变量进行记录，爆炸记为1，不爆炸记为0。然后借助瓦斯浓度与煤尘爆炸下限的关系，利用数据拟合的方法，用多项式函数逼近，得到拟合曲线，从而得到瓦斯浓度与煤尘爆炸下限的关系，然后借助0-1变量对煤尘爆炸情况进行统计，根据瓦斯和煤尘的爆炸情况判断煤矿发生事故的可能性大小约为22.22%。

4、对于问题（3）根据所给数据找出瓦斯、煤尘与风速的关系，通过瓦斯、煤尘的约数条件转化到风速的约束要求，进而得出各处风量的限定条件，据此建立最佳通风量的优化模型，然后利用Lingo求解，得出最佳通风量为1576.621 ，工作面的风量为483，工作面的风量为916，通风机的额定风量为150。本文很好的解决了矿井类型鉴别、发生事故可能性计算及最佳总通风量计算的问题，且本文所述模型适用性较强，稍加修改即可用于其它矿井，大大推广了本模型的适用范围。关键词：瓦斯含量 煤尘浓度 风量 0-1变量 Lingo Matlab一、 问题的重述煤矿开采涉及多方面问题，以下几个问题是煤矿开采中经常遇到的几个问题，试结

5、合实际建立模型说明以下几个问题：（1）当矿井相对瓦斯涌出量小于或等于10 且矿井绝对瓦斯涌出量小于或等于 时为低瓦斯矿井，当矿井相对瓦斯涌出量大于10或矿井绝对瓦斯涌出量大于 时为高瓦斯矿井。结合所给数据判断该煤矿是“低瓦斯矿井”还是“高瓦斯矿井”。（2）瓦斯浓度与煤尘爆炸下限存在一定关系，找出其中的关系，并据此判断该煤矿不安全的程度（即发生爆炸事故的可能性）。 （3）为了保障安全生产，利用两个可控风门调节各采煤工作面的风量，通过一个局部通风机和风筒实现掘进巷的通风。根据各井巷风量的分流情况、对各井巷中风速的要求（见下表），以及瓦斯和煤尘等因素的影响，确定该煤矿所需要的最佳（总）通风量，以及两

6、个采煤工作面所需要的风量和局部通风机的额定风量（在实际中，井巷可能会出现漏风现象）。井巷名称允许风速最低最高主要进、回风巷运输巷，采区进、回风巷采煤工作面、掘进中的煤巷和半煤岩巷0.250.25864二、 问题的分析问题（1）要求确定该煤矿是高瓦斯矿井还是低瓦斯矿井，而经分析知用总回风巷在早中晚班的平均值可近似用来表示这一天的瓦斯涌出量，对此我们利用excel将附表2数据进行计算，得到该煤矿在30天内矿井相对瓦斯涌出量均大于10 ，故可知该煤矿为高瓦斯矿井（具体计算数据见附录1）。问题（2）要求该煤矿发生事故的可能性有多大，瓦斯爆炸情况可根据附表2数据得到。对于煤尘我们借助瓦斯浓度与煤尘爆炸下

7、限的关系，故利用数据拟合的方法，用多项式函数逼近，得到拟合曲线，从而得到瓦斯浓度与煤尘爆炸下限的关系，然后借助0-1变量对瓦斯爆炸情况进行统计，据此判断煤矿发生事故的可能性大小为22.22%。问题（3）要求该煤矿的最佳总通风量，需要求解瓦斯浓度及煤尘浓度与风速的关系，通过建立优化模型利用Lingo即可求解。三、 模型假设与符号说明1 模型假设（1）实际中瓦斯爆炸界限值并不是固定不变的，它还受温度、压力以及煤尘和其它可燃性气体、惰性气体的混入等因素的影响。此模型假设此环境中不受这些因素影响。（2）假设一天内井内环境封闭，不受外界因素的影响，用早中晚三次测量结果的平均值来估计这一天内各个值的分布情

8、况。（3）假设当瓦斯超过警报浓度或当煤尘达到爆炸下限浓度时爆炸事故发生。（4）井巷不可能会出现漏风现象，且因为瓦斯含量很小，所以假设进风口和出风口的风量相同。（5）假设总进风巷和各进风巷中无瓦斯和煤尘掺入，即全部为新鲜风流。（6）假设在各工作面进风巷中的风量与工作面风量相同，即每次进入一个工作面的循环风量达到平衡，不会出现因风量不等产生的局部气压问题。（7）假设从风筒中流出的风在到达掘进面前经过足够的缓冲，使经过掘进面的风匀速。即不存在局部高速风的问题。2符号说明符号说 明发生爆炸事故的有利样本点数样本点总数瓦斯浓度相对涌出量绝对涌出量风筒横截面积风速煤单位时间内的产量工作面I的风量工作面II

9、的风量回风巷I的风量工作面II的风量通风机的风量工作面I的风速工作面II的风速掘进巷的风速掘进巷的风量四、 模型建立与求解4.1问题14.1.1模型建立（1）绝对瓦斯涌出量对总回风巷一天在早中晚班的瓦斯浓度取平均值作为这一天的瓦斯浓度，根据公式 -（1）可求得绝对瓦斯涌出量。（2）相对瓦斯涌出量由（1）可得绝对瓦斯涌出量，然后借助公式 -（2）即可得相对瓦斯涌出量。（3）矿井类型判断 当矿井相对瓦斯涌出量小于或等于10 且矿井绝对瓦斯涌出量小于或等于 时为低瓦斯矿井，当矿井相对瓦斯涌出量大于10或矿井绝对瓦斯涌出量大于 时为高瓦斯矿井。据此建立判定“高瓦斯井”还是“低瓦斯井”的依据如下 为低瓦

10、斯矿井。 -（3）或 为高瓦斯矿井。 -（4）（4）概率模型建立： 4.1.2模型求解因无法得知具体每处工作面的产煤量，所以对总回风巷的数据进行处理，根据所建立判定“高瓦斯井”还是“低瓦斯井”的判据判断次井的类型。对附表2数据进行计算，步骤如下（计算程序见附录）：根据早班、中班、晚班时测得的瓦斯浓度数据算其平均值，算出总回风巷瓦斯平均浓度。对早班、中班、晚班的风速取平均值，得平均风速。通过上面两步得到的结果、，根据公式（1）算出矿井绝对瓦斯涌出量，根据公式（2）算出矿井相对瓦斯涌出量。根据和的数值得该矿井为高瓦斯矿井。对表中所给的30天的数据都进行Step1到Step3的处理，得出30天中每天

11、数据显示的瓦斯井的类型，并计数。根据所建立的概率模型判断次井的类型。通过上述的计算得出，此矿井为“高瓦斯井”的概率为1。即可判断此矿井为“高瓦斯井”。4.2问题24.2.1模型建立（1）瓦斯爆炸判断依次判断工作面、工作面、掘进工作面、回风巷、回风巷、总回风巷是否符合瓦斯爆炸要求并用0-1变量进行记录，如发生爆炸则记为1，否则记为0，当瓦斯浓度大于1%时假定为瓦斯爆炸，否则认为不爆炸，即： （2）瓦斯浓度与煤尘爆炸下限关系确定根据附录2表1，将瓦斯浓度与煤尘爆炸下限关系用图表表示，然后用多项式函数逼近得其逼近曲线为三次多项式函数时拟合度较好，符合计算要求，如图：（3）煤尘爆炸判断在判断瓦斯不爆炸

12、的情况下，根据瓦斯浓度与煤尘爆炸的关系，对应相应瓦斯浓度确定煤尘是否满足爆炸范围，并用0-1变量进行记录，对这组数据进行记录，若满足爆炸范围则记为1，否则记为0。即对记录的数据进行求和，和数即为发生爆炸的样本点数（见附录3）。（4）概率模型建立 -（5）其中各变量满足的条件为：4.2.2模型求解由上述建立的模型，可求得瓦斯爆炸情况及煤尘爆炸情况，求解算法如下（具体代码见附录）：根据已知关系，绘出瓦斯浓度与没尘爆炸下限的散点图，利用数据拟合的方法，得到瓦斯浓度与煤尘爆炸下限的关系。对工作面、工作面、掘进工作面、回风巷、回风巷、总回风巷在一天中三此次量的瓦斯浓度和煤尘浓度算其平均值，针对瓦斯浓度，

13、分别判断30天否达到爆炸范围,并将结果记录下来。根据瓦斯浓度与煤尘爆炸下限的关系，根据矿井中的当日的瓦斯浓度算出此瓦斯浓度下煤尘浓度的爆炸下限。再根据此时煤矿中实际的煤尘含量判断煤尘是否达到爆炸范围，判断30天各自的情况，并记录下计算结果。根据概率模型（5）即得矿井发生事故的可能性大小。通过计算得到结果如下：煤矿发生安全事故的可能性大小为： ，所以此煤矿不安全的程度约为22.22% 。4.3问题34.3.1模型建立该问题是一个最优化问题，因此建立最优化模型，求得最佳总风量。目标函数：1 风量影响的约束条件：2. 风速影响的约束条件：3为确定瓦斯浓度和煤尘浓度所确定的约束条件，首先运用Matla

14、b，对90个样本中6个井巷中的每组数据按风速进行排序，排序的目的是使后面拟合出来的函数更符合客观实际。把排序后风速相同的样本对应的瓦斯浓度和煤尘浓度求平均值，把多个样本点划归为一个样本点，用1次多项式拟合出风量与瓦斯浓度的函数关系F1, 从而确定在风量一定的情况下瓦斯浓度值，再根据煤矿安全规程第一百六十八条关于瓦斯浓度的规定确定瓦斯浓度影响的约束条件：瓦斯浓度影响的约束条件：4用排序拟合后的数据，通过Matlab求出风量与煤尘浓度的函数关系F2，从而确定在风量一定的情况下煤尘浓度值。在根据瓦斯浓度与煤尘爆炸下限的关系f确定煤尘的爆炸下限。通过比较煤尘浓度值与爆炸下限的关系，确定煤尘影响的约束条

15、件。煤尘浓度确定的约束条件：4.3.2模型求解要求最佳总通风量，则需工作面、工作面、掘进工作面、回风巷、回风巷的通风量满足一定要求，具体求解步骤如下：求出瓦斯浓度及煤尘浓度与风速的关系。求出瓦斯需满足的约束条件。求出煤尘需满足的约束条件。根据目标函数由Lingo求解（Lingo程序见附录6），得结果为1576.621，工作面的风量为483，工作面的风量为916，通风机的额定风量为150五、 模型评价与推广5.1优点本文优点在于在问题(1)中利用excel进行计算显然使运算相对简单，免去了编程序的复杂，在问题（2）中用0-1变量进行记录，求和数时明显简化计算。并在利用Matlab操作时进行文件操

16、作，减少了在数据输入时的麻烦，加快了编程效率。5.2不足之处由于数据的有限在问题（1）中我们用总回风巷的数值近似计算矿井绝对瓦斯涌出量及相对瓦斯涌出量，与实际存在一定的误差，而且用均值并不能完全充分反映每一时刻各个数据的实际值。问题（3）中在计算总通风量时忽略了许多外部因素的影响，这可能使运算结果与实际值存在误差，且因为所掌握数据的有限性，不能对误差进行掌握和控制，使得模型结果与实际不符合。5.3推广此模型稍加改进即可应用于其它矿井用于评估其为“高瓦斯矿井”或“低瓦斯矿井”，并能计算其发生事故可能性，借此采取相应的措施，预防事故的发生，保障人员的安全。附录：1总回风巷一天平均风速总回风巷一天平

17、均瓦斯浓度日产量矿井绝对瓦斯涌出量矿井相对瓦斯涌出量5.150.6759710.0193324.167245.360.65975.310.635975.140.66029.82723.506455.390.646025.320.626025.190.6863910.1266722.820665.260.656395.050.636395.170.686169.61093322.467125.140.586165.360.586165.210.616109.57413322.601235.370.586105.030.656105.080.6658810.2077324.998535.010.6

18、85885.070.685885.190.65829.709224.022765.180.625825.290.645825.340.596059.494422.598245.010.626055.130.636055.170.615859.48826723.355735.280.635855.190.585855.180.656209.73646722.613735.230.626205.210.66205.10.626169.56816722.367145.170.616165.080.646165.060.66089.02416721.373035.130.586085.280.5760

19、85.160.6861210.2050724.011925.360.626125.10.666125.130.576069.82323.341785.220.696065.160.646065.320.555929.361722.77175.110.625925.260.625925.20.646099.70922.957245.060.656095.070.616095.140.6860510.4274324.8190150.76055.270.656055.380.65819.79883324.286275.20.635815.310.625815.20.686169.51622.2451

20、95.20.576165.20.586165.280.666129.72422.885.240.596125.080.626125.160.626019.78523.444935.080.666015.210.626015.230.646169.740222.76935.280.616165.20.616165.350.5958810.37425.405715.230.685885.380.685885.120.676259.60556722.131235.050.626255.240.586255.240.575989.6623.261545.160.625985.350.655985.23

21、0.66189.59322.352625.40.616185.270.66185.120.616179.83013322.942295.260.646175.060.666175.110.636059.45146722.496055.190.66055.110.616055.10.715869.77624.022945.360.545865.140.635865.230.656209.912923.023515.230.66205.110.666202表1：瓦斯浓度对煤尘爆炸下限的影响系数 空气中的瓦斯深度%00.500.751.01.502.03.04.0k10.750.600.500.35

22、0.250.10.053Data=xlsread('1.xlsx');Data_1=xlsread('1_1.xlsx');Sz=size(Data)Pos=1:Sz;P=0;x=Data_1(1,:);y=Data_1(2,:);f=polyfit(x,y,3);S=5; for j=1:6for i=1:Sz(1) F(Data(i,3*j),f) if Data(i,1+3*j)<F(Data(i,3*j),f) Data(i,1+3*j)=0; else Data(i,1+3*j)=1; end if Data(i,3*j)<1 Data(i

23、,3*j)=0; else Data(i,3*j)=1; endendendfor i=1:Sz(1) Temp=0; for j=1:6 Temp=Temp+Data(i,1+3*j)+Data(i,3*j); end if Temp=0 Pos(i)=0; Data(i,1)=0; else Pos(i)=1; Data(i,1)=1; end P=P+Pos(i);end Pxlswrite('2.xlsx',Data);fxlswrite('3_3.xlsx',f)4井的爆炸情况工作面1瓦斯工作面1煤尘工作面瓦斯工作面2煤尘挖掘面瓦斯挖掘面煤尘回风巷1瓦

24、斯回风巷1煤尘回风巷2瓦斯回风巷2煤尘总回风巷瓦斯总回风巷煤尘00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000100010000000010000000000000000000000000000000000000000000000000000000100000000100010000000010001000000001000100100000100000000000000000000000000000000000000

25、000000000000000000000000000000100000000100010010000010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000000010000000000000000000000000000000000000000001001000001000000000000000000000000000000000000000

26、000000000000000000000000000000000000000000000000000000010000000010001001000001000100000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000001000000000000000000000000000000000000000000000000000000000000000000

27、001000000001000100000000100000000000000000000000000000000000000000000000000000000000000000000100000000100000000000000000000000000000000000000000000000000000001000000001000000000000000000000000000000000000000000000000000000010010000010000000000000000000000000000010000000010000000000000000000000000000

28、05Data=xlsread('3.xlsx');Fit1=ones(6,2);Fit2=Fit1;Sz=size(Data);for k=1:6 for i=1:Sz(1) C=0; for j=i:Sz(1) if Data(i,k*3-1)=Data(j,k*3-1) C=C+1; else break; end end if C=0 Temp1=0; Temp2=0; for j=i:i+C-1 Temp1=Temp1+Data(j,k*3); Temp2=Temp2+Data(j,k*3+1); end for j=i:i+C-1 Data(j,k*3)=Temp1/

29、C; Data(j,k*3+1)=Temp2/C; end end i=i+C-1; endendfor i=1:6 if i=6 Fit1(i,:)=polyfit(Data(:,i*3-1)*5*60,Data(:,i*3),1); figure; plot(Data(:,i*3-1)*5*60,Data(:,i*3); hold on; plot(Data(:,i*3-1)*5*60,polyval(Fit1(i,:),Data(:,i*3-1)*5*60),'r'); hold off; Fit2(i,:)=polyfit(Data(:,i*3-1)*5*60,Data

30、(:,i*3+1),1); figure; plot(Data(:,i*3-1)*5*60,Data(:,i*3+1); hold on; plot(Data(:,i*3-1),polyval(Fit2(i,:),Data(:,i*3-1)*5*60),'r'); hold off; else Fit1(i,:)=polyfit(Data(:,i*3-1)*4*60,Data(:,i*3),1); figure; plot(Data(:,i*3-1)*4*60,Data(:,i*3); hold on; plot(Data(:,i*3-1)*4*60,polyval(Fit1(

31、i,:),Data(:,i*3-1)*4*60),'r'); hold off; Fit2(i,:)=polyfit(Data(:,i*3-1)*4*60,Data(:,i*3+1),1); figure; plot(Data(:,i*3-1)*4*60,Data(:,i*3+1); hold on; plot(Data(:,i*3-1)*4*60,polyval(Fit2(i,:),Data(:,i*3-1)*4*60),'r'); hold off; endendxlswrite('3_1.xlsx',Fit1);xlswrite('

32、3_2.xlsx',Fit2);xlswrite('3_4.xlsx',Data);6sets:aaa/1.2/:Fit1_1,Fit1_2,Fit1_3,Fit1_4,Fit1_5,Fit1_6,Fit2_1,Fit2_2,Fit2_3,Fit2_4,Fit2_5,Fit2_6;bbb/1.4/:f;endsetsdata:Fit1_1=-0.001088759,1.287519738;Fit1_2=-0.001903513,1.832843745 ;Fit1_3=-0.000335049,0.409851395 ;Fit1_4=-0.00096417,1.211194155 ;Fit1_5=-0.002629594,2.260884817;Fit1_6=-0.000397958,1.24627818 ;Fit2_1=0.008252953,3.216582546; Fit2_2=0.005201876,5.012531394;Fit2_3=0.008