全国可采侏罗系煤炭资源占我国煤炭总储量的67%以上，主要分布在晋陕蒙、新甘宁6个省份，煤层开采普遍面临顶板水威胁[1]。顶板水主要是随工作面开采形成的导水裂缝带进入矿井，因此，工作面涌水随采动变化规律的研究是矿井水害防治最基础的工作，研究结果能够为涌水量预测、矿井防排水系统设计以及水害防治技术方案的选择提供参考[2-3]。

目前，在煤矿顶板水害防治领域，常用以计算工作面涌水过程的方法主要包括解析法和数值法，它们是以地下水渗流运动的基本微分方程为理论基础的[4-9]。受勘探水文地质参数与实际值误差较大原因影响，解析法计算工作面采动涌水量得到的结果偏差往往偏大;边界条件选取的偏差，造成数值法预测结果与实际值也存在较大差距[2,10-12]。并且，两种方法计算的涌水量均是静态量，不能体现涌水过程随采掘工程的变化特征，这与实际涌水规律不符[2]。

涌水量随采掘工程动态变化的研究大多集中在矿井尺度上，主要是结合矿井工作面接续进行的叠加计算[3,13]，真正研究某一工作面水量随推采过程变化特征的是“动—静”储量法，它认为顶板来水由两部分组成，一部分是工作面开采形成的顶板大井空间内的静储量，另一部分是来自于大井外围含水层的动态补给量[14]，这种计算方法在实际生产中得到了广泛应用，但是，它没有考虑两种不同形式的储量随工作面走向长度和采动时间的相互转化，造成预测的工作面采空区涌水量是一个随推采不断增大的过程。实际上，大多数工作面采空区涌水量是一个随推采过程先升高，达到峰值水量，再逐步衰减的过程。众多学者[15-17]普遍认为，峰值水量一般出现在工作面初采初放位置，但还没有找到一种可以描述这一过程的模型，这种能够正确描述工作面涌水量随采掘工程的变化规律的模型是矿井防治水工作所急需的。

此次研究中，笔者从采动涌水过程本身的动态性特征研究出发，进行问题概化，基于系统动力学基本理论，构建能够描述这一过程的动力学模型，选取两个具有典型连续的采空区涌水观测数据工作面，研究工作面采动顶板涌水随推采过程变化特征，建立顶板涌水与采动进尺关系的动态模型。

1 典型工作面选取及采空区涌水过程动力学模型构建

1.1 问题假设及概化

工作面采动过程中，顶板水分别以静态重力释水和动态补给两种形式进入采煤空间，两种形式的水分别被称作静储量和动储量[14]，其中静态释水主要是伴随顶板周期性垮落的瞬间涌入过程。在以往研究中，认为周期性涌入采空区的这部分静态水量，是顶板上覆大井影响空间的增量，并且认为上覆影响范围含水层空间静储量在垮落之前未经提前释放[14]。实际上，在第n+1次顶板垮落，静储量涌入过程中，其上覆静储量已经在前n次垮落过程中被一定程度的逐段疏放。因此，伴随周期性垮落实际进入采空区的静储量随推采过程在不断减少，指示由开切眼向后垮落点上覆含水层水头是一个不断下降的过程，最终达到水头在垂向上的下降与水平推采等速时，周期性涌入采空区的静态水量达到动态稳定(图1)。

工作面推采过程，实际上是采煤空间形成的大井半径不断增加的过程(图2)，根据吴吉春等研究结果，利用Dupuit公式计算井流量有如下现象:对于同样的岩层和降深，小井径状态下(一般<200 mm)，井径增加所引起的流量增长率大，中等井径时增长率减小，大井径时(>500 mm)流量随井径没有表现出明显的增加[18]。对于这种现象，有学者认为是由于井周围的紊流和三维流影响所引起的，或者是与含水层透水能力存在一定的关系，这是稳定流假设条件下的井流量随井径变化的结果。实际上，采煤过程井径的动态增加导致地下水运动形式必定是非稳定流，但目前，用来计算工作面动态补给量的方法大多采用的是地下水向井稳定运动的Dupuit承压—无压公式，这与实际是存在差异的。

