基于改进组合弹簧模型的矿井地应力场计算方法

赵毅鑫, 赵良辰, 杨东辉, 边 华, 杨 哲, 宫智馨

(1.中国矿业大学(北京) 共伴生能源精准开采北京市重点实验室, 北京 100083; 2.中国矿业大学(北京) 煤矿灾害预防与处置应急管理部重点实验室, 北京 100083; 3.中国矿业大学(北京) 能源与矿业学院, 北京 100083; 4.山西大同大学 煤炭工程学院, 山西 大同 037003)

地应力场是影响岩层运移和原位岩体力学性质的重要因素, 对采矿、隧道开挖、油气开采等影响显著。地应力是采矿工程中覆岩运动的力源, 在构造应力强烈的地区进行采矿活动, 可能会引起冲击地压、岩爆等灾害, 影响矿山的安全生产。因此准确获取地应力场对于提高开采设计的科学性及保障矿山生产安全有着重要意义。

要获得区域地应力场需基于地应力测定获取部分点位的地应力值及方向。目前, 常用的地应力测量方法主要包括: 水压致裂法[1]、套孔应力解除法[2]、非弹性恢复法[3]和声发射法[4]等。其中, 声发射法通过钻取的岩芯在实验室开展测试即可获得岩芯所对应位置的三维地应力, 较其他测量方法具有测量深度大[5]、成本低和效率高等特点。声发射法通过识别岩石Kaiser效应点确定原岩应力状态,该方法由KANAGAWA等[6]于1977年首次提出。声发射法测定地应力有单轴测试与三轴测试两种方式。单轴测试方法较为简便, 但试样亦会在出现Kaiser效应点前发生破坏[7]; 而三轴测试方法能相对准确地测定地应力[8–9]。

然而, 地应力测量因测点数或样品数有限, 对矿井区域地应力场的描述仍需结合应力场计算来开展。目前常用的地应力场计算方法包括: 估算法、正演算法及反演算法[10]。估算法是基于地应力计算模型对地应力场进行计算。具有代表性的地应力计算模型包括金尼克模型[11]、Anderson修正模型[12]、黄荣樽模型[13]、组合弹簧模型[14]和葛洪魁模型[15]等。正演算法主要指通过调整边界荷载或位移, 最小化测点处计算值与实测值间的差异, 从而计算出所需区域地应力场的方法[16], 其主要包括边界载荷调整法和边界位移调整法等[17]。反演算法主要指通过假设模型边界荷载类型与分布, 建立边界荷载与测点应力间的回归关系, 通过回归计算, 调整待定系数使回归方程误差最小以获得最优的模型边界载荷, 从而获取所需区域地应力场的方法[18]。反演算法一般包括线性拟合法[19]、非线性拟合法[20]和神经网络法[21–22]等。

笔者以神东矿区布尔台矿42105工作面为工程背景, 结合对布尔台矿5个钻孔中89个地应力测点的Kaiser法实测结果, 提出一种基于改进组合弹簧模型的区域地应力场计算方法。该方法考虑了实测数据的变化规律和岩层物性对区域地应力计算中的控制作用, 提高了反演精度和可靠性。同时,将区域地应力反演方法融入FLAC3D软件建模, 对布尔台矿42105工作面开采过程中矿压显现规律进行了数值模拟并与实测数据进行对比, 进而证明了所提出方法的优越性。

1 研究区域与地应力测量

1.1 研究区域与取样

布尔台矿是神东煤炭集团主力生产矿井之一,位于内蒙古自治区鄂尔多斯市伊金霍洛旗乌兰木伦镇。大地构造分区属于华北地台鄂尔多斯台向斜伊盟隆起中东部, 伊陕斜坡东北部。基本构造形态为一南西倾斜的近水平产状的单斜构造, 井田构造简单, 层位稳定。

试验岩芯取自布尔台矿的白垩系及侏罗系地层, 取芯钻孔包括BK212, BK209, BK213, BK207和BK220, 共计5个钻孔。图1为现场取芯钻孔位置及岩层综合柱状图等。将取芯钻孔岩样沿垂直于其轴线方向再按逆时针0°, 45°和90°钻取尺寸为φ25 mm×50 mm的岩芯。

