赵瑞珏，毛德强，刘再斌，姬中奎，曹祖宝

（1.山东大学土建与水利学院，山东 济南 250061；2.中煤科工集团西安研究院有限公司，陕西 西安 710054；3.长安大学环境科学与工程学院，陕西 西安 710054）

在煤炭开采的过程中，矿井突水事故是影响安全生产的重要因素之一[1]。为了探明同一含水层内部或不同含水层之间的水力联系，深入了解含水层富水性以及补给、径流、排泄条件，通常进行同一矿区大流量、大降深、大范围的放水试验[2-3]。查清本井田及相邻井田的水文地质条件，为煤矿安全生产提供技术保障[4]。

传统放水试验分析方法得到的含水层水文地质参数是含水层的等效参数，不能准确刻画参数的空间异质性。邵红旗等[5]提出了以水量衰减率、单位涌水量的水位降深值及水位恢复速率为主要指标的放水试验分析方法；曹祖宝等[6]利用“大井法”分析放水试验，得到的水文地质参数是含水层的等效参数。

水力层析法是基于抽水或注水试验的改进方法，通过一系列交叉抽水或注水试验，运用水头响应数据，识别含水层水文地质参数的空间分布特征[7]。在同样抽水井数量的情况下，水力层析法能够获得更多有关含水层的水头响应数据，分析结果能够更加准确地反映含水层水文地质参数的空间分布特征，求解结果更加接近实际。蒋立群等[8]基于室内非均质含水层砂箱实验，研究了水力层析法和克立金法刻画非均质含水层渗透系数场的优劣；郝永红等[9-10]利用数值算例系统分析了水力层析法可准确获取含水层的裂隙带分布信息；Illman 等[11]基于两次抽水试验数据，运用水力层析法建立了日本Mizunami 地下实验室水文地质参数的三维分布图，揭示了断层导水特征；董艳辉等[12]应用水力层析法研究了离散裂隙网络的连通性，获取的渗透特性分布可以更好地预测渗流规律；Fischer 等[13]通过结构约束模型推断了岩溶地区水力特性的空间分布。但目前为止，还未有利用煤矿放水试验获取含水层水文地质参数分布信息的相关研究。

在本研究中，运用水力层析法分析煤矿井下叠加放水试验数据，反演计算含水层水文地质参数，并根据互相关分析具体解释水力层析法识别含水层参数空间分布特征的原理。该方法操作简单，适用性广，获取的参数详细具体，可准确刻画含水层的非均质性，精确定位富水异常区的位置及范围，有效降低突水事故风险，为煤矿放水试验数据分析及矿井水文地质条件评价提供一套新方法。

1 水力层析法原理

同步连续线性估计算法（Simultaneous Successive Linear Estimator，简称SimSLE）发展成熟，且已得到广泛应用[14-16]。本研究基于SimSLE 算法，结合煤矿井下叠加放水试验，具体解释水力层析法的基本原理。

1.1 地下水流控制方程

地下水流控制方程用于模拟地下水流场变化，是进行水文地质参数反演的基础。对于饱和带非均质各向同性介质中的非稳定承压水问题，地下水流控制方程为：

边界条件和初始条件分别为：

式中：H—水头/m；

x—二维空间坐标/{m,m}；

Q（xp,t）—xp处t时刻的放水孔流量/（m·s-1）；

T（x）—导水系数/（m2·s-1）；

K（x）—渗透系数/（m·s-1）；

S（x）—储水系数；

H1—Dirichlet 边界Γ1的水头/m；

q—Neumann 边界Γ2单位法向量n方向的流量/（m·s-1）；

H0—含水层初始水头/m。

1.2 同步连续线性估计算法

SimSLE是一种基于地质统计学的参数估计方法[17]。不同于传统放水试验分析中的均质性假设，SimSLE将导水系数和储水系数的自然对数及观测孔水头处理成空间随机变量：

f'、s'、h'—扰动值。

SimSLE 从带有权重的线性估计开始，利用直接测量的含水层水文地质参数值和放水试验过程中的水头测量值，计算得到参数的初始估计值。权重根据水头测量值与水文地质参数的统计属性，如均值、相关长度、水头测量值和参数之间的互协方差和自协方差，利用协同克里格（cokriging）方法计算得到。在估计过程中，参数的自协方差不断更新，代表估计不确定性的变化[18-19]。以导水系数为例，初始估计结果为：

H*—实际水头测量值；

He—通过地下水流控制方程求得的水头模拟值；

Nf—参数测量点的总数；

Np—放水试验的总次数；

Nh（k）—第k次放水试验中水头测量点的总数；

Nt（k,j）—第k次放水试验中第j个水头测量点随时间变化的水头测量数据的总数；

λ0i—协同克里格权重，即第i个参数测量点的参数值f*对x0处的参数估计值的权重；

