一种厚硬顶板破断诱发冲击地压的临界判据方法

本发明属于矿山安全与岩石力学工程，具体涉及一种厚硬顶板破断诱发冲击地压的临界判据方法。

背景技术：

1、冲击地压是深部煤矿开采中极具破坏力的动力灾害，而厚硬顶板的突然破断是诱发该灾害的主要诱因之一。准确判别厚硬顶板何时发生破断以及其破断是否足以引发冲击地压，是采场围岩控制与灾害防治的核心难题。

2、目前，国内外学者针对厚硬顶板破断的判据研究，主要基于以下两类技术路线：

3、（1）基于应力或变形的静态判据：通过监测顶板内部应力（如钻孔应力计）或岩层变形（如离层仪、多点位移计），设定经验阈值。当监测值超过阈值时，判定为危险。例如，文献《坚硬顶板破断的矿压显现与冲击地压关系研究》（采矿与安全工程学报，2018, 35(2):256-262）中，采用顶板极限悬露步距作为破断判据。然而，此类方法忽略了破断瞬间能量释放的动态过程和释放速率，导致预警滞后，且阈值普适性差。

4、（2）基于数值模拟的能量积聚判据：采用有限元或离散元软件（如flac3d、udec）计算顶板内积聚的弹性能总量，当总能量超过某一预设值时即判为危险。例如，专利cn104929664a“一种预测采场冲击地压危险性的方法”通过模拟计算顶板挠曲变形能来评估风险。但此类方法计算复杂，难以实现实时预警，且未考虑能量释放率（单位时间或单位破断面积所释放的能量）这一关键动态指标，而能量释放率的高低才是决定冲击荷载大小的直接因素。

5、综上所述，现有技术存在的根本缺陷在于：缺乏一个能够量化厚硬顶板破断瞬间能量释放率“突变”特征（即从缓慢积聚到瞬间释放的剧烈转变），并以此建立可解析计算的临界判据的方法。这导致对“破断是否诱发冲击”的判别要么滞后，要么不准确，无法满足现场即时、可靠的预警需求。

技术实现思路

1、针对上述现有技术中存在的不足，本发明旨在提供一种基于能量释放率突变特征的厚硬顶板破断诱发冲击地压的临界判据方法。该方法能够定量描述顶板破断前后能量释放率从“缓变”到“突变”的动态过程，并建立明确的数学临界条件，从而实现对冲击地压发生时刻的准确、超前判别，提高矿井动力灾害预测的可靠性与实用性。

2、为了解决上述技术问题，本发明采用了如下的技术方案：

3、一种厚硬顶板破断诱发冲击地压的临界判据方法，该临界判据方法包括如下步骤：

4、s1：采集厚硬顶板关键物理力学参数并建立原位应力环境数据库

5、冲击地压发生的本质是顶板系统从稳态能量积聚到失稳释放的突变过程，该过程高度依赖顶板岩体的本征物理属性及其所处的地应力场。因此，第一步必须通过现场探测与室内试验，获取一组能够表征岩体储能极限、断裂韧性与动态响应的基础参数；与常规方法不同，该步骤额外要求测定损伤局部化系数和非局部影响半径，这两个参数将用于后续步骤s2中的非局部损伤梯度模型。

6、定义物理量及其符号（均采用国际单位制）：

7、

8、s2：建立基于非局部损伤梯度的顶板能量存储模型

9、经典线弹性断裂力学假设裂纹尖端为奇点，无法描述真实岩石中存在的损伤过程区。本发明引入非局部损伤梯度理论，将厚硬顶板视为一个含弥散裂纹的连续介质，其损伤变量d(x)不仅取决于该点的应变，还受邻域损伤状态的影响；顶板内存储的总能量由弹性应变能和损伤表面能两部分组成；

10、非局部损伤梯度下的总势能泛函

11、

