基于多变量监测时序的冲击地压复杂性分析

陶 慧，李 莹，马小平

（1.河南理工大学电气工程与自动化学院，河南 焦作 454000；

2.中国矿业大学信息与电气工程学院，江苏徐州 2211162；3.焦作大学计算机学院，河南 焦作 454000）

基于多变量监测时序的冲击地压复杂性分析

陶 慧1，2，李 莹3，马小平2

（1.河南理工大学电气工程与自动化学院，河南 焦作 454000；

2.中国矿业大学信息与电气工程学院，江苏徐州 2211162；3.焦作大学计算机学院，河南 焦作 454000）

论文基于多变量时间序列相空间重构来计算数据的关联维数，以研究冲击地压监测数据的复杂程度。考虑到冲击地压监测数据含有噪声而且长度有限，对传统G-P算法进行了扩展改进，给出了改进算法求解多变量时间序列的关联维数的原理，并用于Lorenz混沌系统检验了改进算法的有效性。收集了不同冲击情况下多种监测类型的冲击地压时间序列数据，用改进G-P算法求解这些监测数据的关联维数值。研究结果表明：冲击地压监测数据具有混沌特性，而且数据关联维数越大，复杂程度越高，对应矿井的冲击破坏性越强。这为基于混沌理论预测冲击危险性提供了新方法和依据。

冲击地压；多变量时间序列；关联维数；混沌特性

冲击地压是矿井生产中最严重的自然灾害之一，其破坏可比拟为一类非线性微分方程的倍周期分叉而出现的混沌运动现象，谢和平、尹光志等已证实其动力学行为具有混沌特性［1－2］。冲击破坏过程十分复杂，很难建立精确的数学模型，监测方法有微震、电磁辐射、钻屑法等［3－5］，可以从现场获取多种冲击地压监测时间序列，通过观测时间序列来分析其复杂混沌特性具有重要意义。描述系统混沌特性的几何不变量主要指关联维数和Lyapunov指数。其中关联维数d2是系统复杂性的度量，大于关联维数的下一个整数是刻画系统所需的独立变量的个数，为从时间序列恢复原始系统确定了一个框架。大量学者对冲击地压监测时间序列的关联维数d2进行了分析计算，以描述数据的复杂度。如文献［6］对不同尺度的电磁辐射数据的混沌特性进行了分析，文献［7］研究了微震时间序列的关联维数等混沌特性。

但目前对冲击地压监测数据的混沌特性分析都基于单变量时间序列。根据F.Takens的嵌入理论，在数据无限长、无噪声的理想条件下，单变量时间序列足以重构原动力系统［8］，但实际数据一般并不能满足理想条件，实践中并不能保证任何单变量时间序列都能很好的实现重构。另外，由于冲击地压演化过程复杂，用单一指标进行识别本身就存在很多不足［9］，多变量时序重构包含更加丰富的系统信息，可以弥补数据长度的不足，克服噪声的影响，因此基于多变量时序重构的关联维数计算结果更加真实准确［10］。

最常用计算关联维数d2的方法是G-P算法［11］，一般只针对单变量时间序列相空间重构。本文将对G-P算法进行扩展改进后求解多变量时序重构的关联维数d2；接下来将改进的G-P算法用于标准混沌系统以说明其有效性，最后将其用于求解冲击地压监测时序的关联维数d2，以描述数据的复杂程度。

1 改进G-P算法计算多变量时序的关联维数

关联维数是系统复杂性的一个很好的度量，混沌系统的关联维数为分数维数，一般认为，大于关联维数的下一个整数是刻画系统所需的独立变量的个数，这为从时间序列恢复原始系统确定了一个框架［8，10］。对单变量时间序列，最常用的关联维数计算方法是 G-P算法［11］，本文将 G-P算法进行扩展以求解多变量时序相空间重构的关联维数d2。

