时间序列ARIMA模型在煤矿地音监测系统中的应用

2023-11-04孙学波王元杰陈法兵

煤炭工程 2023年10期

孙学波，刘宁，王元杰，陈法兵，李岩

(中煤科工开采研究院有限公司，北京 100013)

在煤炭资源的深入开采过程中，因开采导致煤岩结构变化进而导致一系列的煤岩动力灾害，其中冲击地压带来的灾害后果尤为突出，给煤矿安全生产与人民的生命财产造成严重威胁。到2021年底，全国约有138座煤矿被鉴定为冲击地压矿井，对这一类型煤矿在开采过程中需要引入硬件装备进行监测预警，以降低冲击地压等动力灾害的发生[1，2]。经过科技工作者多年的技术攻关，针对冲击地压的监测预警也形成了多种类型的监测手段与方法，对冲击地压的监测预防起到了一定的效果。比较典型的技术手段有矿压观察法、煤岩应力监测法、电磁辐射监测法、钻屑监测法、地音监测法、微震监测法等。相比于微震监测法[3]，地音监测的振动事件属于低能高频事件，频段一般处于200～3000 Hz，微震法监测的是高能低频事件，频段一般处于0～200 Hz，两者侧重于监测不同的频段区域。地音现象的产生是由于煤岩应力的变化使煤岩体产生弱小微应变(如裂隙产生)，同时微应变以弹性波的形式向外释放而产生的现象。经过大量实际数据表明当有较强煤岩能量释放前，都能收到地音信号，通过对地音信号的有效监测，可以提高对冲击地压的预警作用[4]。

国内外科技工作者对地音信号在煤矿冲击地压监测方面做了大量研究，韩玉鉴从能量释放方式角度深入分析了煤岩体破坏过程中的几种演化可能，再通过计算地音事件加权平均能量值对地音事件进行分析。姚明远借助于冲击启动理论，提出了根据不同的冲击启动及冲击显现的位置，将地音预警信号划分为三种类型，并对三种类型进行细化，研究地音监测技术在预测煤矿冲击地压的地音信号特点。笔者通过应用时间序列ARIMA模型研究地音数据中的潜在规律，提取监测区域煤岩活动的相关信息，进而推断被监测区域未来的发展变化，为冲击地压的监测预防提供有力支撑。本文分析了时间序列ARIMA模块的原理与工作流程，同时在煤矿现场安装一套ARES-5/E监测系统并实施连续监测，通过对地音多通道监测数据进行数据清洗，剔除随机与干扰数据，然后对数据进行平稳性检验与数据差分分析，接着建立ARIMA数据模型并对模型进行相关性估计，识别模型参数后再对模型的有效性进行检验，通过检验后的模型对地音数据进行预测与识别。通过对某矿地音设备监测的长期数据进行多层次提取并预测后一时间段的地音监测值，发现时间序列ARIMA模型在中短期地音数据预测方面具有显著优势[5-9]。

1 时间序列ARIMA理论基础

1.1 ARIMA模型

ARIMA模型全称为自回归移动平均模型(Autoregressive Integrated Moving Average Model，简称ARIMA)[10]。

AR为自回归模型，该模型描述当前数据和历史数据之间的关系，用变量的历史数据对未来的数据进行预测，不受模型变量相互独立的假设条件约束，若要应用AR自回归模型则要求数据具有平稳性特点。

(1)

式中，yt为当前参与计算的数据；μ为模型的常数项；p为自回归项数；γi为数据自相关系数；εt为误差系数。

MA为移动平均模型，用过去各个时期的随机干扰或预测误差的线性组合来表达当前预测值，关注的是模型中的误差项的累加。

(2)

式中，q为滑动平均项数；θi为误差自相关系数。

ARIMA(p，d，q)为自回归移动平均模型，模型中参数p、q为上述的自回归与滑动平均项数，d为使数列成为平稳数列的差分阶数。该模型经过差分转换把数据转换成平稳数列，再将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型，ARIMA模型是在ARMA模型的基础上加上数据差分化[11]。