12、式中：ω为顶板占据的二维或三维连续域；σij为应力张量，εij为应变张量，满足本构关系σij=cijklεkl，其中cijkl为弹性刚度张量，εkl为二阶应变张量的分量形式；d(x)为损伤变量，取值范围0≤d≤1，d=0表示完整岩石，d=1表示完全破坏；β为梯度损伤模量（n），为材料常数，通过拟合单轴拉伸下的应变软化曲线确定；为非局部影响半径（即步骤s1中的参数），控制损伤非局部相互作用的距离；为损伤场梯度的平方，表征损伤的空间不均匀程度；为局部损伤耗能函数，通常取（线性软化）或（二次软化），其中为i型断裂能；本发明采用双线性形式以更好地匹配硬岩脆性行为：

13、

14、该式确保当d=1时，单位体积耗散的断裂能为，积分后总断裂能与一致；

15、γt为边界上施加外力的部分；ti为面力分量；ui为位移分量；

16、该势能泛函首次将梯度损伤理论与顶板破断诱发冲击地压问题结合，其中的项能够模拟裂纹尖端损伤带的宽度效应，避免了经典线弹性断裂力学中的应力奇异性，为后续计算能量释放率的连续演化提供了物理基础；

17、步骤s2的输出量：顶板在任意载荷状态下的总势能π及其对裂纹长度的变分导数，该导数将在步骤s5中与能量释放率建立联系。

18、s3：推导弯曲变形下顶板内部的弹性能密度分布解析式

19、厚硬顶板在上覆岩层载荷及水平构造应力作用下发生弯曲变形，为准确计算存储的弹性能，需要获得板内任一点的应力状态；本发明采用弹性薄板理论，并考虑双向应力状态，推导出能量密度随板厚h和悬露长度l的分布函数；不同于kirchhoff经典假设，本发明额外引入了剪切变形修正系数以适用于厚跨比较大的坚硬顶板（h/l>0.1）;

20、考虑剪切变形的能量密度函数：

21、设顶板为两端固支的矩形板（宽度方向取单位长度），上表面受均布载荷q（mpa），同时存在沿悬露方向的水平构造应力σh和沿宽度方向的水平构造应力σh2，通常取σh2=λhσh，λh为侧压系数；取坐标：x沿悬露方向，原点在固支端；z垂直向下，z=0为中面；

22、应力分量：

23、

24、弯曲正应力：

25、

26、其中，为弹性模量，为中面法向坐标，为泊松比（无量纲），w(x)为挠度函数，对于均布载荷下的固支梁，近似取：

27、

28、式中：为面荷载，为顶板的悬露长度，为顶板的厚度，为单位宽度惯性矩，为剪切修正项为剪切修正项；

29、总正应力：

30、

31、

32、剪应力（由横向荷载引起）：

33、

34、式中：为剪力分布；

35、弹性能密度（单位体积储存的能量）：

36、

37、式中：为剪切模量，为x 方向（悬露方向）正应力，为y 方向（宽度方向）正应力，为x-z 平面内切应力的平方。

38、s4：定义等效裂纹扩展驱动力并引入非局部断裂准则

39、顶板中的原生裂隙或采动裂隙在应力作用下扩展；经典断裂力学采用应力强度因子ki与断裂韧度kic的比较作为起裂准则，但该准则无法描述裂纹扩展过程中的能量耗散历史；本发明提出等效裂纹扩展驱动力jeq，它基于j-积分但进行了非局部平均，以消除尖端奇异场的影响；

40、非局部等效驱动力

41、

42、式中： 为当前主裂纹半长或全长（根据几何模型），单位为m；j()为经典的ricej-积分，对于平面应变条件下受远场应力的裂纹，有：

43、

44、式中：，e′为平面应变等效弹性模量，为i 型（张开型）应力强度因子，为主裂纹长度，ν为泊松比，jsh为剪切分量贡献；

45、为步骤s1的非局部影响半径，作为指数权重函数的衰减长度；为指数权重函数，赋予靠近当前裂纹尖端（r小）更高的权重，而远处的影响呈指数衰减；r为积分变量，表示距离当前裂纹尖端的偏移距离；

46、物理意义：

47、jeq可以理解为裂纹尖端前方一个损伤过程区内平均的能量释放率，当jeq达到材料的非局部断裂韧度 jicnl时，裂纹开始稳态扩展；其中：

48、