图1 布尔台矿区域位置及取芯情况Fig.1 Location of Buertai Mine in Shendong mining area

1.2 基于Kaiser效应的地应力测量

试验采用WDW–100E型试验机加载系统对岩样施加载荷, 采用美国物理声学公司(PAC)生产的PCI–II型AE监测系统监测声发射信号。WDW–100E试验机提供的最大轴向载荷量程为100 kN, 试验力测量范围为0.4%~100%FN, 变形测量范围为2%~100%FS。试验力测量精度、位移测量准确度及变形示值相对误差均为±1%以内。PCI–II型AE监测系统具有40 MHz的采样率和18位的A/D转换率, 其频带范围为1 kHz~3 MHz。AE数据由AE win for PCI–2软件进行采集、存储和分析。在试验过程中, AE监测系统的前置放大器增益设置为40 dB,阈值设定为45 dB, 采样率为1 MHz, 并采用双通道传感器进行数据采集, 谐振频率设定为20~400 kHz。使用双通道传感器采集数据, 传感器沿试样两侧轴向中部布置。试验时控制室内温度为26 ℃, 湿度为58RH%, 试验中监测应力、应变、AE计数等变量。循环加载模式为载荷控制与位移控制相结合的多次循环加载。前4次循环采用载荷控制, 末循环采用位移控制。其中前2次循环主要用于AE法分析, 之后的循环用于DRA法分析。首循环峰值载荷设定为σhmax的80%。之后3次循环既要减少首循环加载摩擦型AE的干扰, 又要避免岩石扩容损伤产生的影响, 其循环峰值载荷应大于σhmax, 且不超过岩石抗压强度σc的63%~70%。试验加载速率设置为0.05 kN/s, 卸载速率为0.25 kN/s。加载方式如图2所示。

图2 循环加载模式示意Fig.2 Schematic diagram of cyclic loading mode

采用AE–DRA法[23]识别Kaiser点。在岩石试件进行单轴压缩试验时, 记录声发射事件的数量、能量和频率等参数, 然后利用数据重构算法(DRA)对声发射信号进行分析, 找出声发射响应曲线的最大曲率点, 并将其作为Kaiser点。进一步根据Kaiser点对应的应力值计算地应力。

1.3 古地磁岩芯定向

由于钻取岩芯时采用非定向取芯方法, 故采用古地磁岩芯定向法确定岩芯在地下的原始方位。其基本原理为: 岩芯在形成时会记录当时的地磁场方向, 形成稳定的剩磁。通过实验室的高灵敏度磁力仪系统, 分离和测定岩芯的剩磁方向, 进而确定岩芯在地下的空间状态[24]。

试验使用美国2G超导磁力仪和英国MMDT80热退磁仪。在退磁后进行剩磁测量, 获得岩芯黏滞剩磁的平均磁偏角Da与磁倾角Ia, 结合标志线与古地磁方向, 恢复岩芯在地下的原始方位; 进一步计算出最大水平主应力方向。

布尔台矿区部分钻孔岩芯Kaiser法测定的地应力值及方向详见表1。

表1 布尔台矿区部分钻孔岩芯Kaiser法地应力测试结果Table 1 Stress measurement results in Buertai mine area based on Kaiser method

1.4 地应力测量结果验证

由表1可知, 同一钻孔中的地应力测点, 随着深度的增加, 最大水平主应力σH、最小水平主应力σh和竖直主应力σv均逐渐增加。使用最小二乘法线性拟合各点地应力数据, 结果如图3所示。可以发现: 竖直主应力拟合曲线斜率最大, 即随埋深增加竖直方向主应力增加最快; 埋深127 m以浅的地层,地应力呈σv<σh<σH的逆断型应力状态; 埋深1 725 m以深的地层, 地应力呈σh<σH<σv的正断型应力状态;埋深处于127~1 725 m的地层, 地应力呈σh<σv<σH的走滑应力状态。该结果与其他针对神东矿区构造应力场的研究结果相一致[25–26]。需要说明: 采用线性函数对水平方向地应力数据点进行拟合时, 离散性较大。因此, 尝试将岩石弹性模量E引入对地应力的线性回归拟合中, 进行多元线性回归。多元线性回归拟合结果如图4所示。