对长度为N的有限时间序列，重构的状态变量个数Nm＝N－N0＋1，是有限的，而关联积分的值只有在N充分大的情况下接近理论值。常用的计算关联维数方法是让嵌入维数m变化，来观测双对数lnCN（r）-lnr图的斜率。当m＜d2时，由于重构相空间的维数不够高，一般不能分辨吸引子结构，此时lnCN（r）-lnr图的斜率近似为嵌入维数；当m增加到m＞d2时，吸引子的分辨能力改善，斜率达到一个平稳时期，在双对数lnCN（r）-lnr图中找到平行直线区域，这些直线的斜率与m的依赖性消失，该斜率即为关联维数d［12］。

2

考虑到多变量重构时，嵌入维数m增加对应于多种mi的变化，本文将依次增加各个重构变量的嵌入维数，让嵌入维数m在一定的范围内变化，来确定双对数lnCN（r）-lnr图的线性部分，再通过斜率计算对应的关联维数。在实际应用中，确定了线性部分后，为了消除计算误差，一般采用最小二乘法来拟合双对数lnCN（r）-lnr图的线性部分的斜率。

取双对数线的线性区段内的K个点r1，r2，…rk，由最小二乘法得到式（6）。

2 改进G-P算法应用实例

下面先通过一个具体的例子来说明改进的G－P算法求解多变量时间序列关联维数的计算效果。对文献［12］的Lorenz系统，取数据长度为3000并进行归一化，分别利用第一个变量或者前两个变量进行相空间重构，并在原始数据的基础上添加5%的噪声，得到的双对数lnCN（r）-lnr图见图1。

在图1中，单变量重构时，重构参数m＝2∶1∶10，τ＝3。图1（a）表示采用无噪声x1进行重构的结果；图1（b）为附加5%噪声的x1进行重构的结果；双变量重构时，令τ1＝τ2＝3，从m1＝2，m2＝1依次增加m1和m2使m＝3∶2∶19；图1（c）表示采用无噪声的x1和x2进行双变量重构的结果；图1（d）表示采用附加5%噪声的x1和无噪声的x2进行双变量重构的结果。

图1 求解Lorenz的关联维数的lnCN（r）-lnr图

分析图1可以看出，即使数据含有一定的噪声，在相同的嵌入维数下，双对数曲线的变化也不大，能比较容易找到双对数曲线的线性部分。在m＝5时取20个点进行最小二乘拟合，求得四种重构方式下对应的关联维数如表1所示。根据参考文献［13］，将d2＝2.04作为Lorenz系统的理论关联维数值。由于本文数据长度较短，采用单变量重构求得的嵌入维数有一定误差，但采用多变量时间序列重构时计算结果非常精确，但本文所用数据的长度比参考文献［13］中要短得多；另外，噪声的存在使求得的关联维数d2的值略有增加，但增加不大。

表1 Lorenz系统的关联维数

总之，采用改进的G-P算法估计多变量时序的关联维数，对噪声干扰的鲁棒性好，对数据长度要求比单变量时序低，能够方便且准确地求出基于多变量时序重构的混沌系统的关联维数，描述系统的复杂程度。

3 冲击地压数据来源

为了全面分析冲击地压监测数据的混沌特性，本文收集了不同冲击情况下的多种监测数据。

3.1 微震时间序列数据

煤岩体发生破裂、滑动过程中，会产生一定频率的向周围传播的地震波，利用微震监测系统，可确定各种震动参数。微震事件与冲击地压的孕育和发生密不可分［3］。微震能量释放值能够表现出微震活动的阶段性和微震系统整体特征指标，因此将一段时间区间的累计微震能量指标作为一个微震监测变量。另外，微震时间分布具有层次结构和自相似性，因此将一段时间区间的微震频次指标，作为第二个微震监测变量。根据微震监控系统我们可以获得累计能量、累计频次2个微震监测时间序列。需要说明的是，由于微震能量的范围从101到108，变化幅度较大，为了更好的反映其变化规律，对其做了取对数处理。