49、式中，为初始裂纹尺寸，该修正项反映了小尺度屈服或损伤区对断裂韧性的提高效应；

50、原创性说明：传统j-积分仅在裂纹尖端满足路径无关，但真实岩石中由于损伤区的存在，j值依赖于路径。本发明引入指数加权非局部平均，首次将（非局部半径）嵌入驱动力定义中，使得判据具有内在长度尺度，解决了网格依赖性难题，并自然地引出后续的能量释放率突变分析。

51、s5：构建能量释放率geff()随裂纹长度的动态演化方程

52、能量释放率g定义为裂纹扩展单位面积所释放的系统势能，在非局部损伤模型中，g不是瞬时局部量，而是一个依赖于裂纹扩展历史的累积效应；本发明推导出有效能量释放率geff()的积分微分方程，它同时包含弹性释放部分和损伤耗散部分；

53、有效能量释放率的演化方程

54、

55、式中：为裂纹扩展中间长度，为权重函数，g0()为弹性名义能量释放率，假设材料无损伤时由弹性给出；对于受弯曲和压缩的顶板，推导得：

56、

57、式中:为水平构造应力，为裂纹闭合影响系数， y(l/)为几何修正因子，对于有限宽板（经验拟合多项式）；为累积损伤耗散，表示裂纹从初始长度0扩展到′过程中消耗于损伤过程区的能量，其表达式为：

58、

59、其中，w(d)已在步骤s2中定义，为裂纹尖端位置，为裂纹尖端的厚度方向坐标；为记忆核函数，为记忆核函数，表征历史裂纹扩展对当前能量释放率的影响；本发明取为：

60、

61、其中，为记忆长度（取为0.5m，与顶板厚度同量级），ξ为相对距离变量，；0为初始裂纹长度（步骤s1中实测或假设为0.05m）；

62、该方程表明，真实的能量释放率geff小于弹性解g0，因为一部分能量被损伤区耗散；但更关键的是，当裂纹扩展加速时，历史耗散项对当前geff的“阻碍”作用减弱，导致geff在临界点附近出现跳跃；这正是能量释放率突变特征的理论源头；

63、创新点：本步骤首次将记忆核函数引入能量释放率演化，借鉴了粘弹性力学中的本构方程形式，并转移应用于断裂过程。这使得模型能够描述裂纹扩展速率的历史依赖性，为后续突变因子的计算提供了时间（或路径）累积效应。

64、s6：计算能量释放率对裂纹长度的二阶导数，揭示突变前兆信号

65、在经典突变理论中，系统失稳的临界点对应势能函数的三阶导数为零或状态变量的一阶导数趋于无穷。本发明不直接使用势能，而是分析能量释放率geff()作为裂纹长度的函数；当系统接近失稳时，dgeff/d急剧增大，而d2geff/d2会从负值变为正值，产生一个特征性的“拐点”，该拐点即为突变的前兆信号；

66、理论公式：能量释放率的一阶与二阶灵敏度

67、首先，对步骤s5中的演化方程关于求导：

68、

69、式中，dg0/d为弹性名义能量释放率的增长率，由步骤s5中g0显式求导；利用莱布尼茨法则，得到：

70、

71、其中：为裂纹尖端自身位置的权重系数，

72、再求二阶导数：

73、

74、

75、式中 为记忆核函数的二阶导数，,其中；由于(ξ)是偶函数且在ξ=0处取最大值，∂(0)/∂a=0；因此上式简化为：

76、

77、式中，为损伤耗散能对裂纹长度的二阶变化率，通过数值差分获得;

78、突变前兆判据：

79、当 =0时，对应geff()曲线的拐点；在该点之前，d2geff/d2<0，曲线上凸，能量释放率增长趋缓；过该点后，d2geff/d2>0，曲线下凹，能量释放率加速上升，系统进入非稳定扩展阶段，该拐点对应的裂纹长度记为。

80、符号补充说明：

81、

82、s7：定义能量释放率突变因子λ的新表达式（包含高阶梯度项）