(3)

1.2 ARIMA模型数列平稳性与白噪音检验

ARIMA模型需要数据序列是平稳非白噪声型序列，在使用该模型进行数据预测前需要对两者进行检验。若通不过检验则需要进行数据差分化以满足平稳性等相关要求。

平稳性检验一般有三种方法：一是观察序列时序图是否存在波动或向上向下的趋势；二是对数列进行自相关分析，若随时间间隔增大，自相关系数衰减率高则是平稳序列，否则是非平稳序列；三是对数列进行单位根检验，p<α时，则拒绝存在单位根的原假设，序列平稳[11-15]。

1.3 ARIMA模型定阶

数列通过了平稳性与白噪音检验后，需要确定模型中p、q参数。确定p、q参数有多种方法，具体如下：

1)通过分析数列的自相关函数图(ACF)和偏自相关函数图(PACF)。

2)通过信息准则AIC和BIC来确定p、q值。

赤池信息准则AIC：是一种衡量统计模型拟合优良性的标准，在一组可供选择的参数值中一般选择AIC值最小的[16]。

AIC=2k-2ln(L)

(6)

式中，k为模型中参数个数；L为极大似然函数值。

贝叶斯信息准则BIC：与AIC类似，也是一种判别统计模型拟合优良性的标准。

BIC=kln(n)-2ln(L)

(7)

式中，n为样本数量。

AIC和BIC相比较，后半部分一样，当n较大时kln(n)≥2k，所以一般BIC用于选择参数较少的模型[17]。

1.4 ARIMA模型检验

定阶后的ARIMA是否满足数列处理需要，需要验证所选择参数的模型残差是否为白噪声，若模型残差不属于白噪声，说明选择的参数或模型还不是最优的。ARIMA模型残差的白噪声检验常用的方法有以下几种：

1)分析残差是否是符合均值为0、方差为常数的正态分布，同时也观察连续残差自相关系数的相关性。

2)分析数列的自相关函数图(ACF)，若所有系数值在置信区间内，则通过白噪声检验。

3)NormalQ-Q：通过Q-Q图判断残差的正态性，若基本和蓝线重合，头尾部可以有偏移，只要头尾部的outlier不是过分大，则满足正态性。

4)德宾-沃森(D-W)检验，称D-W检验，是目前检验自相关性最常用的方法，但它只使用于检验一阶自相关性。当DW值显著的接近于2时，则不存在(一阶)自相关性。

式中，DW为D-W检验量；ut为t时刻序列值；ut-1为t-1时刻序列值。

5)Ljung-Box检验，是对时间序列是否存在滞后相关的一种统计检验，检验的结果一般是看最后一列前十二行的检验概率，如果检验概率小于给定的显著性水平就拒绝原假设，其原假设是相关系数为零。

1.5 ARIMA模型预测

ARIMA经过定阶与模型检验之后，就可以进行数据训练并进行数据预测，然后再对预测的数据进行误差分析，判断选择的模型是否满足实际需要。

2 现场应用

2.1 工程概况

以石拉乌素矿地音监测数据为例，应用ARIMA模型进行数据预测。石拉乌素矿采区面积约26.6 km2，煤层底板标高分别为+691～+765 m、+660～+710 m、+647～+694 m、+632～+684 m。本文监测的数据来自于103工作面巷道，工作面标高+685.7～+696.6 m，平均+695.15 m，工作面范围内地面标高+1366.5～+1337.3 m，平均+1350.0 m。

地音监测系统安装的位置为103工作面，煤层上方直接顶为粉砂岩，厚度约为15.67 m，上方赋存有一层中砾岩，硬度较大，厚度约为11.26 m。老顶为细粒砂岩，厚度约为64.37 m，直接底为砂质泥岩，厚度约为12.89 m。机械化综掘施工工艺，顶板采用锚网锚索联合支付方式。根据103工作面顶底板冲击倾向性结果显示，该工作面具有强冲击倾向性，其顶、底板经鉴定具有弱冲击倾向性。但上覆岩层中细粒砾岩由于厚度达64.37 m，且岩体抗压强度较高，一旦发生垮断将释放大量能量，造成冲击事故。