根据众多工作面采空区实际涌水量观测结果:推采前期，随采空区大井井径不断增加，动态补给逐步增加。推采后期，有限延展含水层补给条件下采空区动态补给增加水量逐渐减少，或基本趋于稳定。可近似于无限延展含水层补给条件下，随井径不断增大，采空区动态水量将会一定程度的增加，但这种增量有限，按照解析法计算，动态增量与伴随顶板垮落瞬间进入采空区的静储量相比，特别是与基本顶初次完全垮落产生的瞬时静储量相比，存在较大差距，而且随推采大井井径不断增加，影响半径不断向外发展至隔水边界，原本可近似假定为无限延展含水层补给的条件也不复存在，动态补给水量将不再增加。

1.2 工作面推采采空区涌水过程动力学模型构建

工作面推采顶板静储量释放过程可以用单指数衰减模型(一阶动力系统)刻画，它是基于一系列离散模型代替并模拟系统动态连续变化过程的，也称作系统的弛豫过程。多用于土壤养分流失[19-21]、有机荧光分子寿命测量[22]以及放射性元素能量衰减动态变化等过程模拟。

工作面静储量Qj随推采释放过程模型表述如下

Qj=A1exp(-k1at/b)+C1

式中，at/b表示推采t时刻已回采工作面走向长度与倾向宽度的比值(简称“走倾比”);Qj为at/b位置顶板静储量释放水量，m3/h;k1为静储量释水衰减系数，m3/h·L;C1为顶板静储量释放水量动态稳定值(Qje),m3/h，根据模型意义，容易推导A1=Qj0-Qje;Qj0为推采直接顶初次垮落时刻位置释放静储量[23]，可表述为Qj0=bvhc;v为直接顶初次垮落位置推采速度，m/h。

同样，工作面采动动态补给涌水随推采过程的变化可用一阶指数恢复系统描述，它也是一种系统最终归于平衡态的弛豫过程，核磁共振中质子受到磁场作用后的纵向磁化强度向平衡态的恢复过程就是典型的系统恢复弛豫过程[24]。

工作面动态补给随推采逐步增加的过程模型为

Qd=C2-A2exp(-k2at/b)

式中，Qd为at/b位置顶板动态补给水量，m3/h;k2为动态补给水量增长系数，m3/h·L;C2为顶板动态补给水量达到的动态稳定值(Qde),m3/h，根据模型意义，初始时刻至直接顶初次垮落时刻位置，顶板等效于无破坏，动态补给水量为0，at/b趋于无穷大时，Qd=C2，据此，容易推导A2>C2。

利用系统弛豫过程基本原理刻画了工作面采空区两种涌水形式随推采过程的模型，据此，工作面采空区涌水量随推采过程的一般变化规律如图3所示。

采空区涌水过程可用双指数衰减模型表述为

Q=Qj+Qd=A1exp(-k1at/b)-

A2exp(-k2at/b)+C

式中，C=C1+C2。

需要进一步说明的是，一般情况下，系统弛豫过程，状态变量是时间要素，但具体到本文研究的问题，用走向推采长度与倾向长度的比值(走倾比)作为状态变量，其意义远比时间要素大，这主要有两点原因:首先，煤矿防治水领域对于工作面采空区涌水最为关心的是，正常推采条件下的峰值涌水量的发生位置，而不是发生时间;其次，工作面推采速度并不是匀速推进，利用模型找到的峰值水量发生时间以及由此计算出来的推采长度不具备统计学意义，对于工作面采动涌水规律研究意义不大。

1.3 典型工作面案例选取

宁夏宁东煤田鸳鸯湖矿区位于鄂尔多斯盆地西缘，是典型的大水矿区，主要开采的侏罗系延安组2煤面临上覆直罗组下段巨厚强富水粗砂岩含水层威胁[25]，该含水层平均厚度(M)近100 m，渗透系数(K)介于0.001～0.3 m/d，大部分地段单位涌水量(q)在0.1 L/(s·m)以上，矿井水文地质条件复杂—极复杂。