图3 布尔台矿区地应力随埋深变化规律Fig.3 In-situ stress variation law of the Buertai mine area with the change of depth

图4 布尔台矿区地应力随埋深与弹性模量变化规律Fig.4 In-situ stress variation law of the Buertai mine area with the change of depth and modulus of elasticity

由图4可知, 引入弹性模量E后, 拟合优度R2明显上升。对最大水平主应力的多元线性回归R2为0.612 01, 较线性回归的0.501 41提升22.06%; 最小水平主应力的多元线性回归R2为0.631 05, 较线性回归的0.506 27提升24.65%。结果表明, 引入弹性模量E对于改善水平方向地应力拟合效果具有显著作用。同时也说明水平方向地应力与弹性模量间具有一定关联性。

图5分析了侧压系数, 即(σH+σh)/2与σv的比值k,随埋深变化的规律。能够看出, 随着埋深的不断增加, 侧压系数k的分布范围逐渐收窄。当埋深大于300 m时,k值收敛于1附近。k值随埋深的变化规律为k=88/H+0.72, 且35/H+0.29

图5 布尔台矿区侧压系数随埋深变化规律Fig.5 Lateral pressure coefficient variation law of the Buertai mine area with the change of depth

由于钻孔位置集中在布尔台矿三、四盘区, 为验证地应力测量结果在全矿井的适用性, 选取布尔台矿一、二盘区3个钻孔的空心包体法地应力测试结果[27]与三、四盘区埋深相近的Kaiser法效应测点测量得到的3组地应力值进行对比。对比结果见表2, 表2中BK209–21, BK209–22与BK209–23为位于三盘区的Kaiser法地应力测点。BK1号, BK4号与BK5号为位于一、二盘区的空心包体法地应力测点。

表2 一、二盘区与三、四盘区地应力测试结果对比Table 2 Comparison of in-situ stress test results between the first and second plates, and those from the third and fourth plates

由表2可知, 一、二盘区σH测量值与三、四盘区σH测量值间相对误差平均值为14.90%,σh相对误差平均值为19.67%,σv相对误差平均值为15.47%。因此, 一、二盘区地应力测试结果与三、四盘区地应力测试结果误差较小, 同时考虑到测试方法的不同, 误差在工程合理范围内, 能够证明Kaiser法地应力测量结果在全矿井的适用性。

2 基于改进组合弹簧模型的地应力计算方法

尽管经典组合弹簧模型考虑了岩层物理力学参数对地应力的影响, 能较为准确反映地应力场的分布规律。但该方法假设地层在水平方向最大和最小主应变为恒定值。因矿区地域范围广, 地质构造复杂, 地应力场受岩体成因、岩体结构和地质构造运动等多因素的影响, 表现出显著的非均质性。故组合弹簧模型无法准确描述区域地应力场分布。

为提高地应力场反演计算精度, 将克里金插值方法引入到经典组合弹簧模型的应变计算过程中,提出一种基于Kaiser法实测地应力数据与改进组合弹簧模型的地应力场计算方法。该方法在考虑岩层物性的同时克服了组合弹簧模型理论中应变参量无法反映围岩变形非均质性的问题。

2.1 地应力测点处应变计算

基于Kaiser法测定的地应力数据, 使用组合弹簧模型, 结合实测地应力值与岩石物理力学参数求取各地应力测点处的应变。最大水平主应变与最小水平主应变的计算公式为

式中,Hε为该点最大水平主应变;hε为该点最小水平主应变;ν为该点岩石泊松比;Hσ为该点最大水平主应力, GPa;hσ为该点最小水平主应力, GPa;vσ为该点竖直方向主应力, GPa;E为该点岩石杨氏模量, GPa。