在103工作面巷道安装了ARES-5/E地音监测系统，该系统主要有三部分组成：ARES-5/E地面中心站、N/TSA-5.28/E发射器和SP-5.28/E探头。探头实时监测被监测区域内的地音信号，并将地音信号转化为电压信号，再传送给N/TSA-5.28/E发射器，发射器把信号传递给地面中心站进行数据处理。

地音探头的安装一般在超前工作面约300 m、掘进迎头后方约150 m范围或其它重点区域进行实时监测。探头间距为50 m：对于回采工作面，地音探头与工作面最小距离为30 m，当工作面推至距离最近地音探头30 m时，将最近探头挪至最外侧探头以外50 m处。对于掘进工作面，地音探头与迎头最小距离为30 m，当迎头掘进至距离最近地音探头180 m时将最远探头挪至距离迎头30 m位置；在103工作面的上巷与下巷共布置四个SP-5.28/E地音探头，上巷为D3、D4，下巷为D1、D2，探头之间距离50 m，离切眼最近探头的距离在150 m范围内，103工作面地音探头布置如图1所示。

图1 103工作面地音探头布置

2.2 数据分析

2.2.1 数据获取

从系统安装并成功投入监测运行，每间隔1 h统计一次能量及相关频次，记录时间从2020年11月10日12时—2021年11月16日19时，经过数据处理与清洗共获取有效数据8913条。数据整体包含了工作班和检修班期间的数据，为分析该工作面整体地音活动运行连续数据记录。

在103工作面的监测区域，区域相对较小，综合四个地音传感器的监测数据，对四个通道数据在相同时间点求均值，去除多余的数据并对确定的数据进行插值补充，同时ARIMA模型要求参与运算的数列为平稳数据，观察单通道每小时的数据，不具备平稳性特征，故本文对四个通道数据求平均后又进行了降采样，降采样率为2D(2d)，并对使用降采样后的数据进行分析预测研究，数据曲线如图2所示。

图2 103工作面地音每日均值能量

从图2可以看出，从2020年12月15日—2021年1月15日，2021年04月03日—2021年05月10日，2021年08月初—2021年8月底时间段三个时间段103工作面地音能量集中释放，其余时间段相对平稳，周期大约为4个月，体现了103工作面煤岩动力地音活动规律。图2中，以2020年11月10日—2021年08月17日数据为训练集，以2021年08月19日—2021年09月06日共10条数据为测试集进行建模预测，并做预测误差分析，训练数据见表1。

表1 103工作面四通道地音能量训练数据

2.2.2 平稳性检验

ARIMA模型需要参与测试的数据为平稳性数列，平稳序列本身就存在特定的分布规律，前后数据具有一定的相关性和连续性，通过这些数列数据的相关性进而可以预测未来的数据；而非平稳数据凌乱无章，基本无规律可循，应用ARIMA模型就收效甚微。

103工作面地音数据相关性如图3所示，数列自相关系数下降非常缓慢，且结合图2发现曲线均值、方差波动较大，没有确定趋势，可以判断为非平稳序列。若要在地音序列数据中应用ARIMA模型，则必须对非平稳序列进行差分化操作，实现数据的平稳化。当数列自变量x变化到x+1时，函数y=f(x)的改变量为：

图3 103工作面地音数据相关性

Δyx=y(x+1)-y(x)，(x=0，1，2，…)

(10)

Δyx=yx+1-yx，(x=0，1，2，…)

(11)

经过差分化处理后，对数据再次进行平稳性检验，采用的方法为单位根和白噪声检验。

ARIMA测试数据ADF检验见表2，从表2中看出，t值同时小于1%、5%和10%临界值，P-value小于0.05，所以差分后的数据是属于平稳性数列。

表2 ARIMA测试数据ADF检验

