基于混沌时间序列分析方法的矿山塌陷区范围预测

doi:10.11872/j.issn.1005-2518.2019.02.249

基于混沌时间序列分析方法的矿山塌陷区范围预测

曾俊晖^,, 李夕兵^,

中南大学资源与安全工程学院，湖南长沙 410083

Prediction of Mine Subsidence Area Based on Chaotic Time Series Analysis

ZENG Junhui^,, LI Xibing^,

School of Resources and Safety Engineering，Central South University，Changsha 410083，Hunan，China

收稿日期: 2018-01-23 修回日期: 2018-03-28 网络出版日期: 2019-04-29

基金资助:

国家自然科学基金重点项目“深部资源开采诱发岩体动力灾害机理与防控方法研究”. 编号：41630642

Received: 2018-01-23 Revised: 2018-03-28 Online: 2019-04-29

作者简介 About authors

曾俊晖（1993-），男，甘肃白银人，硕士研究生，从事采矿稳定性和岩层移动方面的研究工作shvezjh@163.com , E-mail：shvezjh@163.com

李夕兵(1962-)，男，湖南宁乡人，教授，从事硬岩开采及岩石动力学方面的研究工作xbli@csu.edu.cn , E-mail：xbli@csu.edu.cn

摘要

地下开采所引发的地表变形对矿山工业生产造成了严重危害，因此准确预测矿山塌陷区范围对于矿山安全生产具有重要意义。将混沌时间序列分析方法应用于矿山塌陷区范围预测中，以相空间重构理论为基础，采用小数据量算法计算得到时间序列的关键指标——最大Lyapunov指数，研究了塌陷区在相空间相点距离的演变规律，建立了塌陷区范围边界预测模型，并应用该模型对红岭铅锌矿塌陷区范围进行了分析预测。结果表明，矿山塌陷区范围变化具有混沌特性，时间序列分析方法能够很好地反映塌陷区范围的内在规律，通过计算得出红岭铅锌矿塌陷区范围时间序列的最大Lyapunov指数大于0，该矿山塌陷区范围的预测值与实际值基本吻合，误差大小不超过0.1%，验证了该方法的可靠性，为矿山塌陷区预测提供了一种新思路。

关键词： 地下采矿 ; 矿山监测 ; 混沌时间序列 ; 相空间重构 ; Lyapunov指数 ; 地表变形规律 ; 塌陷区预测

Abstract

With the advance of mining to deep in domestic underground metal mines,surface subsidence caused by mining has seriously threatened the industrial production of the mine and the ecological environment of the mining area.Therefore，accurately predicting the extent of mine subsidence area is of great significance to mine safety production and ecological protection of mining areas.Because the mine subsidence area is a complex system，it contains various random factors，and the evolution process of the system is often accompanied by the exchange of energy，showing nonlinear characteristics.This paper demonstrates the mine subsidence area as a nonlinear dissipative dynamic system.Combined with the time series analysis method in chaos theory，the range of mine subsidence area is predicted.The specific method studied in this paper is based on the phase space reconstruction theory，and the phase space reconstruction of the deformation time series of the mine subsidence area is carried out.A small data amount algorithm is used to calculate the key index of the time series-the largest Lyapunov exponent.According to the calculation results，the chaos of the mine subsidence zone system is identified.The evolution law and energy variation law of the phase distance of the collapse zone in the phase space are studied.The prediction model of the boundary of the subsidence zone is established.The model is used to analyze and predict the subsidence zone of the Hongling lead-zinc mine.The results show that the formation and expansion of the subsidence area of the Hongling lead-zinc mine is the result of the comprehensive effects of all kinds of factors such as geological conditions，mining engineering operations and inherent nonlinear characteristics.The maximum Lyapunov exponent and associated dimension of the time series of the subsidence area of the Hongling lead-zinc mine are calculated.It is verified that the variation of the subsidence area of Hongling lead-zinc mine has chaotic characteristics.The time series analysis method can well reflect the inherent law of the range change of the subsidence area.The chaotic characteristics of the different locations in the collapsed zone are different in the phase space due to different geological conditions and mechanical properties.The established mine subsidence area prediction system can better predict the change of the subsidence area.The predicted value is basically consistent with the actual value，and the error size does not exceed 0.1%，which verifies the reliability of the method.It provides a new idea for the prediction of mine subsidence area and guides mine safety production.