选取的区内典型工作面(A和B)水平距离较大的断层构造在5 km以上，回采前均进行了一定程度的顶板水预疏放，采前工作面顶板含水层残余水头约15 m，工作面回采初期涌水过程可近似假定为无限延展含水层补给条件下，地下水向采空区大井的非稳定运动;回采后期，随大井井径不断扩展至水文地质边界(隔水断层)，补给水量将不再增加。A,B工作面走向长度分别为1 603.9和4 118.7 m，倾向宽度分别为305和210 m，采厚分别为5.09和2.80 m，煤层顶板距离直罗组下段粗砂岩含水层平均分别为4.95和22.20 m，相对稳定，经验公式计算导水裂缝带发育高度约分别为100和56 m。由此可知，两个工作面开采，导水裂缝带均能波及上覆直罗组下段含水层。同时，选择的两个工作面采动时间相近。A,B工作面推采过程中采空区涌水量观测频率分别为15和10 d(图4)。

2 模型求解及讨论

2.1 工作面推采采空区涌水量动态特征

A,B工作面涌水量随推采过程均表现出很好的双指数变化规律(图5)，利用Matlab 2016a软件拟合得到的随机曲线决定系数分别高达0.89和0.92，涌水量弛豫过程达到动态平衡的水量分别为54.92和38.62 m3/h，模型中其他参数值(A1,A2,k1,k2)并不惟一(表1)，主要是因为以上4个模型参数是相互关联的，其中一个参量的变化，可以通过其他参量的变化作补偿，而不改变拟合曲线结果，也就是说，在不考虑其代表的实际物理意义时，参数值可以在很大的空间范围内选取，而曲线形状及拟合效果基本没有变化[26]，但由于模型参数在本文中是有确切物理意义的，其值应是定值。实际上，双指数模型参数求解问题是较为困难的[27-28]，多个领域均在尝试利用不同的方法，以期得到更为逼近实际的值，例如，陈文驹等[22]利用迭代卷积法和仪器函数求解了双指数荧光衰减模型参数，林久令等[28]基于步长加速法完成了多指数模型拟合和参数的求解。因此，采空区涌水过程模型其他参数的求解是一个难点。

2.2 基于麦克劳林公式的模型参数求解

针对2.1节提出的采空区涌水双指数衰减过程面临的求参问题，尝试利用麦克劳林公式指数多项式在零点化解为非线性多项式近似求解的基本思想，求解采空区涌水过程双指数模型参数。模型化解过程如下:根据A,B工作面实际尺寸，模型自变量at/b最大为5.26和19.21，为了使自变量满足麦克劳林公式限定的自变量趋近零点的要求，笔者另取自变量x，令x=(at/b)/1 000，这样A和B工作面采空区涌水双指数模型可变为

注:Qe为动态平衡水量，Qm为峰值水量。

QA=A11exp(-k111 000x1)-

A21exp(-k211 000x1)+54.92

(1)

QB=A12exp(-k121 000x2)-

A22exp(-k221 000x2)+38.62

(2)

在x=0附近，方程(1)和(2)可用麦克劳林公式化为

/n!+Rn(0)

(3)

/n!+Rn(0)

(4)

实际上，按照麦克劳林公式展开时，3次项以后的余项值对涌水量Q值的影响可以忽略，但利用Excel软件进行多项式拟合时，3次多项式拟合效果并不理想(A工作面拟合决定系数仅为0.34)，利用麦克劳林公式和多项式对应项相等求参，难以得到较为理想的结果。但利用6次多项式进行拟合，A,B工作面采空区涌水量拟合曲线决定系数均在0.9以上，效果较好(图6)。

根据参数表达的物理意义，利用实测数据求解了A1值(表2)。再根据6次多项式拟合得到模型各项系数与麦克劳林公式对应项系数相等的关系，基于最小二乘法原理，求解了模型参数A2,k1,k2的最佳近似解，并据此求解出工作面推采过程峰值涌水量及发生位置，动态平衡水量及发生位置(表3)。

2.3 关键参数讨论