2.2 研究区域任意点应变插值计算

按测定的最大水平主应力方向将应变分解为x方向上的正应变εx,y方向上的正应变εy及剪应变γxy。基于地应力测点处的应变, 使用克里金插值方法求取应变场空间分布。克里金插值法可考虑空间数据的各向异性、非线性和多变量等特征, 其计算式如下:

使用克里金插值法得到研究区域内任一点x方向上的正应变εx,y方向上的正应变εy及剪应变γxy。

2.3 研究区域内任意点地应力计算

根据研究区域内任一点的应变张量, 进而可计算出该点最大水平主应变与最小水平主应变及方向, 计算公式为

式中,θ为最大主应变方向;xε为x方向上的正应变;yε为y方向上的正应变;γxy为xy方向上的剪应变。

在得到任一点的水平方向最大和最小主应变大小和方向后, 可结合该点的杨氏模量E(GPa)、泊松比ν等岩石物理力学性质, 根据组合弹簧模型计算出该点的地应力, 具体算公式为

3 结果及分析

3.1 地应力反演结果分析

为验证改进的组合弹簧模型地应力场反演方法的有效性, 笔者将其与侧压系数法[28]、线性拟合法[29]、分段拟合法[30]及经典组合弹簧模型方法[31]进行了对比, 图6为地应力反演方法对比评估流程。

图6 地应力反演方法对比评估流程Fig.6 Comparative evaluation process of in-situ stress inversion methods

基于经古地磁岩芯定向的15个地应力测点数据, 使用各方法进行地应力反演, 随后测试各方法对71个非定向地应力测点的反演计算精度, 从而对比各拟合方法的计算准确性。

3.1.1 侧压系数拟合结果

侧压系数法是一种简单的地应力反演方法, 该方法假设岩石在加载过程中保持弹性, 水平方向地应力大小与竖直方向地应力大小成正比, 比例系数为侧压系数λ。侧压系数计算公式为

式中,n为地应力测点数量;σHi为第i个测点的最大水平主应力, GPa;σhi为第i个测点的最小水平主应力, GPa;σvi为第i个测点的竖直方向主应力, GPa。

基于实测地应力计算得到布尔台矿侧压系数为1.034。侧压系数拟合得到研究区域任意点的水平方向地应力为

式中,H为埋深, m;ρ为覆岩密度, kg/m3;g为重力加速度, 取9.8 m/s2。

3.1.2 线性函数拟合结果

线性拟合法是一种常用的地应力反演方法, 它利用最小二乘法对实测数据进行线性拟合, 得到地应力的空间分布, 其拟合结果如图7所示。

图7 布尔台矿地应力线性拟合结果Fig.7 Linear fitting results of in-situ stress in Buertai Mine

使用线性拟合方法得到的地应力拟合公式为

3.1.3 分段函数拟合结果

分段函数拟合法对地应力数据进行分段回归,每100 m埋深设定为一个分段, 回归表达式为

式中,β0,…,β6为函数系数;xi1,…,xi5为虚拟变量。

当H<100 m时,xi1,…,xi5均为0; 当100

图8 布尔台矿地应力分段拟合结果Fig.8 Piecewise fitting results of the in-situ stress in Buertai Mine

使用分段拟合方法得到的拟合式为

3.1.4 经典组合弹簧模型反演结果

经典组合弹簧模型地应力反演方法通过求取最佳水平方向地应力系数确定地应力分布, 基于地应力实测数据, 使用经典组合弹簧模型求得的平均最大主应变与平均最小主应变分别为1.04×10–3与4.56×10–4。使用经典组合弹簧模型方法得到的地应力拟合式为

3.1.5 改进组合弹簧模型反演结果

使用笔者所提出的改进组合弹簧模型得出的研究区域地应力各分量随埋深变化曲线如图9所示。

图9 布尔台矿地应力改进组合弹簧模型反演结果Fig.9 Inversion results of the improved composite spring model for the in-situ stress in Buertai Mine