本文微震监控时间序列来源于3个矿井，分别为平煤十一矿，义马千秋煤矿、义马跃进煤矿。平煤十一矿由于存在多个孤岛工作面，局部应力高度集中，从2008年至今发生40余次冲击地压事故，但微震能量都不高，冲击破坏性很小，本文收集了该矿2009／10／01至2012／12／30期间每天微震累计能量和累计频次两个时间序列作为第一组微震数据。义马千秋煤矿主采煤层为侏罗系二煤，具有中等冲击倾向性，基本顶为厚度205m的巨厚砾岩层，该矿井的冲击地压威胁十分严重，共计发生破坏性冲击地压20多次，其中破坏巷道超过60m的冲击次数占60%以上，相比于平煤十一矿，冲击破坏性更强。本文整理得到的2010／04／06至2012／09／23每天累计能量和频次微震时间序列作为第二组微震数据。义马跃进煤矿与千秋煤矿情况相似，至2004年首次发生冲击地压以来，该矿已发生冲击地压30多次，整理的2010／06／01到2012／9／30日的每天累计能量、累计频次两个时间序列作为第3组微震数据。

3.2 电磁辐射时间序列

煤岩体受载发生破裂，因裂缝的形成及扩展过程中电荷的迁移及煤岩颗粒的摩擦而产生电磁辐射。电磁辐射的变化反映了煤岩破坏发生发展的过程，电磁辐射强度（幅值）反映了煤岩体受载及变形破裂强度，脉冲数反映了煤体应力集中程度［4］。一般采用KBD5型便携式电磁辐射监测仪，在工作面及其上巷、下巷等不同位置每隔15m或20m设置一个测点，一个工作面会设置十几至几十个测点，来测量不同位置的应力情况。可以收集每天整个工作面的电磁辐射最大强度、平均强度和脉冲数3个时间序列，来反映其应力变化。

本文电磁辐射数据来源于跃进煤矿23070工作面。该工作面在2011年2月到2012年1月进行上、巷掘进，2012年2月至12月进行切眼掘进，在切眼掘进期间共发生13次冲击地压，而上下巷掘进期间基本未发生冲击事故。由于不同时期煤岩系统的动力学特性发生了变化，因此将其数据分成两组，第一组为2011／02／25至2012／01／31上下巷掘进期间的电磁辐射数据，第二组为2012／02／01至2012／12／31切眼掘进期间的数据。

3.3 结合多种监测方式获得的监测数据

跃进煤矿25110工作面采深达1040m，所采煤层具有强冲击倾向性，直接顶为厚30m块状易碎的泥岩，老顶为砾岩，采用综放开采方式易导致容易导致应力集中。该工作面冲击事故多发，仅回采期间就发生10余次冲击破坏事故。跃进煤矿安装有ARAMIS微震监测系统，该工作面微震数据完整。另外采用KBD5电磁辐射仪监测获得了该工作面的电磁辐射数据［5］。作者收集了 2011／07／01 至2012／10／31该工作面回采期间的微震累计能量和频次、电磁辐射最大强度、平均强度和脉冲数5个时间序列数据。

跃进煤矿23130工作面顶板覆400余米巨厚砾岩，也是一个冲击地压事故多发的工作面。本文收集了该工作面2008／10／05至2009／03／20的每天电磁辐射最大强度、平均强度和脉冲数，和0～5m平均煤粉量、6～10m平均煤粉量5个时间序列。

4 冲击地压数据的d2计算及结果分析

4.1 微震时间序列的d2计算及结果分析

首先求解平煤十一矿和千秋煤矿微震时间序列的关联维数。对原始数据归一化后，采用微震累计能量、频次2个变量进行相空间重构，重构参数选择m1＝3，m2＝2，τ1＝1，τ2＝2，然后依次增加m1和m2使m＝5：2：19增加，通过改进的G-P算法得到的十一矿数据的双对数lnCN（r）-lnr见图2（a）。同样的，得到千秋煤矿微震数据的lnCN（r）-lnr曲线，见图2（b）。