83、现有技术中的突变因子通常定义为峰值与平均值的比值（如背景技术中提到的λ=gpeak/gpre），该定义忽略了能量释放率变化速率的突变特征。本发明提出一个无量纲突变因子λ，它同时包含能量释放率的幅值跃升和其变化率的阶跃，且引入非局部梯度项以反映空间突变的剧烈程度;

84、理论公式：多尺度突变因子λ

85、

86、式中，gpeak为裂纹失稳扩展（破断）瞬时的峰值有效能量释放率，它不等于静态断裂能gf，而是动态效应导致的放大值，将在步骤s8中具体计算；gpre为失稳前稳态扩展阶段的平均能量释放率，取裂纹长度从0.9c到0.99c区间内geff的算术平均值，其中c为临界裂纹长度（满足jeq(c)=jicnl）；β1为幅值跃升敏感指数，由岩性决定，取值范围，初始取β1=1.0，后续通过s12的自适应学习调整；dgeff/d为能量释放率对裂纹长度的导数，c−表示从小于c的方向趋近，c+表示从大于c的方向趋近（由于破断发生后裂纹长度瞬时跳跃，该右导数在经典模型中无定义，但本发明采用极限外推）；β2为速率比敏感指数，取值范围，通常取β2=0.5；为能量释放率场的拉普拉斯算子，在裂纹尖端附近计算；对于二维问题，，其值反映了geff空间分布的不均匀程度；为能量释放率梯度的模，;ℓ为非局部影响半径（同步骤s1）;第三项 称为空间突变增强因子；当裂纹尖端前方出现高度集中的能量释放率梯度时，为较大的正值，该因子大于1，显著提升λ;

87、物理意义：

88、λ综合了三个维度的突变特征：

89、（1）能量幅值的跃升；

90、（2）能量释放率变化率的逆转（因为在失稳前 dg/d>0 且递增，失稳瞬间由于裂纹跳跃，dg/d 理论上趋于无穷，所以比值极大）；

91、（3）空间梯度场的奇异性增强。

92、s8：建立考虑动态破裂速度的峰值能量释放率gpeak计算模型

93、顶板破断瞬间，裂纹以高速扩展（可达材料瑞利波速的0.3-0.9倍），此时惯性效应导致动态应力强度因子显著高于静态值，从而释放的能量大于静态断裂能；本发明基于动态断裂力学中freund公式，提出一个适用于岩体的动态放大系数，并引入裂纹扩展速度v的估算方法；

94、动态峰值能量释放率

95、

96、式中：gstatic(c)为在裂纹长度达到临界值c时的静态有效能量释放率，即geff(c)；f(vcr)为动态放大系数，无量纲；根据freund对于半无限裂纹的解析解并修正适用于有限几何体，本发明提出以下经验-理论混合公式：

97、

98、式中，v为裂纹扩展速度（m/s）；cr为顶板岩层的瑞利波速（m/s），，其中μ为密度（kg/m3）；γ为材料常数，反映岩体脆性程度；，量纲为1，对于硬岩系砂岩，γ通常在0.5~1.2之间；

99、裂纹扩展速度vcr的估算：

100、在临界失稳点，裂纹扩展速度由能量释放率与断裂能的平衡决定；本发明采用如下能量平衡方程求解v：

101、

102、式中，为纵波波速，右边为动态断裂能，考虑了效应；该方程为超越方程，可通过数值迭代法（如牛顿法）求解v，代入f(v)即得gpeak；

103、创新点：传统方法在定量预测完全忽略裂纹扩展速度的影响，本发明首次将动态断裂力学的放大效应定量引入峰值能量释放率的计算中，使得λ中具有真实的物理意义。

104、s9：推导基于最小耗能原理的临界突变阈值λc

105、临界阈值λc是判别冲击地压是否发生的分界线，现有技术多采用经验值或统计回归，缺乏理论依据。本发明基于非平衡热力学中的最小耗能原理，推导出λc关于岩体基本参数的显式表达式；该原理指出：系统在失稳时总是选择能量释放率最大的路径，而临界状态对应于耗能对控制参数的二阶变分为零；

106、临界突变因子的理论表达式

107、