μ0kjl—实际水头测量值H*（k,xj,tl）对参数估计值的权重。

将参数估计值代入地下水流控制方程，计算在水头测量点不同时刻的模拟值，判断模拟值与实际测量值之间的差异，若差异没有达到指定精度，上述计算过程会重复迭代进行。

1.3 互相关分析

运用水力层析法时，不同位置的放水孔和观测孔水头响应数据反映了含水层不同区域的水文地质参数信息。通过分析随时间变化的水头测量值与导水系数和储水系数的互相关性，可以进一步解释水力层析法识别含水层空间异质性的原理[20]。处的一阶泰勒公式为扰动项的表达式为：

式中：J—雅可比（Jacobian）系数。

假设导水系数和储水系数相互独立，式（9）两边同时右乘f′，并求期望，且利用矩阵形式表达：

式中：Rhf—h与f的互协方差矩阵；

Rff—f的自协方差矩阵。含水层参数的空间变异性通过自协方差矩阵中的方差和相关长度表示；

Rhs—h与s的互协方差矩阵；

Rss—s的自协方差矩阵。

h′

式（9）两边同时右乘，并求期望，代入式（10）、（11）：

式中：Rhh—h的自协方差矩阵。

由相关系数的定义，代入式（10）、（11）：

任意位置的水头H（x,t）在均值

式中：ρ—互相关系数矩阵；

σ2h、σ2f、σ2s—h、f、s的方差，即自协方差矩阵Rhh、Rff、Rss的对角线元素。

互相关系数矩阵可表征空间任意位置水文地质参数的不确定性对观测孔水头测量值的影响。互相关分析是在随机背景下进行的敏感性分析，采用了地质统计学的概念，通过期望和方差代入水文地质参数的先验信息（prior information），综合分析水文地质参数的空间异质性对水头的影响。

互相关系数的取值范围为[-1,1]。互相关系数大于零，即为正相关，表现为当空间某一位置水文地质参数增大时，水头测量值相应增大；反之，互相关系数小于零，即为负相关，表现为水文地质参数增大时，水头测量值相应减小。互相关系数的绝对值越大，相关性越强，表明若赋予水文地质参数一个微小的扰动值，水头测量值相应地改变较大，说明对于该处水文地质参数的变化比较敏感，该观测孔的水头响应数据对于反演水文地质参数效果显著。

2 研究区概况

研究区柠条塔井田位于陕西省榆林市境内。煤系地层地表露头近乎水平展布，下伏煤层底板向西倾斜，坡降为4‰～10‰。在以往勘探及矿井建设和到目前为止的生产过程中，均未发现较大的断层（落差大于15 m）或褶皱，仅煤层底板有一些宽缓的波状起伏[21-22]。

据地质填图及钻孔揭露，柠条塔井田地层由老至新（自下而上）分别为：侏罗系中统延安组（J2y）、直罗组（J2z），新近系上新统保德组（N2b），第四系中更新统离 石组（Q2l），第 四系全新统风 积层（Q4el）和冲 积层（Q4al）。延安组为主要含煤地层；直罗组为砂岩含水层；保德组为相对隔水层，结构较致密，渗透性较弱。开采延安组2-2煤层的S1210 工作面涌水量接近1 000 m3/h，维持了三个多月，涌水来源判断为直罗组砂岩含水层。涌水期间，直罗组上覆和下覆地层的观测孔未发生相应的水头变化，因此本研究将直罗组含水层概化为二维含水层。

3 放水试验与数值建模

以柠条塔煤矿水文地质勘探过程中的井下叠加放水试验为研究内容，依据水力层析技术，利用基于地质统计学的同步连续线性估计算法SimSLE，对水头数据进行定量解译，并对侏罗系中统直罗组（J2z）的含水层结构进行模拟计算。

3.1 井下叠加放水试验

在井田南翼较大范围内，工作面推采前方的水文地质条件尚未完全探明，需进一步查明是否存在富水异常区。工作面井下施工钻孔数量较多，自切眼往外约2 000 m 范围内工作面顶板有几十个钻孔（图1），斜向上全部贯穿直罗组含水层直至保德组黏土层。利用井下钻孔进行采前大流量放水试验。第四系地面观测孔（图1Q5）测量放水试验过程中离石组含水层的水头变化，探查工作面涌水与第四系水的水力联系，测量数据表明放水过程中没有明显的水头变化，故进一步确定保德组为相对隔水层，放水试验与第四系水无水力联系。

放水试验共持续9 d，放水次序为自南向北，采用叠加放水的方式。根据放水前对井下所有钻孔涌水量的排查情况，共选取工作面顶板上30个水量较大的钻孔作为放水孔。分区域选取28个井下钻孔（其中工作面顶板钻孔17个，巷道外侧钻孔11个）和6个地面钻孔作为水头观测孔。