工作面采动进入采空区的动、静储量互相转化，此消彼长，A,B工作面静储量释水速率(k1)明显小于动态补给增长速率(k2)，并且A1值小于A2值(表3)，表示采煤扰动顶板进入采空区的水量变化曲线形态主要受控于动态补给量。从两个工作面曲线变化过程来看，涌水量峰值以及动态平衡值出现位置并没有表现出很好的规律性，峰值水量出现位置分别在走倾比约为0.80和3.49的位置，这与部分学者研究成果并不一致[15-17]。一般认为，峰值涌水量发生位置一般在工作面推采至走倾比约为1的位置，实际上，在上覆含水层厚度、富水性和导水性均一且构造条件影响可以忽略的条件下(按大井法推算，走倾比0.8大井影响半径尚未波及断层边界)，峰值水量出现位置主要受控于推采速度和顶板导水裂缝带发育高度，从A,B工作面计算结果来看，B工作面平均推采速度是A工作面的2.5倍，但峰值水量发生位置B较A明显滞后，指示推采速度与峰值水量发生位置成正比，导水裂缝带发育高度主要受控于顶板岩性条件[29]，同时，又决定顶板含水层波及厚度[30]。因此，实际涌水过程峰值水量及相应的出现位置是地质、水文地质以及采动条件综合作用的结果。

表2 采空区涌水过程双指数模型参数A1计算结果
Table 2 Results of A1 calculated using the measured data

表3 采空区涌水过程双指数模型最终计算结果
Table 3 Results calculated using the McLaughlin formula

动态平衡水量及发生位置的研究对于工作面采动后期防排水设计具有重要意义，与峰值水量发生位置类似，推采速度较快的B工作面较A工作面需要的推采稳定距离更长。但最终的稳定水量却出现相反的现象，这主要与煤层开采顶板覆岩破坏波及的含水层厚度有关，尽管采前两个工作面上覆残余水头基本一致(约15 m)，但采动后期，工作面水量主要受控于动态补给量，A,B工作面上覆含水层基本特征基本一致，但A工作面采厚为B工作面的1.82倍，且A工作面上覆隔水层厚度较B小约17 m，经验公式计算的A工作面采动顶板导水裂缝带波及含水层厚度(95 m)仍远大于B工作面(34 m)，意味着采动后期A工作面采空区大井接受的垂向补给范围较B工作面更大。

需要说明的是，利用上述方法最终计算的模型参数也仅仅是相对逼近于实际值的数字，仍存在一定程度的误差。最终形成的曲线方程，实现了对工作面推采全过程采空区涌水量的拟合，拟合实际观测数据效果较好，拟合曲线均方差分别为9.55和23.39，决定系数分别为0.88和0.90(图7)。与初始随机拟合曲线相比，最终构建的模型曲线拟合精度虽略有下降，但更能够从全局角度准确表达推采涌水过程的物理意义。

从典型水量拐点位置来看，峰值水量和动态平衡水量计算误差相对较小(表4);而大量统计资料表明，以往计算方法对峰值水量或者平均水量计算误差小于30%的矿井也仅有10%，误差远大于双指数模型计算结果。从全局角度来看，双指数模型计算误差在10%左右(表4)，而以往计算方法针对工作面推采全过程涌水量变化趋势预测与实际存在的偏差远大于这一数字[3]。对比结果指示双指数预测模型计算精度更高。

表4 工作面采空区涌水过程模型计算相对误差
Table 4 Relative error of water-inflow calculated using the McLaughlin formula

注:全局误差/n;Oi表示第i次涌水量观测值;Ci表示对应模型拟合值;—表示缺乏观测值。

3 结 论

(1)建立的宁东鸳鸯湖矿区A,B工作面采动顶板涌水双指数衰减模型，拟合决定系数分别达到0.88和0.90，表示该模型能够用来描述煤矿工作面顶板动、静两种形式的水随推采叠加进入采空区的变化过程。

(2)利用麦克劳林公式等方法求解了A,B工作面采空区涌水双指数衰减模型参数，确定了A,B工作面采动过程峰值水量及发生位置、动态平衡水量及发生位置，与实测值进行的误差分析结果，指示模型计算精度较高。

(3)峰值水量和动态平衡水量发生位置主要受控于推采速度;峰值水量和动态平衡水量受采厚、顶板岩性特征、含水层特征等条件共同影响。

致谢 感谢中煤科工集团西安研究院有限公司李德彬、张雁博士提供的大量基础数据，感谢西北农林科技大学陈实副教授在模型构建过程中给予的大量建议，感谢中煤科工集团西安研究院有限公司乔伟、尚宏波，兰州大学董军、朱高峰教授以及华中科技大学杜远博博士在模型求解过程中给予的帮助，还要感谢现场工作人员在井下基础数据记录中付出的辛苦劳动！