Keywords： underground mining ; mine monitoring ; chaotic time series ; phase space reconstruction ; Lyapunov index ; surface deformation law ; subsidence area prediction

PDF (5418KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

曾俊晖, 李夕兵. 基于混沌时间序列分析方法的矿山塌陷区范围预测[J]. 黄金科学技术, 2019, 27(2): 249-256 doi:10.11872/j.issn.1005-2518.2019.02.249

ZENG Junhui, LI Xibing. Prediction of Mine Subsidence Area Based on Chaotic Time Series Analysis[J]. Gold Science and Technology, 2019, 27(2): 249-256 doi:10.11872/j.issn.1005-2518.2019.02.249

随着浅部矿产资源的逐渐枯竭，越来越多的矿山已进入深部开采阶段，整体规模逐步扩大，地下开采所引发的地表变形和岩层位移对矿山工业生产造成了巨大影响。有关煤矿塌陷区的研究得到了长足发展，相关成果已经应用在矿山实践中。但是，金属矿山与煤矿的差异巨大，相比之下，金属矿山的地质条件更复杂，影响因素更多，因此对于金属矿山地下开采所带来的地表变形和岩层位移的研究难度很大，研究得到的规律性较差。塌陷区范围预测对于矿山有着十分重要的现实意义，若预测结果与实际值之间的偏差较大，则会影响矿山的安全生产，使矿山经济效益受损。因此，准确预测塌陷区范围是当前很多矿山亟待解决的难题之一。

地表塌陷是一种很严重的地质灾害，主要原因是采空区失稳无法支撑上覆岩石，导致地表大面积垮塌，其影响因素包括地质条件和采矿方法2个方面。研究发现，地质灾害的致灾机理具有明显的非线性动力学特征，根据已有文献资料^[1]可知，边坡是一个典型的非线性耗散动力系统，通过对比研究发现，塌陷区也是受地质条件控制的非线性耗散动力系统，且其受地质条件和采矿活动等众多因素的影响，比边坡更复杂。

通过时间序列来研究混沌行为由来已久，吕金虎^[2]提出的重构相空间理论开启了相关研究，并一直运用至今。对时间序列的分析已成为当前数据挖掘领域的研究热点之一，未来必将广泛应用于诸多领域，如金融、农业、工程和气象等领域。系统的混沌行为包含所有变量的长期演变，但显然对所有变量进行研究是不现实的，通常采取的方法是确定某一个单变量的特性来代替系统整体的混沌特性，该变量的特性通过时间序列得到。奇异吸引子作为研究对象，其本身指混沌系统的轨迹，且需要相空间的稳定，通常认为吸引子的不变量有若干个，选取常用的Lyapunov指数来研究整个系统的混沌量水平是合理的，因此，选取适当的方法计算得到Lyapunov指数对混沌系统的研究十分关键。复杂系统内部的关联和演变能够由外在的时间序列表现，在混沌系统中，一般来说奇异吸引子是相空间运动轨迹收缩得到的。

在混沌时间序列的研究中，当前使用的预测方法主要有神经网络法^[3]、小波网络法^[4]、遗传算法^[5]和最大Lyapunov指数预测方法，每种预测方法都有各自的特点和优势。根据塌陷区实际情况进行分析，神经网络法、小波网络法和遗传算法存在以下问题：神经网络法的缺点是所需数据量较大，小波网络法的不足是难以区别异常值，遗传算法仅仅搜索最优解。相比之下，最大Lyapunov指数预测方法具有计算简单、适用性强且更适用于混沌系统的特点，根据混沌系统的规律性和特殊性，基于监测塌陷区获得的历史数据，按照混沌时间序列分析方法，在相空间重构理论的基础上，计算得出Lyapunov指数，在适当的时间范围内研究矿山塌陷区的位移及塌陷情况，从而预测塌陷区范围的变化^[6]。

通过长期安全监测得到的塌陷区相关数据，运用时间序列预测方法建立的预警系统，能够克服以往主观分析结果的不确定性。长期以来矿山往往通过主观分析来判断塌陷区未来的范围，本文根据矿山塌陷区安全监测得到的相关数据采用时间序列预测方法建立的预警系统，能够克服以往主观分析结果的不确定性。本文选取合适的时间序列，采用Lyapunov指数预测方法，来预测矿山塌陷区的范围，进而分析该方法在塌陷区预测中的可行性。

1 基于Lyapunvo指数的预测方法

1.1 相空间重构

相空间重构的方法最早是由Ruelle^[7]提出的，后来经过不断发展，该方法更加完善。其基本原理是混沌系统经过不断演变最终会形成一个特定的轨迹，将这种轨迹称为吸引子，相空间重构的目的在于恢复高维相空间的混沌吸引子^[8]。相空间重构方法中系统的某一分量与其他分量之间是相互影响的，某一分量的演化进程可由其他分量确定，因此在一个分量的不断发展中包含了其他相关分量的内容，通过延迟处理等方法，将一些数值转变为新的观测值，从而构造出一个相空间与原来的空间等价，这样就能够研究其吸引子的性质^[9,10]，选择适当的动态模型（本文选择局域法）进行预测。

1.2 最大 Lyapunov指数计算

Lyapunov指数是判断系统是否具有混沌特性的指标，若其计算结果大于0则说明系统具有混沌性。同时，Lyapunov指数还能说明系统轨道的幅散程度，系统本身的可预测性问题则与最大Lyapunov 指数有关，它决定着初始信息的影响大小、邻近轨道的辐散程度和吸引子的发散快慢等，因此准确得到最大Lyapunov指数会有一个很好的预测结果。计算Lyapunov指数的方法有多种，由于本文数据量较小，对结果精度的要求较高，因此选取小数据量方法来计算最大Lyapunov指数。

（1）塌陷区变形序列重构相空间。由于塌陷区本身是个闭合的不规则曲线，将其近似为一个矩形，那么对于水平轴线来说上下边的位移属于单变量移动，选取了一些能够确定塌陷区范围的关键点，即x1，x2，…，x16。对于x1来说，设单变量的位移为{x1（ti），i=1，2，…，n}，通过G-P算法（遗传编程算法）得到关联维数d，嵌入维数m则可以通过m≥2d+1来确定。选择自相关函数达到极小值时确定时间间隔为λ=kd，可得到时间序列{x1（ti），i=1，2，…，n}的相空间表示为

\begin{array}{l} X_{1 i} (t) = \{x_{1} (t), x_{1} (t + λ), \dots, \\ x_{1} [t + (m - 1) λ]\} \end{array} i = 1,2, \dots, M

（1）

式中： $X_{1 i} (t)$ 为相空间中的点；M应具备的条件为 $M = n - (m - 1) k$ ；d为关联维数；m为嵌入维数；λ为时间延迟。

（2）小数据量算法计算最大Lyapunov指数。小数据量算法具有计算量少且快速准确等优点，十分适合本研究，其计算步骤如下：

①将该系统的时间序列嵌入欧式空间 $R$ 中，得到一个向量集： $X_{1 i} (t) = \{x_{1} (t), x_{1} (t + λ), \dots, x_{1} [t + (m - 1) λ]\}$ ，运用快速傅立叶变换方法，计算得到该组数据的时间延迟 $λ$ 和平均周期t；

②关联维数d由G-P方法计算得出，根据Takens嵌入定律计算系统的嵌入维数m，根据时间延迟λ和嵌入维数m重构相空间^[11]，即 $X_{1 j}, j = 1,2, \dots, M$ ；

③找到相空间中与 $X_{1 j}$ 相对应的最近邻点 $X_{1 \overset{⌢}{j}}$ ，将其分离开来^[12]，即：

d_{1 j} (0) = \underset{j}{m i n} ‖X_{1 j} - X_{1 \overset{⌢}{j}}‖, |j - \overset{⌢}{j}| > T

（2）

④相空间中的每个点 $X_{1 j}$ ，计算出该邻点对的第i个离散时间后的距离为 $d_{J} (i) = |X_{j + i} - X_{\overset{⌢}{j} + i} |$

i = 1,2, \cdot \cdot \cdot, m i n (M - j, M - \overset{⌢}{j})

（3）

对每个i求出对应的 $y (i)$ ：

y (i) = \frac{1}{p h} \overset{p}{\sum_{j = 1}} l n d_{j} (i)

（4）

式中： $p$ 为非零 $d_{j} (i)$ 的数目。

（3）Lyapunov指数预测模式及预测时间。设 $X_{M}$ 为预测模式的交点，可得到其在相空间中的最邻近点为 $X_{k}$ ，两点的距离为

d_{M} (0) = m i n ‖X_{M}⌋ - X_{j}‖ = ‖X_{M}⌋ - X_{k}‖

（5）

而已知最大Lyapunov指数为 $λ_{1}$ ，则系统的预报模式可表示为

‖X_{M}⌋ - X_{M + 1}‖ = ‖X_{k}⌋ - X_{k + 1}‖ e^{λ_{1}}

（6）

式中：e为自然常数。

根据上式可知，点 $X_{M + 1}$ 只有最后一个分量 $x (t_{n + 1})$ 是未知数，因此 $x (t_{n + 1})$ 是可计算得到的，是可预测的。

通常，将最长的预报时间定为

T_{m} = \frac{1}{λ_{1}}

（7）

$T_{m}$ 能够说明系统误差增加一倍时系统所需的时间，因此可成为短期系统的可靠性预测指标之一^[13]。

2 工程实例

2.1 工程地质和水文地质概况

红岭铅锌矿位于内蒙古自治区赤峰市巴林左旗林东镇乌兰达坝苏木浩布高嗄查（村）。矿区全长 5 700 m，矿体赋存于标高630～1 110 m之间，最大延深480 m，矿体总体走向NE59°，倾向NW，倾角为60°～80°，平均倾角为75°，矿体走向长度为1 350 m。主要矿体为磁铁矿和铅锌矿，矿区地质构造主要是NE向构造，与区域构造方向基本一致。

红岭铅锌矿属于中山地形，坡度较陡，海拔为890～1 360 m，附近没有地表水体，在北西1 500 m处有浩布高小河，当地最低侵蚀基准面标高为 890 m，位于浩布高小河的下游。区内主要含水层是第四系孔隙潜水含水层和基岩裂隙含水层。矿坑主要充水因素为基岩裂隙含水层的地下水，通过构造裂隙和断层等直接充水方式进行充水，水文地质条件属于简单类型。

2.2 地表塌陷情况

红岭铅锌矿采用的开采方法是崩落法转空场法，故矿区内已形成了一定范围的塌陷区，如图1所示。通过塌陷区周界的实测发现，目前塌陷区最深约320 m，长度100 m，宽度50 m，周边可见拉裂缝最长在6 m以上，深度达2 m，如图2所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 红岭铅锌矿部分塌陷区

Fig.1 Partial subsidence area of Hongling lead zinc mine

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 红岭铅锌矿塌陷区拉裂缝

Fig.2 Pull cracks in the subsidence area of Hongling lead zinc mine

2.3 塌陷区动力系统维数计算

红岭铅锌矿塌陷区范围监测工作始于2015年1月，近2年的监测数据为寻找规律提供了丰富的资料。由于塌陷区存在近似对称性，并根据矿体上下盘可划分为塌陷区上部和塌陷区下部，在A₁~A₁₆16个监测点中选择了A₂、A₇、A₉和A₁₅共4个监测点进行分析，其中监测点A₂和A₇处于塌陷区上部，监测点A₉和A₁₅处于塌陷区下部，同时监测点A₂和A₉的经度相同，监测点A₇和A₁₅的经度相同。根据式（1）采用重构相空间方法，并结合G-P方法^[14]计算出关联维数d₂、d₇、d₉和d₁₅分别为1.3965、1.3876、1.2738和1.2854，根据Takens嵌入定律^[15]计算得到A₂、A₇、A₉和A₁₅的嵌入维数均为4。由于A₂、A₇、A₉和A₁₅的关联维数均不是整数，表明系统具有混沌特征。

2.4 不同监测点塌陷区变形在相空间中的演变规律

根据相空间重构理论的观点，相空间中的每一个相点都对应系统的一个状态^[16]，塌陷区系统重构意味着对其变形时间序列重构，其相点和状态也是一一对应的，需要计算任意相点最近邻相点的距离，应满足条件：

Z = m i n (‖X_{1} - X_{i}‖)

，

i = 1,2, \dots, 8

（8）

若相点 $X_{k}$ 的最邻近相点为 $X_{b}$ ，则 $X_{k}$ 与 $X_{b}$ 之间的距离为

Z_{k} = {\{{[x (t_{k}) - x (t_{b})]}^{2} + {[x (t_{k} - 1) - x (t_{b} - 1)]}^{2} + \dots + {[x (t_{k} - 4) - x (t_{b} - 4)]}^{2}\}}^{\frac{1}{2}}

（9）

根据式（8），在塌陷区系统重构相空间后，获得了系统的最邻近点距离的演变序列Z₂、Z₇、Z₉和Z₁₅，并得到了该序列的曲线随时间的变化情况，如图3所示。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 塌陷区变形的相空间距离演变曲线

Fig.3 Phase space distance evolution curve of deformation in subsidence area

由图3可知：当d $<$ 350时，Z₂ $>$ Z₉，Z₇ $>$ Z₁₅，表明塌陷区上部变形量大于塌陷区下部，当d $>$ 350时，Z₂ $<$ Z₉，Z₇ $<$ Z₁₅，表明塌陷区上部变形量小于塌陷区下部；同时观察曲线可知，Z₉和Z₁₅的起伏程度明显大于Z₂和Z₇，表明监测点A₉和A₁₅的变形过程比监测点A₂和A₇更加复杂。

根据以上分析得出，采用相空间重构方法之后，系统原有的细微变化特征能够得到充分展示，能够将细小的变化放大。若采用传统的概率统计方法，则不能发现监测点A₂、A₇、A₉和A₁₅变形过程的区别，也无法得知监测点A₂、A₇、A₉和A₁₅在d=350时变形量发生大小转换，这无疑会导致研究结果出现较大误差。

2.5 塌陷区变形时频分析

对塌陷区系统进行能量分析，通常采用的方法是将序列Z₂、Z₇、Z₉和Z₁₅进行傅立叶变换^[17]，用 $Z_{j}$ 表示 $Z$ 的时间功率，根据式（10）计算得到监测点A₂、A₇、A₉和A₁₅的功率谱曲线，如图4所示，该曲线表示不同监测点的变形能量随时间的变化规律^[18]。

Z_{j} = \frac{1}{\sqrt[]{n - m + 1}} \sum_{j = 1}^{n - m + 1} Z_{k} e^{\frac{2 π j k}{n - m + 1}}

（10）

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 塌陷区变形的时间—功率谱曲线

Fig.4 Time power spectrum curve of deformation in subsidence area

式中： $n = 92, m = 5, k = 1,2, \dots, 88$ 。

由图3可知，监测点A₂与A₇、A₉与A₁₅所处位置的地质条件相似，所以其相空间距离演变曲线相近，起伏相同；而监测点A₂、A₇和A₉、A₁₅所处位置的地质条件不同，所以其相空间距离演变曲线差异明显，因监测点A₉和A₁₅的地质条件更为复杂，相比之下其相空间距离起伏更大。由图4可知，监测点A₂和A₇曲线所形成的面积稍大于监测点A₉和A₁₅，说明监测点A₂和A₇所在位置承受破坏的能力稍弱。

以上分析表明：在塌陷区的不同位置，其变形表现出不同的动力学特征，地质条件、外部随机影响因素等共同决定着塌陷区变形的动力学行为。

2.6 塌陷区范围预测

典型监测点A₂和A₉在2015年1月至2017年12月的累计位移曲线如图5所示，取 $h = 30 d$ ，得到 $h$ 时间段内的净时间位移序列 $\{u (t_{1}), i = 1,2, \cdot \cdot \cdot, 8\}$ 。现对8个该时序进行洛伦茨变换，得到平均周期 $T = 1 000 d$ ，取 $m = 4$ ， $τ = 2 h$ ，通过小数据量方法计算得到最大Lyapunov指数大于0，表明该序列为混沌序列。为了更精确地得到Lyapunov指数，通过数据累加的方法减弱了监测数据的随机性对计算结果的干扰。采用小数据量算法^[19,20]分别得到 $λ_{1} = 0.00301 >$ 0、 $λ_{2} = 0.0074 >$ 0。根据式（7）计算的结果进行预测，预测步数为40^[21]，预测数据与实际监测数据对比情况见图5、图6和表1,2，预测数据与实际监测数据误差较小，获得了较理想的预测结果。通过上述方法计算出16个典型观测点的预测值，最终得到塌陷区预测范围与实际范围对比情况如图7所示。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 实测位移曲线

Fig.5 Measured displacement curves

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 塌陷区预测范围与实际范围比较

Fig.7 Comparison between prediction and actual scope in subsidence area

表1 A₂监测点位移预测结果

Table 1 Prediction results of displacement at A₂ monitoring point

时间/d	实测值/m	预测值/m	相对误差/%
390	29.42	29.61	0.65
420	31.43	31.68	0.80
450	35.06	35.01	0.14
480	40.85	41.09	0.59
510	49.88	50.21	0.66
540	57.02	56.72	0.53
570	66.14	66.57	0.65
600	75.05	74.45	0.80

新窗口打开| 下载CSV

表2 A₉监测点位移预测结果

Table 2 Prediction results of displacement at A₉ monitoring point

时间/d	实测值/m	预测值/m	相对误差/%
390	25.06	25.25	0.76
420	26.08	26.32	0.92
450	28.71	28.65	0.21
480	33.49	33.73	0.72
510	41.55	41.86	0.75
540	48.62	48.36	0.53
570	59.83	59.31	0.87
600	71.69	71.11	0.81

新窗口打开| 下载CSV

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 实测与预测位移曲线比较

Fig.6 Comparison curves of measured and predicted displacement

3 结论

本文将塌陷区引申论证为一个非线性耗散动力系统，结合相空间重构理论，采用Lyapunov指数构造预测模式进行塌陷区范围数据的预测，得出以下结论：

（1）在矿山生产中，塌陷区的形成及扩大是其地质条件、采矿工程作业及内在非线性特性等多种因素综合作用的结果，研究发现该系统具有混沌特征，在相空间重构理论基础上，综合分析得到系统混沌特性的内在基本规律。

（2）通过研究不同位置、不同地质条件下的监测点变形情况，发现其表现出不同的力学特征，同时在相空间中也表现出不同的混沌特性。

（3）研究结果表明，重构相空间方法能够对塌陷区范围大小进行短期预测，在红岭铅锌矿塌陷区范围预测中取得了良好的结果，计算得出红岭铅锌矿塌陷区范围时间序列的最大Lyapunov指数大于0，该塌陷区范围的预测值与实际值基本吻合，误差大小不超过0.1%。这是研究塌陷区混沌现象内在规律的一次有益尝试，为地表变形范围的圈定和矿山安全生产提供了指导。

本文将混沌时间序列方法用于矿山塌陷区变形研究及范围预测尚处于初步阶段，还需开展进一步深入研究，下一步将通过大量工程实例验证该方法的通用性，并尝试不同算法改善其精确性和可靠性。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Packard

N H

，Crutchfield

，Farmer

，et al.

Geometry from a time series

［J］.Physical Review Letters，1980，45（9）：712.