对差分后的数列进行Ljung-Box检验，该检验的LB_p值为0.01141559，小于0.05，所以也证明该差分后的数列为平稳性数列，可以对该数列进行建模操作。

2.2.3 模型ARIMA参数p和q的确定与建模

在判定差分后的数列为平稳数列后，需要确定模型ARIMA(p，1，q)中参数p和q。

一种方法是根据数列的自相关图(ACF)和偏自相关图(PACF)的截尾与拖尾性来判断p值和q值。p值可从偏自相关系数(PACF)图的最大滞后点来大致判断，q值可从自相关系数(ACF)图的最大滞后点来大致判断，但这个定阶策略比较困难，有时不出现截尾，有时相关系数的曲线还存在小幅震荡等现象，这要确定相对精确的p，q值比较困难。

另一种方法是通过遍历搜索信息准则AIC和BIC最小参数组合来确定p、q值。本文数据通过对AIC和BIC测试并绘制热力图，在热力图中最低值的AIC或BIC组合中以最深的颜色显示，以这组颜色的值作为AR与MA系数，实现ARIMA模型的建模。本文数据的AIC的热力图如图4所示，BIC的热力图如图5所示。从图4中可以看到，AIC最小值为13452.03，在AR为3、MA为3处，所以ARIMA(p、d、q)选用的参数为(3，1，3)，因为已对数据进行了一次差分运算，所以d=1。从BIC热力图分析，ARIMA(p、d、q)选用的参数为(1，1，1)，现对AIC参数进行分析，同时比照BIC参数建模数据。

图4 103工作面地音数据AIC热力图

图5 103工作面地音数据BIC热力图

ARIMA模型应用参数(3，1，3)、(1，1，1)进行建模，模型参数输出见表3。由表3可以看出，ARIMA模型在AIC和BIC信息准则参数下的建模数据，在AIC参数表中，除const值外，其它的ar和ma系统的P值皆小于0.05，在参数误差范围之内。

表3 ARIMA模型参数计算

2.2.4 ARIMA模型参数检验

对已建模的ARIMA模型是否是能够满足生产需要的数据处理模型，或者说建模取得的参数是否能正确反应已知数据的潜在规律，需要给已建模的ARIMA模型进行参数检验[19]，通常的ARIMA模型参数检验主要针对残差进行正态性检验和自相关性检验，主要用的方法有通过自相关系数图、通过NormalQ-Q检验模型残差的正态特性，残差满足正态性，主要是因为残差集中于某一个数值，如果该值与0很接近，则它实际服从均值为0的正态分布，即它是一个白噪声。通过德宾-沃森(D-W)和Ljung-Box检验来检验模型残差的自相关性，如果残差满足非自相关性，则说明残差中不再包括AR或者MA过程产生的序列信号，有用的信号已被提取完毕。

ARIMA模型残差自相关系数如图6所示，可以看出，残差自相关系数置信区间在5%的置信水平下，所有的系数都是显著的，模型通过检验。残差数据来自正态分布，如图7所示。

图6 ARIMA模型残差自相关系数

图7 ARIMA模型Normal Q-Q

Ljung-Box检验数据见表4。检验的结果Prob(>Q)前12行的检验概率(一般观察滞后1～12阶)如果小于0.05、0.10等就拒绝原假设，就结果来看，前12阶的P值都是大于0.05，所以在0.05的显著性水平下，不拒绝原假设，即残差为白噪声序列。

通过ARIMA模型参数检验结果显示ARIMA模型残差符合正态分布和白噪声特性，所建模型参数符合ARIMA模型应用特性[20]。

2.2.5 ARIMA模型数据预测与结果误差分析

ARIMA模型预测数据曲线如图8所示，测试数据曲线延续了监测区域地音整体能量的后期走势。2021年8月19日—2021年09月06日的测试数据与预测数据对比见表5，对于ARIMA模型，在式(3)中μ参数进行统计已测试数据均值进行误差调整，通过μ调整后的表4数据计算误差为4.46%，在数据允许范围之内。用BIC生成的参数(1，1，1)参与计算，数据计算误差为5.82%。