参考文献(References) :

[1] 虎维岳,田干.我国煤矿水害类型及其防治对策[J].煤炭科学技术,2010,38(1):92-96.

HU Weiyue,TIAN Gan.Mine water disaster type and prevention and control countermeasures in China[J].Coal Science & Technology,2010,38(1):92-96.

[2] 虎维岳,闫丽.对矿井涌水量预测问题的分析与思考[J].煤炭科学技术,2016,44(1):13-18.

HU Weiyue,YAN Li.Analysis and consideration on prediction problems of mine water inflow volume[J].Coal Science & Technology,2016,44(1):13-18.

[3] 虎维岳.浅埋煤层回采中顶板含水层涌水量的时空动态预测技术[J].煤田地质与勘探,2016,44(5):91-96.

HU Weiyue.Water inflows prediction technique of water inflow from roof aquifer during extraction of shallow seam[J].Coal Geology & Exploration,2016,44(5):91-96.

[4] 薛禹群.地下水动力学原理[M].北京:地质出版社,1986.

[5] 陈崇希,唐仲华.地下水流动问题数值解法[M].武汉:中国地质大学出版社,1990.

[6] 虎维岳.矿井水害防治理论与方法[M].北京:煤炭工业出版社,2005.

[7] 周如禄,戴振学.矿井涌水量预测的理论与实践[J].煤炭科学技术,1998,26(6):47-49.

ZHOU Rulu,DAI Zhenxue.The theory and practice of coal mine inflow prediction[J].Coal Science & Technology,1998,26(6):47-49.

[8] 骆祖江,杨林,周强,等.煤田矿井涌水量预测三维数值模型[J].水动力学研究与进展,2010,25(6):737-742.

LUO Zujiang,YANG Lin,ZHOU Qiang,et al.Three-dimensional numerical model for predicting groundwater inflow of mine in coal field[J].Chinese Journal of Hydrodynamics,2010,25(6):737-742.

[9] 骆祖江,王琰,陆顺,等.基于矿井生产过程的涌水量预测三维数值模拟模型[J].煤炭学报,2010,35(S1):145-149.

LUO Zujiang,WANG Yan,LU Shun,et al.Three-dimensional numerical simulation model for predicting groundwater inflow of mine[J].Journal of China Coal Society,2010,35:145-149.

[10] 杜敏铭,邓英尔,许模,等.矿井涌水量预测方法综述[J].四川地质学报,2009,29(1):70-73.

DU Minming,DENG Yinger,XU Mo,et al.Review of methodology for prediction of water yield of mine[J].Acta Geologica Sinaica,2009,29(1):70-73.

[11] 连会青,夏向学,徐斌,等.矿井涌水量预测方法及适用性评价[J].华北科技学院学报,2014,11(2):22-27.

LIAN Huiqing,XIA Xiangxue,XU Bin,et al.Evaluation and applicability study on prediction methods of water inflow in mines[J].Journal of North China Institute of Science and Technology,2014,11(2):22-27.

[12] POLUBARINOVA-KOCHINA P Y.Theory of ground water movement[M].Princeton:Princeton University Press,1962.

[13] 乐建,李文平,乔伟,等.采动影响下多工作面顶板涌水量预计方法[J].矿业安全与环保,2009,36(3):22-24.

LE Jian,LI Wenping,QIAO Wei.Prediction method of roof water inflow on multi-working faces under mining effect[J].Mining Safety & Environmental Protection,2009,36(3):22-24.

[14] 刘洋,张幼振.浅埋煤层工作面涌水量预测方法研究[J].采矿与安全工程学报,2010,27(1):120-124.

LIU Yang,ZHANG Youzhen.Forecast method for water inflow from working face in shallowly buried coal seam[J].Journal of Mining & Safety Engineering,2010,27(1):120-124.

[15] 朱宏军.鸳鸯湖矿区矿井涌水量预测方法研究[D].北京:煤炭科学研究总院,2014.

ZHU Hongjun.The prediction method of mine water in Yuanyanghu mining area[D].Beijing:China Coal Research Institute,2014.

[16] 郑纲,何渊.深埋矿井超大工作面顶板砂岩含水层涌水量预测方法[J].煤矿安全,2015,46(S1):76-80.