图1 钻孔分区布置平面图Fig.1 Plan of the boreholes zoning scheme

34个观测孔全程测量水头。试验正式开始时，第一阶段开启F3-1、F5、F5-1、F6、F6-1、F7 共6个放水孔阀门（图1）。第一阶段放水持续2 d。第二阶段开启T6、F8、F9、F10、F11、F12、F12-2 共7个放水孔阀门（图1）。此时共有13个钻孔联合放水。第二阶段放水持续2 d。第三阶段开启T9、F15、F16、F17-2、F19、F19-1、F20-1、F21-2 共8个放水孔阀门（图1）。此时共有21个钻孔联合放水。第三阶段放水持续3 d。第四阶段开启YZ10-2、YZ11-1、F22、F22-1、F24、F26-1、F26-2、F28、F29 共9个放水孔阀门（图1）。此时共有30个钻孔联合放水。第四阶段放水持续2 d。每个阶段的稳定放水量分别为270，214，167，96 m3/h，最后稳定放水总量为747 m3/h。按照1 h时间步长，将采集到的流量数据分段累加，并保持水量不变，换算至整小时时刻处。

在放水的四个阶段分别选取一个井下观测孔，及一个地面观测孔，共5个观测孔表示放水试验期间的降深随时间变化的动态趋势（图2）。位于第一阶段放水孔群附近的观测孔SK14-2 在放水的第一阶段便有水头响应，降深随着放水的进行不断增大；位于第二阶段放水孔群附近的观测孔F9-1，由于第一阶段放水影响的含水层范围有限，故在第一阶段该处水头响应较弱，降深较小，第二阶段初期降深迅速增大，之后稳步上升；第三阶段放水孔群附近的观测孔F18的水头响应降深见图2；第四阶段放水孔群附近的观测孔YZ10-1的水头响应降深见图2。在整个放水试验过程中，直罗组地面观测孔J12、J17、SK6、SK8、SK11、SYJ的降深较小，均小于1 m，J17的降深见图2。

图2 降深随时间变化的动态图Fig.2 Representative temporal drawdown data

3.2 模型建立

概化数值模型为1 800 m × 3 200 m的二维非均质饱和含水层，将含水层模型划分为90 × 80个单元格，每个单元格为20 m × 40 m，四周为1 000 m的定水头边界，初始水头为1 000 m。定水头边界及初始水头设置为1 000 m可保证含水层始终维持饱和。水力层析反演计算时需给定水文地质参数的初始场，设定lnT=1.3，方差为1.0，lnS=-5.0，方差为0.5，x和y方向的相关长度均为100 m。利用放水试验获得的每小时观测孔水头数据，迭代估计含水层水文地质参数的空间分布特征。选取的18个时刻的水头数据分别位于t=1，6，20，40，50，56，70，80，97，100，110，130，150，169，174，185，195，216 h，保留了水头的整体趋势和关键转折点[23-24]。水力层析法的反演算法SimSLE 由亚利桑那大学Yeh 教授课题组开发的软件VSAFT2 实现，下载链接为http://tian.hwr.arizona.edu/downloads。

4 结果与分析

4.1 水力层析反演结果分析

水文地质参数的反演结果见图3。图3（a）、（b）是仅由第一阶段的观测孔水头数据反演得到的导水系数和储水系数的空间分布，一个靠近涌水位置的高T高S区初步显现；图3（c）、（d）是由第一、二阶段的水头数据联合反演得到的参数空间分布；图3（e）、（f）是由前三个阶段的水头数据联合反演得到的参数空间分布；图3（g）、（h）是由完整的四个阶段的水头数据联合反演得到的参数空间分布。

随着放水孔阀门的逐阶段开启，T和S的分布模式逐渐形成且轮廓清晰。工作面的南部涵盖强导水区域，且有可能延伸至模型区域之外，与大流量、长时间的涌水密切相关。在巷道内侧，强导水区域向北延伸，继而出现相对弱导水区域，并与弱储水区域对应。

实施第四阶段放水试验的目的是刻画工作面北部甚至更远区域参数的空间分布。由于放水试验影响范围有限，且钻孔集中分布于巷道附近，估计参数的空间异质性仅体现在巷道内侧有限区域内，巷道外侧大部分区域基本按照参数均值均匀分布。

图3中标注了T=2 m2/d、T=8 m2/d、S=0.005和S=0.011时的等值线，等值线内的区域分别对应相对低T区、高T区、低S区和高S区。由图3可推知，工作面南北富水性差异大，涌水区位于强导水带上；从涌水区往北，导水系数和储水系数整体逐渐减小；南部导水系数和储水系数均较大，属于富水异常区。据文献[25]所述，工作面南北富水性差异大与直罗河冲刷沉积有关，工作面位于冲刷带的尖灭部位，因此南北富水性差异大。