108、式中：gf为材料的断裂能（j/m2），由步骤s1获得；g0(0)为初始裂纹长度下的弹性名义能量释放率，按步骤s5中公式计算；比值gf/g0(0) 表征了材料抵抗裂纹起裂的能力相对于当前驱动力的富裕程度，该比值越大，需要越高的突变因子才能触发冲击，因此λc与之正相关；为非局部影响半径（m）；为裂纹从初始长度扩展到临界长度的增量；比值反映了损伤过程区尺度相对于裂纹扩展距离的重要性，当 较大时，表示损伤区域宽，系统需要更大的能量释放率突变才能克服耗散，因此λc增大；exp(σt/σc)为指数项，来源于最小耗能原理推导过程中的玻尔兹曼因子；物理意义：抗拉强度σt相对于抗压强度σc越大，岩石越显脆性，越容易发生突变，因此指数项大于1，且随σt/σc增大而迅速增大整个表达式量纲为 1，且所有参数均可通过步骤s1试验获得；

109、推导要点：从系统总势能π出发，考虑裂纹扩展的不可逆性，建立耗散；最小耗能原理要求在临界点成立；经过变分运算，并假设损伤演化满足形式，最终得到上述λc表达式；λc是一个仅依赖材料本征属性的常数，而非经验拟合值。

110、s10：提出考虑围岩刚度和支护反力的修正临界条件

111、厚硬顶板并非孤立结构，其上下方有煤层、直接顶及底板，且采场存在支架支护；围岩刚度（即顶板弯曲时周边岩体对其的约束能力）和支护反力会吸收一部分能量，从而改变破断瞬间的能量释放率突变特征；对步骤s9中的λc进行修正，引入围岩刚度系数ksurr和支护等效阻尼csupp；

112、修正后的临界突变因子（初始修正版）

113、

114、式中：kplate为顶板本身的弯曲刚度（n/m），对于单位宽度悬臂或固支板，（固支梁中点刚度）；ksurr为围岩等效刚度（n/m），通过现场原位压板试验或数值反演获得；硬岩顶板；第一修正因子；

115、物理意义：围岩刚度越大，顶板系统整体越“硬”，能量释放越集中，因此临界突变因子应减小（即更容易达到临界）；但注意公式中乘以小于1的因子，使得λcmod<λc，表示更易触发冲击； csupp为支护系统的等效阻尼系数，由支架的让压特性、立柱安全阀流量决定，可以通过液压支架的力-缩量滞回曲线积分得到；v为裂纹临界扩展速度（步骤s8中计算）；第二修正因子：指数项小于1，表示支护系统通过阻尼耗散能量，提高了临界阈值，因为λcmod更小，实际需要的λ更易达到，原判据是λ≥λc触发冲击；修正后λcmod变小，意味着更小的λ即可触发，这似乎与支护提高安全性的直觉相反；应修正为：支护吸收能量，需要更大的突变才能诱发冲击，因此临界因子应增大；为此，本发明重新定义

116、修正后正确的判据应为：

117、

118、这样，围岩刚度越大（ksurr大），λcmod增大；支护阻尼越大，λcmod也增大，意味着需要更强的能量释放率突变才能诱发冲击，符合物理实际；

119、最终临界条件：

120、；

121、创新点：本步骤首次将围岩刚度和支护阻尼以指数及线性形式嵌入临界判据，使得该方法能够适应不同采区地质条件和支护参数，具有高度的工程适应性。

122、s11：现场实时计算流程与动态预警规则

123、将上述理论转化为可操作的现场预警系统，需要明确数据采集、参数更新、λ计算、与修正临界值比较的时序流程；本步骤设计一套五级预警机制，并提供简化计算流程以降低实时计算负担。

124、s11.1 数据采集与预处理（每30分钟一次）

125、通过微震监测系统定位顶板内的微震事件，采用聚类算法识别裂纹群，估算当前等效裂纹长度；

126、通过矿压监测系统获取顶板悬露长度l；

127、通过钻孔应力计获取水平构造应力σh；

128、s11.2 参数更新

129、根据当前l，，调用步骤s5中的g0()公式；

130、调用步骤s6中的dgeff/d和d2geff/d2（采用有限差分近似，无需实时求解积分方程，可预先计算不同/l下的导数表）。