3.2 地应力反演效果比较

使用均方根误差(RMSE)与平均绝对百分比误差(SMAPE)指标对各类地应力反演方法得到的σH与σh的绝对误差与相对误差进行评估。RMSE与SMAPE的计算方法为

式中,n为样本个数;yi为第i个样本的实际值;为第i个样本的预测值。

表3给出了各反演计算模型对比分析部分结果。

表3 各反演方法计算误差Table 3 Errors of different inversion methods.

由表3可知, 改进组合弹簧模型反演得到的Hσ与实测值间的均方根误差为1.849 8, 较侧压系数法的3.153 2降低41.3%; 较线性拟合法的3.272 8降低43.5%; 较分段拟合法的3.066 3降低39.7%, 较经典组合弹簧模型反演的2.305 1降低19.8%。实测数据反演法得到的hσ与实测值间的均方根误差为1.811 9, 较侧压系数法的3.183 8降低43.1%; 较线性拟合法的2.607 4降低30.5%; 较分段拟合法的2.448 5降低26.0%, 较经典组合弹簧模型反演的1.981 5降低8.6%。

此外在评价反演结果与实测数据间的相对误差时, 改进组合弹簧模型反演得到的Hσ与实测值间的SMAPE值为17.27%, 较侧压系数法的26.86%降低35.7%, 较线性拟合法的23.52%降低26.6%, 较分段拟合法的25.12%降低31.3%, 较经典组合弹簧模型的18.14%降低4.8%。改进组合弹簧模型反演得到的hσ与实测值间的SMAPE值为21.62%, 较侧压系数法的34.77%降低37.8%, 较线性拟合法的31.14%降低30.6%, 较分段拟合法的28.76%降低24.8%, 较经典组合弹簧模型的22.14%降低2.3%。

综上, 改进组合弹簧模型计算的最大水平主应力及最小水平主应力结果与实测值间的均方根误差为1.830 9, 相较其他方法减少14.6%~42.2%, 平均绝对百分比误差为19.45%，相较其他减少3.4%~36.9%。所提出的改进组合弹簧模型方法计算误差最小, RMSE和SMAPE值均最小, 说明其拟合效果和稳定性最好, 能有效反映水平地应力与应变的非线性关系。除改进组合弹簧模型外, 经典组合弹簧模型的RMSE与SMAPE最小。相比之下, 侧压系数法和线性拟合法的RMSE与SMAPE较大, 说明其拟合效果与稳定性较差。

3.3 数据敏感性分析

为分析数据量对改进组合弹簧模型拟合精度的影响, 进行了改进组合弹簧模型的数据敏感性分析。由于克里金插值需要最少3个数据点, 故分析了数据点数量从3增加到15的过程中拟合精度的变化。结果显示, 随着数据点的增多, 拟合精度逐渐上升。其中, 最大水平主应力的RMSE值从4.88下降至2.56, SMAPE值从26.91%下降至25.83%; 最小水平主应力的RMSE值从6.13下降至1.51, SMAPE值从49.75%下降至18.93%。同时, 增加同样数量数据点对拟合精度的提升逐渐下降。在数据量增加至12~15时, 各拟合精度指标均不再发生显著下降。因此继续增加数据点对精度提升作用有限, 故15个数据点能够较好地拟合该区域地应力。各精度指标随数据量变化的曲线如图10所示。

图10 改进组合弹簧模型拟合误差随数据量变化曲线Fig.10 Error curves of the improved combined spring model varies with the amount of fitted data

3.4 模型应用效果评价

为进一步验证所提出方法的有效性和优越性,将地应力反演计算方法融入FLAC3D软件建模。分别以改进组合弹簧模型方法与侧压系数方法计算得到的结果作为数值模型应力边界条件进行模拟,并与实测矿压数据进行对比。从而评估不同地应力计算方法对模拟准确性的影响。

使用FLAC3D软件模拟布尔台矿42105工作面的开采过程, 在45号及85号支架位置处设置支架阻力测线。