图2 微震数据的关联维数求解

在m＝9时，取图2（a）的线性段lnr＝－0.425～－0.153和图2（b）的线性段lnr＝－0.362 ～－0.113中的20个点进行最小二乘拟合，得到两矿微震数据的关联维数d2分别为2.48和3.92。也就是说，至少需要3个变量才能描述平煤十一矿的微震监控信息，4个变量才能描述跃进矿的微震监控信息，跃进矿的监控信息更加复杂。

在多变量微震时间序列相空间重构的基础上，利用改进的G-P算法计算得到的跃进煤矿及其25110工作面的关联维数见表2。可以看出，多变量微震时间序列的关联维数d2均是分数维数，定性说明微震时间序列具有混沌特性。

比较表2数据，发现平煤十一矿的关联维数明显小于千秋和跃进煤矿的结果，表明千秋和跃进煤矿的监测数据复杂程度要强于平煤十一矿。实际上，在监测期间，平煤十一矿发生冲击事故仅十余次，且破坏性都不大；千秋煤矿和跃进煤矿发生的破坏性冲击地压达20多次，其中多次事故造成严重的损失和人员伤亡。对比跃进煤矿及其25110工作面的关联维数发现，全矿微震监测数据的关联维数明显大于其下属工作面的结果，这说明全矿范围内的冲击地压孕育发生过程的复杂程度明显强于其子工作面。因此，可以得出结论，破坏性冲击事故越多，破坏程度越强，其关联维数越大。

表2 多变量微震时间序列的关联维数

4.2 电磁辐射数据的d2计算及结果分析

23070的工作面的两组电磁辐射数据的关联维数时，采用电磁辐射最大强度、平均强度和脉冲数三个变量进行相空间重构后，改进的G-P算法求得两组数据的关联维数d2分别为2.43和4.03，都是分数维数，这说明电磁辐射时序也具有混沌特性。对同一个工作面，由于冲击情况不同，我们看到计算出的关联维数有明显差别，上下巷掘进期间，基本没有冲击危险，第一组数据的复杂程度不高，关联维数较小。而在切眼掘进期间，发生了10余次冲击地压破坏事故，明显第二组数据的更加复杂，对应的关联维数较大。23070工作面两组数据的结果再次说明，冲击地压破坏事故次数越多，破坏程度越强，其关联维数越大。

4.3 结合多种监测数据的d2计算及结果分析

25110工作面的监测数据包括微震和电磁辐射数据，采用2个微震监测变量重构可以计算微震时序的关联维数；采用3个电磁辐射时序重构后，得到辐射时序的关联维数；采用微震和电磁辐射5个时序重构后，得到所有监测时序的关联维数。结果见表3。类似的，23130工作面的监测时序包括电磁辐射和转屑量两类数据，仅采用电磁辐射或钻屑量一类数据和结合两类数据计算的关联维数也见表3。

表3 多类监测数据的关联维数

比较表3结果我们发现，对同一工作面，不同类型数据计算得到的关联维数差别不大，都是分数维数。这说明冲击地压监测时序具有混沌特性，而关联维数基本一致的原因是这些监测数据都是该工作面冲击地压孕育发展和产生过程中各种信息的反映，同一工作面其过程相同，因此数据复杂程度接近，对应的关联维数也接近。

5 总结

基于多变量重构理论，对求解单变量时序关联维数的G-P算法进行扩展和改进以求解多变量时序的关联维数。以Lorenz系统为例，发现改进G－P算法求解多变量时序的关联维数时，对数据长度的要求降低、鲁棒性强、求解结果更加准确。求得冲击地压监测时序的关联维数均为分数维，说明其具有混沌特性。另外，不同监测数据的关联维数表明，冲击地压事故破坏性越强，其关联维数越大，监测数据的复杂程度越高。同一工作面的不同种类监测数据的关联维数接近，说明反映同一工作面冲击地压孕育发展、产生过程中各类数据复杂程度接近。这些结论说明可以利用计算的关联维数值间接判断冲击危险性，为冲击危险预测提供了一种新方法。另外，由于证明监测数据具有混沌特性，这为利用混沌理论预测冲击危险性提供了依据。