ZHENG Gang,HE Yuan.Prediction of roof sandstone water discharge at super large working face in deep mine[J].Safety in Coal Mines,2015,46(S1):76-80.

[17] 刘洋,王振荣,牛建立.工作面涌水量预测方法的确定[J].矿业安全与环保,2010,37(5):29-30.

LIU Yang,WANG Zhenrong,NIU Jianli.Determination of prediction method of water inrush from working face[J].Mining Safety & Environmental Protection,2010,37(5):29-30.

[18] 吴吉春,薛禹群.地下水动力学[M].北京:中国水利水电出版社,2009:67-68.

[19] WOODRUFF C M.Estimating the nitrogen delivery of soil from the organic matter determination as reflected by Sanborn Field[J].Proceedings Soil Science Society of America,1949(14):208-212.

[20] WEI Guoxiao,ZHOU Zhenfang,GUO Ying,et al.Long-term effects of tillage on soil aggregates and the distribution of soil organic carbon,total nitrogen,and other nutrients in aggregates on the Semi-Arid Loess Plateau,China[J].Arid Land Research and Management,2014(28):291-310.

[21] YANG Lixia,PAN Jianjun,YUAN Shaofeng.Predicting dynamics of soil organic carbon mineralization with a double exponential model in different forest belts of China[J].Journal of Forestry Research,2006,17(1):39-43.

[22] 陈文驹,林美荣,姜宏丽.双指数荧光衰减动力学模型分析[J].光学学报,1966,6(12):1124-1129.

CHEN Wenju,LIN Meirong,JIANG Hongli.Analysis of Kinetic model of double exponential fluorencence decay[J].Acta Optica Sinica,1966,6(12):1124-1129.

[23] 王庆康,霍振奇.采场直接顶的初次垮落规律及其计算方法[J].采矿与安全工程学报,1987(1):9-15,66.

WANG Qingkang,HUO Zhenqi.The first caving law of stope direct roof and its calculation method[J].Journal of Mining & Safety Engineering,1987(1):9-15,66

[24] 熊国欣,李立本.核磁共振成像原理[M].北京:科学出版社,2007.

[25] 赵宝峰.基于AHP-物元可拓模型的顶板水疏放程度评价[J].煤田地质与勘探,2017,45(3):82-86.

ZHAO Baofeng.Evaluation on drainage degree in roof based on AHP-matter element extension model[J].Coal Geology & Exploration,2017,45(3):82-86.

[26] 俞文 ,黄晓.多指数衰减曲线的拟合问题[J].应用数学与计算数学学报,1990,4(1):57-69.

YU Wenji,HUANG Xiao.Curve fitting problems of multi-exponential decay expressions[J].Comm.on Appl.Math.and Comput.,1990,4(1):57-69.

[27] 郭纲,王树勋,孙晓颖,等.超声信号的双指数模型及参数确定方法[J].电子学报,2009,37(7):1501-1504.

GUO Gang,WANG Shuxun,SUN Xiaoying,et al.Ultrasonic double exponential model and its parameter determination[J].Acta Electronica Sinica,2009,37(7):1501-1504.

[28] 林久令,董志山,张秀峰,等.荧光衰减曲线的一种多指数拟合方法[J].发光学报,2002,23(1):93-96.

LIN Jiuling,DONG Zhishan,ZHANG Xiufeng,et al.A method for multi-exponential fitting to the fluorescence decay[J].Chinese Journal of Luminescence,2002,23(1):93-96.

[29] 胡小娟,李文平,曹丁涛,等.综采导水裂隙带多因素影响指标研究与高度预计[J].煤炭学报,2012,37(4):613-620.

HU Xiaojuan,LI Wenping,CAO Dingtao,et al.Index of multiple factors and expected height of fully mechanized water flowing fractured zone[J].Journal of China Coal Society,2012,37(4):613-620.

[30] 刘英锋,郭小铭.导水裂缝带部分波及顶板含水层条件下涌水量预测[J].煤田地质与勘探,2016,44(5):97-101.

LIU Yingfeng,GUO Xiaoming.Prediction of water inflow in roof aquifer affected by water-flowing fracture zone[J].Coal Geology & Exploration,2016,44(5):97-101.