布尔台矿42105工作面布置如图11所示。

图11 布尔台矿42105工作面平面布置Fig.11 Layout of No.42105 working face in Buertai Mine

3.4.1 模型建立及岩石力学参数赋值

布尔台矿42105工作面位于4–2煤一盘区, 走向长5 231 m, 倾向长230 m。基于SJ–3–1钻孔与SJ–3–2钻孔信息构建布尔台42105工作面数值模型, 钻孔位置如图11所示。数值模型长900 m, 宽400 m, 高487.62 m, 共包含2 380 800个单元, 2 440 125个网格点。根据布尔台矿42105工作面钻孔岩芯的物理力学性质测试数据及以往文献[32]中的神东矿区岩石物理力学性质, 确定了数值模型的岩石物理力学参数, 具体见表4。

表4 布尔台矿42105工作面岩石物理力学参数Table 4 Physical and mechanical parameters of No.42105 working face in Buertai Mine

3.4.2 地应力分块精细化加载

为将地应力精确加载到数值模型上, 将工作面的侧面应力边界划分为4 312个矩形区域, 划分情况如图12所示。

图12 数值模型应力精细加载方法示意Fig.12 Diagram of numerical model based on the stress fine loading method

在FLAC3D中, 应力边界条件可表示为

式中,Cσ为坐标点(0,0)处的应力值; ∂σ/∂x为应力沿x方向的梯度; ∂σ/∂y为应力沿y方向的梯度。

任意矩形区域的应力边界如图13所示。

图13 FLAC3D应力边界条件示意Fig.13 Diagram of the stress boundary condition in the model of FLAC3D

已知矩形区域4个边界点处的应力, 求取该区域的应力边界条件等价于求解以下矩阵:

式中, (xi,yi)为边界点Pi的坐标;σi为边界点Pi处的应力。

对于任意一个矩形区域, 其应力边界条件可以由3条FLAC3D指令确定, 分别用于指定坐标原点处σx,σy与τxy的值及其沿x,y,z方向的梯度。使用Python编写边界条件求解程序, 运行后共生成12 936条边界条件加载指令。最后使用FLAC3D软件逐条运行边界加载指令, 以此实现应力边界条件的精细化加载。

3.4.3 工作面支架支撑力分布

采用FLAC3D模拟布尔台矿42105工作面开采过程, 从开切眼位置开始每次向工作面推进方向开挖3.5 m, 共开采70 m。对邻近工作面的一列顶板单元施加支撑力从而体现液压支架对顶板的支撑作用。FLAC3D模型每运算50步后, 根据顶板下沉量更新支撑力, 直到模型最大不平衡力比率小于10–5。记录每一步开挖运算至平衡时支撑力的值。为便于与现场实测数据比较, 根据布尔台矿42105工作面使用的ZFY18000/25/39D型液压支架参数[33]将支撑力换算为液压支架立柱压强。

随后将基于改进组合弹簧模型法、经典组合弹簧模型法与侧压系数法计算得到的3种应力边界条件下的数值模拟结果与现场实测数据进行对比,结果如图14所示。

图14 支架支柱压强实测值与不同应力边界条件模拟结果Fig.14 Measured values of support pressure and the simulated results of different stress boundary conditions

由图14可知, 使用侧压系数法得到的应力场作为数值模型的水平应力边界条件时, 最大支架阻力在开采至35 m处出现, 来压时45号支架立柱压强为40.8 MPa, 85号支架立柱压强为40.6 MPa。使用经典组合弹簧模型方法得到的应力场作为数值模型的水平应力边界条件时, 最大支架阻力在开采至35 m处出现, 来压时45号支架立柱压强为64.7 MPa,85号支架立柱压强为63.0 MPa。使用改进组合弹簧模型法得到的应力场作为数值模型的水平应力边界条件时, 最大支架阻力出现在开采至56 m处时,来压时45号液压支架立柱压强为47.8 MPa, 85号支架立柱压强为46.6 MPa。根据布尔台矿42105工作面液压支架实测数据[34], 首次基本顶破断引起的矿压显现发生于工作面开采到52.0~61.5 m时, 来压期间45号支架峰值立柱压强为49.9 MPa, 85号支架峰值立柱压强约为48.0 MPa。因此, 相较于使用侧压系数计算结果作为数值模型的水平应力边界条件,利用改进组合弹簧模型的计算结果能更准确地模拟现场实际情况。