［1］ Xie H P，Chen Z D.Fractal geometry and fracture of rock［J］.Acta Mechanica Sinica，1988，4（3）：255-264.

［2］ 黄滚，尹光志.冲击地压粘滑失稳的混沌特性［J］.重庆大学学报，2009，2（6）：633-637.

［3］ 夏永学，潘俊锋，王元杰，等.基于高精度微震监测的煤岩破裂与应力分布特征研究［J］.煤炭学报，2011，36（2）：239-243.

［4］ 陈国祥，窦林名，曹安业，等.电磁辐射法评定冲击矿压危险等级及应用［J］.煤炭学报，2008，33（8）：866-870.

［5］ 陆振裕，窦林名，徐学锋，等.钻屑法探测巷道围岩应力及预测冲击危险新探究［J］.煤炭工程，2011（1）：71-74.

［6］ 刘贞堂，赵恩来，王恩元，等.不同尺度电磁辐射时间序列的混沌特征初步分析［J］.煤炭学报，2009，34（2）：224-227.

［7］ 陶慧，马小平，乔美英.基于微震时间序列的冲击地压混沌特性分析［J］.煤炭安全，2012，3：52-55.

［8］ Takens F.On the numerical determination of the dimension of an attractor［J］.Lecture notes in Mathematics，1981（898）：361-381.

［9］ 赵善坤，李宏艳，刘军，等.深部冲击危险矿井多参量预测预报及解危技术研究［J］.煤炭学报，2011，36（S2）：339-345.

［10］ 王海燕，卢山.非线性时间序列分析及其应用［M］.北京：科学出版社，2006.

［11］ Grassberger P，Procaccia I.Characterization of strange attractors［J］.Physical review letters，1983，50（5）：346-49.

［12］ 王海燕，盛昭瀚，张进.多变量时间序列复杂系统的相空间重构［J］.东南大学学报：自然科学版，2003，33（1）：115-118.

［13］ Lai Y C，Lemer D.Efective scaling regime for computing the correlation dimension from chaotic time series［J］.Physica D，1998（115）：1-18.

Rock burst complexity analysis based on multivariate monitoring time series

TAO Hui1，2，LI Ying3，MA Xiao-ping2
（1.School of Electrical Engineering and Automation，Henan Polytechnic University，Jiaozuo 454000，China；
2.School of Information and Electrical Engineering，China University of Mining and Technology，Xuzhou 221116，China；
3.School of Information Engineering，Jiaozuo University，Jiaozuo 454000，China）

This paper studies the complexity of the data to monitor Rock burst through computing correlation dimension based on phase-space reconstruction of multivariate time series.Given that Rock burst monitoring data had limited-length and contained noise，traditional G-P algorithm is extended and improved.The principle of improved G-P algorithm of solving correlation dimension of multivariate time series was provided，and algorithm was verified through employing it to Lorenz chaotic system.Then a mass of timeseries data to monitor Rock burst were collected by diverse equipment under different burst degree，and their correlation dimensions were computed through improved G-P algorithm.The results demonstrate that the data have chaotic characteristic，and the larger correlation dimension is，the more complex monitoring data is，the stronger Rock-burst damage of corresponding coal mine is.Our achievement can give a novelty approach and basis to predict Rock burst risk based on chaos method.

rock burst；multivariate time-series；correlation dimension；chaotic characteristic

TD324，O415.5

A

1004-4051（2015）10-0118-05

2014-04-23

国家自然科学基金项目“基于多传感器信息融合技术的矿井瓦斯突出预警研究”资助（编号：60974126）

陶慧（1979－），女，汉族，河南信阳人，副教授，博士，主要研究领域为混沌预测和智能计算。E-mail：65724138＠qq.com。