4.2 水头与水文地质参数的互相关分析

为了直观地将互相关系数的变化规律展现出来，在放水试验的各阶段分别选取一个观测孔，计算了观测孔在各阶段初期、中期、末期的水头与导水系数和储水系数的互相关系数。

图3 水文地质参数的反演结果Fig.3 Inversion results of the hydrogeological parameters

图4 水头与导水系数的互相关系数ρhf时空分布Fig.4 Spatial and temporal distribution of cross-correlation between head and transmissivity

由图4（a）、（b）、（c）可知，在放水试验初期（t=10 h），直罗组地面观测孔SYJ 与放水孔群之间存在明显的负相关性，该负相关区域在整个放水阶段持续存在，但随着时间的延续，范围逐渐缩小，负相关性逐渐减弱。在放水试验中期（t=110 h），该观测孔的上游（根据图中流线所示方向判断）及放水孔群的上游存在正相关性，随着时间的延续，该区域范围逐渐扩展，最终连为一个整体，环绕在负相关区域周围，正相关性逐渐加强。

由图4（d）、（e）、（f）可知，在放水试验的第一阶段，井下观测孔F9-1水头未发生明显变化，故互相关分析的初期选取t=80 h。互相关系数空间分布的发展趋势与观测孔SYJ 类似，由初期的负相关发展为末期的正相关。

由图4（g）、（h）、（i）可知，在第三阶段放水试验初期（t=100 h），井下观测孔F14-2 下游及上游分别存在明显的负相关区域及正相关区域，这两个区域在第三、四阶段持续存在，但随着时间的延续，负相关区域范围逐渐缩小，并逐渐减弱；正相关区域范围逐渐扩大，且逐渐增强。互相关系数空间分布的发展趋势与从第一阶段开始的放水类似，但没有显现出放水孔群上游的正相关区域。在通常的水头与导水系数的互相关分析中，一般针对一个观测孔及一个放水孔，由于这里是放水孔群，观测孔与放水孔之间的负相关区域及放水孔上游的正相关分布区域较大，最终导致正相关区域仅出现在观测孔上游。

由图4（j）、（k）、（l）可知，井下观测孔W1 与观测孔F14-2的互相关系数空间分布的发展趋势类似。

由图4可知，水头与导水系数的互相关系数空间分布的发展趋势类似，在初期负相关性较强，末期正相关性较强，且相关区域范围扩大。故放水试验末期的水头数据更利于估计导水系数，且不同位置观测孔的相关区域不同，相近时间的互相关系数分布类似，因此可联合不同观测孔在不同时间的水头数据，求取含水层导水系数的高精度空间分布。

对于相同观测孔，在同一时间段的水头与储水系数的互相关系数矩阵的计算结果见图5。水头与储水系数的互相关性规律相对简单，整个放水过程均呈正相关性，且靠近观测孔附近的相关性较强。随着时间的延续，正相关性逐渐减弱。由图5可知，放水试验初期的水头数据对于估计储水系数比较有效。

图4、图5中的虚线为正、负相关性较强区域的轮廓线，直观表征了该观测孔的水头数据可精细刻画的参数空间范围。对于同一观测孔，图中标注了对应时刻的同条流线，以判断观测孔和放水孔群各自的上、下游区域。

综合水头与导水系数和储水系数的相关性分析结果，各个阶段的相关系数分布具有相似性，正相关和负相关区域按照一定的规律变化。本研究针对每个观测孔选择的18个水头数据，集中分布在每个阶段的初期和末期，可以更加有效的获取研究区的水文地质参数分布信息。

图5 水头与储水系数的互相关系数ρhs时空分布Fig.5 Spatial and temporal distribution of cross-correlation between head and storativity

5 结论

通过本研究可见，利用水力层析法分析井下叠加放水试验，是在传统放水试验分析基础上的创新，根据本次井下叠加放水试验，可以得出以下结论：

（1）水力层析法是一种有效的非均质含水层参数识别的新技术，相较于传统放水试验分析方法，它在识别含水层非均质特性方面具有优势。

（2）互相关分析结果表明每个观测孔数据与水文地质参数都存在相关性较强的区域，且初期及末期数据分别与储水系数和导水系数相关性较强，可增加观测数据，中期数据适当稀疏选取，将多个观测孔数据联合起来结合水力层析技术，可以较准确地反演出符合实际情况的参数分布信息。

（3）本研究将井下放水视为针对含水层的刺激源，联合多个观测孔的水头响应数据，反演研究区域水文地质参数的空间分布，捕捉到了研究区域的主要地质结构特征，刻画较精细。针对有突水风险的矿井，预先采用水力层析法分析矿井水文地质条件，推采时避开富水异常区，或者提前采取有针对性的防治水工程措施，可以有效降低突水事故风险。