采用均方根误差(RMSE)与平均绝对百分比误差(SMAPE)指标对使用3种应力场计算结果作为应力边界条件的模拟结果与实测矿压数据间的差异进行定量评估, 指标定义见式(14), 评估结果见表5。

表5 各计算方法数值模拟结果与实测数据间误差Table 5 Error between numerical simulation results and measured data

由表5可以发现, 使用改进组合弹簧模型法得到的地应力场作为模型水平方向应力边界时, 矿压模拟值与实测值间的误差较小。

对于45号支架的监测数据而言, 改进组合弹簧模型法与实测值间RMSE为68.10, 较侧压系数法的75.53低9.84%, 较经典组合弹簧模型的103.16低33.99%。SMAPE为16.95%, 较侧压系数法的18.37%低7.73%, 较经典组合弹簧模型的20.90%低18.90%。此外, 对于85号支架的监测数据而言, 改进组合弹簧模型法与实测值间RMSE为65.00, 较侧压系数法的68.76低5.47%, 较经典组合弹簧模型的93.89低30.77%。SMAPE为14.54%, 较侧压系数法的15.35%低5.28%, 较经典组合弹簧模型的18.81%低22.70%。

综上, 模拟得到的工作面矿压数据与实测数据间均方根误差平均值为66.55, 较使用侧压系数法降低7.75%, 较使用经典组合弹簧模型法降低32.45%。平均绝对百分比误差平均值为15.75%, 较使用侧压系数法降低6.61%, 较使用经典组合弹簧模型法降低20.67%。通过对比以改进组合弹簧模型与侧压系数法的计算结果作为数值模拟边界条件的模拟结果与实测值间的误差, 证明了所提出的方法较其他模型能更准确地模拟现场实际情况。从而为使用数值模拟方法进行工作面矿压显现分析, 强矿压防治以及巷道优化布置等提供一种新的数值模型应力边界条件确定方法, 具有重要的理论指导意义和工程应用价值。

4 结 论

(1)以布尔台矿为研究对象, 基于Kaiser法测量地应力数据。测量结果显示: 同一钻孔中的地应力测点, 随深度的增加, 最大水平主应力σH、最小水平主应力σh和竖直主应力σv均逐渐增加。随着埋深不断增加, 侧压系数k的分布范围逐渐变小。当埋深大于300 m时,k值收敛于1附近。k值随埋深的变化规律为k=88/H+0.72, 且35/H+0.29

(2)将克里金插值方法引入经典组合弹簧模型,提出了一种同时考虑岩层物性与围岩变形非均质性的改进组合弹簧模型。通过比较改进组合弹簧模型、侧压系数法、线性拟合法、分段函数拟合法和经典组合弹簧模型与实测数据间的误差, 评估了各方法的计算精度。结果表明改进组合弹簧模型计算的最大水平主应力及最小水平主应力结果与实测值间的均方根误差为1.8309, 相较其他方法减少14.6%~42.2%, 平均绝对百分比误差为19.45%,相较其他减少3.4%~36.9%。

(3)以布尔台矿42015工作面为工程背景, 使用改进组合弹簧模型方法计算得到的结果作为FLAC3D数值模型应力边界条件进行模拟。模拟得到的工作面矿压数据与实测数据间均方根误差平均值为66.55, 较使用侧压系数法降低7.75%, 较使用经典组合弹簧模型法降低32.45%。平均绝对百分比误差平均值为15.75%, 较使用侧压系数法降低6.61%, 较使用经典组合弹簧模型法降低20.67%。结果表明采用改进组合弹簧模型的计算结果作为数值模型边界条件较其他模型能更准确地模拟现场实际情况。