131、调用步骤s7中的λ公式，其中gpre取过去24小时内geff的滑动平均值（窗口长度5个采样点）；

132、s11.3 计算修正临界值

133、根据步骤s1测定的本征参数，计算λc（一次性离线计算）；

134、根据当前采场支护状态（阻尼csupp由支架压力监测动态评估），计算λcmod；

135、s11.4 预警等级判定

136、设比值；

137、

138、s11.5 动态反馈修正

139、每次发生顶板破断（无论是否诱发冲击），记录实际破断时的λ实际，通过卡尔曼滤波更新β1,β2等敏感指数（见步骤s12），实现自适应学习。

140、s12：验证与自适应参数调整方法

141、理论模型中的β1,β2,γ,等经验或半经验参数，在初始阶段采用经验值；但随着现场推进和多次顶板破断事件的记录，可以通过贝叶斯反演或递推最小二乘法进行自适应修正，以提高预报准确率；本步骤给出具体的参数在线更新算法；

142、s12.1 目标函数定义

143、对于第k次顶板破断事件（k=1,2,...），设模型预测的临界突变因子为λc(k)（理论值），实际触发冲击时的λ实际(k)（由步骤s11中实测数据反演得到）；定义误差：

144、

145、s12.2 待修正参数向量

146、

147、各参数初值：β1(0)=1.0；β2(0)=0.5；γ(0)=0.8（典型砂岩）；；

148、s12.3 递推最小二乘更新公式

149、每次获得一新的事件数据，更新参数：

150、

151、式中：为第k步更新后的待估参数向量，为第(k-1)步的参数估计值，为第k步的观测残差，，为灵敏度向量，可通过有限差分计算，其中，为破裂阈值强度函数；为灵敏度向量的转置；pk为协方差矩阵，更新公式，初值；为初值向量；该递推公式能够在线实时调整参数，无需存储全部历史数据；

152、s12.4 收敛判据

153、当连续3次事件后，参数变化量时，认为参数已收敛，可固定用于后续预测；

154、s12.5 验证试验方案

155、在实验室尺度，使用相似材料模型（水泥、石膏、云母粉模拟顶板）进行三点弯曲或均布加载试验；布置声发射探头监测裂纹扩展过程，同时测量载荷—位移曲线，计算实验的geff和λ，与理论预测对比；现场验证则选择两个地质条件相近的工作面，一个采用本判据预警并采取解危措施，另一个作为对照，统计冲击地压发生率。

156、与现有技术相比，本发明提出的一种厚硬顶板破断诱发冲击地压的临界判据方法具有以下技术效果：

157、1）物理机制明确，科学性更强：本发明不再简单依赖应力或变形的经验阈值，而是从厚硬顶板破断诱发冲击的本质——能量释放率的突变——出发，基于断裂力学建立了严格的数学模型，揭示了冲击发生的力学机制。

158、2）预警超前性与准确性高：本方法通过监测裂隙稳定扩展阶段的能量释放率变化趋势，并在破断瞬间的突变因子超过临界值之前即可做出判断，实现了从“事后监测”到“超前预警”的转变。同时，量化了动态增强效应，避免了误报和漏报。

159、3）参数物理意义清晰，适应性强：临界判据λc中的参数均可通过标准岩石力学试验获得，具有明确的物理意义。不同矿区的顶板条件差异直接体现在参数取值上，使得本判据方法具有良好的地质适应性和推广价值。

160、4）计算简便，易于现场实施：本发明提供了解析计算公式，现场技术人员仅需输入少量关键参数（顶板厚度、悬露长度、岩石强度等）即可快速计算λ值并与λc比较，无需复杂耗时的数值模拟，便于开发为便携式计算软件或集成到矿山安全监控系统中。

161、5）有效指导工程实践：本判据方法不仅能预警，还能指导卸压措施的效果评价。例如，采取爆破切顶后，顶板的有效悬露长度l减少，计算出的λ值会显著下降，当降至λc以下时，可认为解危措施有效，从而形成“监测-预警-解危-检验”的闭环控制。