基于V／S分析的瓦斯涌出量分形特性研究＊

乔美英 陈 鑫 兰建义

（1.河南理工大学电气工程与自动化学院，河南省焦作市，454000；2.河南理工大学能源科学与工程学院，河南省焦作市，454000）

1 概述

我国是煤与瓦斯突出灾害最严重的国家之一，由于突出是一种复杂非线性动力学行为，且瓦斯突出的形成条件、演化过程及诱发因素具有多样性、复杂性及随机性等特点，最终体现为确定性与随机性相统一的特点。基于此特点，当前煤矿工作面瓦斯涌出量预测方法均存在中长期预测方面的不足。

Hurst指数早已被用于许多领域做中长期趋势分析研究，目前，在大多数求Hurst指数的文献中，常用的方法是经典R／S法、Lo法及V／S法。而国内的文献在研究分形时，大多集中在利用经典R／S分析上，仅经济学领域对Lo法与V／S法有涉及。近年来Hurst指数在煤炭领域也做了相关的研究：黄存捍等将R／S分析用于矿井涌水量的预测，李业学等将其用于围岩变形趋势分析，李远耀等将其用于滑坡变形趋势预测。

本文在已有成果的研究基础之上，将改进的R／S分析法引入到掘进工作面的瓦斯涌出量的中长期趋势预测中。文中在详细介绍了R／S 分析基本算法原理基础上，给出了算法实现步骤并利用MATLAB2009b对其进行实现，同时利用两种噪声序列进行了验证。最后选用鹤壁十矿1113掘进工作面瓦斯涌出量时间序列进行分析研究，分别计算了该工作面3次发生突出时间序列的Hurst值及两次未突出的时间序列的Hurst指数值，对所研究的几个时间序列的Hurst值与长程相关的时间限度进行比较分析。

2 基本原理

R／S分析法是一种已经被广泛使用的非参数统计方法，其最大的优点是不必假定时间序列测度的分布特征，而且不论是正态分布还是非正态分布，R／S分析的结果都是稳健的。R／S分析法是通过计算非线性时间序列的Hurst指数H 值来判断时序信号的分形特性及长程相关性，以此区分时间序列的随机性与非随机性，最终确定趋势性持续及强度。

2.1 R／S分析法

经典R／S分析法基本原理：设长度为N 的时间序列数据 x（1），x（2），…，x（N

｛｝） ，将其分割为长度为n （2≤n≤L）的A 个相邻不重叠子区间Ma（a=1，2，…，A），其中L 为最大子区间长度，由上面的分析有A×n=N。记Ma中的每个元素为mk，a（k=1，2，…，n），则可得每个子区间的数学期望为：

式中：Xa——每个子区间的数学期望值；

mk，a——Ma中的每个元素。

计算每一个子序列的偏离均值差值序列：

式中：ΔXk，a——每个子序列的偏离均值差值。

ΔXk，a的均值为零，这一步为重标度或标准化。

由此可知，累积均值离差Yk，a、标准差SMa以及极差RMa分别为：

式中：Yk，a——累积均值离差；

SMa——标准差；

RMa——极差。

而经典R／S 分析法的统计量，计算每个子区间的重标度极差值RMa／SMa及A 个区间平均重标度极差：

Mandelbrot证明R／S（n）与n有幂律关系，即R／S ∝nH。所以lg（R／S（n））与lgn之间存在线性关系：

实际求取过程中是在双对数坐标中绘制出点（lgn，lg（R／S）），并进行回归分析，如果坐标点在双对数坐标中表现出很好的线性关系，则直线的斜率即为H 值，即Hurst指数。经典R／S分析法在不需要任何假设前提下所得结果最稳定，也是最接近时间序列的实际特性，特别是经典R／S 分析法能够从分形时间序列中区分出随机时间序列，并能够计算出分形时间序列内在的非周期循环长度，可以更深刻地揭示非线性时间序列内在统计规律。

2.2 V／S分析法

经典R／S 分析方法存在一定的局限性，若一个时间序列数据表现出很强的短期相关性，此时经典R／S分析方法会出现偏差，更易得出具有长期相关性结论。鉴于这方面的原因，Lo法在经典R／S分析方法的基础上提出了修正的R／S分析方法。Lo法通过引入样本协方差来对经典的R／S分析法进行修正。Lo法使用如下的统计量：

γj——j阶样本自协方差。在此方法中q的选择非常重要，其直接影响到检验效果。Lo法给出选择的最优q值如下：

式中：e——子样本长度；

＾ρ——一阶自相关系数的估计值。

当q=0时，Lo法成为经典R／S法。Lo法比经典法具有一个明显的优点就是剔除了序列短期相关性，能更有效地检测出时间序列的长期相关性。然而Vadim Teverovsky的研究表明：Lo法对某些具有显著长期相关性的时间序列更倾向于得出没有长期相关性的结果。所以，Lo法可能会低估序列的长期相关性。L.Giraitis（2003）在Lo法的基础上提出V／S法，其使用的统计量为：

2.3 Hurst指数判定长程相关特性分析

长程相关性 （长记忆性）通常被定义为时间序列数据的具有一种持久性和长期信赖关系，可用Hurst指数的相应函数来描述。间隔为τ的时间序列数据间的关联函数C（τ）为H 的函数，其表达式为：C（τ）=τ2H-1-1。当H=0.5时，C（τ）=0表明时间序列数据间不相关；当H ＞0.5时，C（τ）＞0表明时间序列数据间正相关；当H ＜0.5 时，C（τ）＜0表明时间序列数据间负相关。时间序列Hurst指数值H 与时间序列的长程相关性的关系具体可表述为：

（1）如果H 值接近0.5，则此时间序列可以用随机游走描述，如利用计算机随机生成白噪声序列数据，或是把尼罗河流量时间序列打乱，再进行R／S分析，得到的Hurst指数值也接近0.5。

（2）若0.5 ＜H ＜1，表明存在状态持续性即长期记忆性，变量间具有相关性，表现为趋势追随倾向，暗示着时间序列是一个持久性的或趋势增强的序列，例如若序列在前一段时期是向上 （下），那么在下一段时期将越有可能继续是正 （负），H 越接近于1，相关性越强，H越接近于0.5，其噪声越大，趋势也越不确定。从理论上来说，这种长期记忆性表示当前发生的事件对以后发生的事件会有影响，用混沌动力学的语言来表述这一现象就是系统对初始条件的敏感信赖性，需要说明的一点是这种长期记忆性并不会随着时间标度的改变而改变，例如若以时间增量为分钟间隔来观察，则未来的分钟变化与过去分钟变化相关，若改为以日来观察，则未来的日变化与过去的日变化相关。

（3）当0＜H ＜0.5时，存在逆状态持续性，时间序列是反持久性的，表现为均值回复倾向，此时序列比随机序列更具有突变性与波动性。若序列在前一个期间向上 （下）走，则它在下一个期间就越有可能向下 （上）走，反持久性强度随着H 接近于0逐步加强。

总之，只要H ≠0.5，就可以用有偏的布朗运动 （分形布朗运动）来描述该时间序列数据。

2.4 长程相关性的时间限度分析

绝大多数复杂非线性时间序列的长程相关性是有时间尺度做为界限的，若超过这一时间尺度将表现出不相关的随机特性，自相似行为丧失，而其统计分布也随之变为接近正态分布，因此时间序列的长程相关性是有时间限度的，而这个时间限度常用无标度区来衡量。本文中用临界点Tm（Crossover point）来表示这一时间尺度，在求几种R／S分析中绘制log（R／S（n））-log（n）回归图中，这一时间临界特性表现为拟合直线斜率的变化。可由各斜率出现明显转折点处的值可推算出长程相关的最大时间尺度Tm。这一时间尺度也就是时间序列非周期循环 （Non-periodic cycle）的平均循环长度 （Average cycle length）。

E.Peters引入了下面的统计量来判断平均循环周期：

当R／S（n）与n 的平方根同步增长时，Vn-log（n）的散点图将分布在一条水平线上，若比n的平方根增长较快时，即持续性存在，则出现上升趋势，如果比n的平方根增长较慢时，即存在反持续性，则出现下降趋势。但是在实际中Vn-log（n）可能有好几个拐点，这些拐点可以作为重要参考。在本文中将同时利用Vn-log（n）与log（R／S（n））-log（n）曲线图来确定较明显的转折斜率点。如果在某一处出现明显转折，说明此时出现临界时间尺度的对数值，再通过相应的换算即可求出有效的时间序列长程相关性的限度。

而且在Vn-log（n）图中有另一个显著的特点，若0.5＜H ≤1时，Vn-log（n）曲线图是一条向上倾斜的直线，若0＜H ＜0.5时，则Vn-log（n）是一条向下倾斜的直线。确定长程相关性的时间限度为利用时间序列趋势预测提供许多定量参考。

3 实例分析

3.1 数据来源

鹤壁煤业集团十矿为煤与瓦斯突出矿井，到目前为止十矿共发生5次突出，其中5次突出中有4次发生在1113工作面，所以本文以1113工作面做为研究对象。

本文利用每隔5min的工作面瓦斯涌出量时间序列来研究工作面瓦斯涌出的复杂行为特性，对于已发生突出的时间序列数据，本文选取包括突出发生当天与前一天的检测数据，即共48h时间序列数据，所以本文时间序列至少为576个。

文中给出的数据是已经对原始数据经过初步分析与预处理，利用前一监测值与后一监测值之和平均来代替。首先给出该工作未突出时的两次时间序列数据。未突出的时间序列数据均选取了3d 时间，数据长度为864个。未突出时瓦斯时间序列数据见图1所示。

图1 未突出时瓦斯时间序列数据

由图1可知，未发生煤与瓦斯突出时瓦斯浓度变化趋势均匀且平缓。掘进中无突出现象发生，瓦斯浓度变化平稳，变化幅度小。

1113工作面3 次已突出的瓦斯涌出量时间序列数据波形见图2所示。由图2可以看出，这3次煤与瓦斯突出前或是突出过程中均出现了瓦斯含量即涌出量变化异常。而这种异常现象用普通的统计方法很难深刻刻画。本文引入V／S分析法对1113掘进工作面的瓦斯涌出时间序列数据进行分析。

3.2 算法实现步骤

（1）导入数据；

（2）按照时间序列数据的长度分割子序列，子序列的最短长度为10，最长长度为总时间序列长度的1／2，在程序过程中子序列的长度由最短依次增加，共生成的子序列数为50；

（3）利用autocorr函数求出时间序列的自相关函数并求q值；

（4）设置子序列分割数进行循环计算，计算每一个子序列下R 值与S 值；

图2 发生突出时瓦斯时间序列数据

（5）根据式 （1），利用reshape函数将依次生成子序列矩阵，利用mean函数得均值。据式 （2）和 （3）求子序列的差值序列及每个差值序列的累积和；

（6）据式（4）求子序列标准差；

（7）据式（6）计算每个子序列的重标度极差，并求其平均后的对数值；

（8）利用式 （16）求Vn值；

（9）重复（4）～ （7）步后计算完所有子序列参数 时，用 最 小 二 乘 法 拟 合 求log（R／S（n））-log（n）的斜率值，并绘制log（R／S（n））-log（n）与Vn-log（n）曲线，求曲线的转折点。并求出明显的转折点处的横坐标值。将此值换算可得出所研究时间序列的长程相关性的限度。

3.3 实例分析

根据程序计算出不同的分割长度n 下的log10（VS）-log10（n）曲线，在此基础上利用最小二乘拟合其斜率即为H 值，如图3 （a）、图4（a）、图5 （a）、图6 （a）、图7 （a）所示。同时在同一个横坐标下，再绘制出Vn-log10（n）关系曲线，利用这一曲线可以较明显地确定所分析时间序列的长程相关性的转折点，如图3 （b）、图4 （b）、图5 （b）、图6 （b）、图7 （b）所示，其中几个图中的横、纵坐标无量纲。瓦斯时间序列的Hurst曲线特征情况见表1所示。

图3 未突出时数据I长程趋势图

表1 瓦斯时间序列的Hurst曲线特征

图3～图7及表1表明，1113工作面的几组瓦斯时间序列的Hurst指数值均大于0.5，说明这几个时间序列数据服从分形布朗运动且具有正相关长程趋势特性。从其Vn曲线的走向上也可以看出，几种情况的Vn曲线均向上倾斜，进一步说明其趋势具有可持续性，在所测定的时间周期内表现的特征将得以延续。

从未发生突出的两组数据 （第I组与第II组）的特征参数中可以看出，未发生突出时其Hurst指数值均较大，第I组的Hurst指数为0.9137，第II组的Hurst指数为0.8358，说明未发生瓦斯突出时的数据其正相关特征较强，即用前面的数据来预测后面数据时可预测的有效时间要长。

图4 未突出时数据II长程趋势图

图5 突出发生时数据III的长程趋势图

图6 突出发生时数据IV 的长程趋势图

从图3中看出，Vn曲线的第一个小幅变化的横坐标处的值为1.699，而在较明显发生转折处的横坐标值为2，通过计算其有效关联长度为102。从图4中可以得出，Vn曲线的第一个小幅变化的横坐标处的值为1.833，而发生较大幅度的转折时的横坐标值为2.134，同样其对应的有效关联长度为102.134。由于有效关联长度是通过Vn曲线发生明显转折处的横坐标值计算得出，所以计算结果可能稍有出入。

图7 突出发生时数据V 的长程趋势图

对已发生突出的时间序列数据 （即第III、IV与V 组）的特征参数中可以看出，第III组数据的Hurst指数为0.7199，从图5中可以看出，其Vn曲线的第一个小幅变化的横坐标处的值为1.462，第二个发生转折的坐标值为1.556，文中选取这一点作为计算其关联长度的值，之后一个大的转折点处的坐标值为1.949。第IV 组数据的Hurst指数为0.6556，从图6中可以看出，Vn曲线的第一个小幅变化的横坐标处的值为1.462，第二个发生转折的坐标值为1.556，文中选取这一点作为计算其关联长度的值，之后一个大的转折点处的坐标值为1.949。第V 组数据的Hurst指数为0.6301，从图7中可以看出，Vn曲线的第一个较明显转折处的横坐标值为1.415，之后的Vn曲线均出现较大幅度的波动，特别是在2.161点之后出现了较大幅度的反趋势，说明此时间序列的突变特性特别强，同样从第III组数据的时间序列曲线中也可以看出，此次突出时瓦斯涌出量总体是呈现增长的趋势，趋势性较明显，所以所求得的Hurst指数值较大；从第IV 组数据曲线中可以看出发生突出前瓦斯涌出量有两次小的异常出现，但是在突出发生前却表现的比较正常，所以其趋势性没有第III组强，但是通过Hurst指数的求取，发现其序列仍然具有正向趋势性；从第V 组数据曲线中可以看出此次突出时其瓦斯含量的超限值较多，而且突出发生了两次，且涌出的瓦斯量也是最大的，所以其非线性特性特别复杂、随机性比上面几种情况都要复杂，所以其Hurst指数值最接近0.5，但是由于其是实际的动力学行为，所以并非纯粹的随机行为，表现的趋势行为仍然为正趋势。

4 结论

利用V／S法分析后发现，无论是发生突出时的时间序列数据，还是未突出的时间序列数据均表现出一定长程正相关性，说明未来的时间序列数据与之前的时间序列具有较强的正相关性。所以利用V／S法对掘进工作面瓦斯涌出量时间序列的长程相关特性的研究结果可作为工作面瓦斯突出预报的一种辅助手段。

对于发生突出时的瓦斯涌出时间序列的动态分形特性的深入研究是建立在突出实例基础之上，而这又与煤矿安全生产将相违背，所以利用V／S 分析法深入分析突出发生时的分形特性为突出时工作面瓦斯涌出量时间序列的分形特性寻找普遍规律还需要做大量的研究工作。

［1］ Hua Su，Lin Yang.RS Analysis of China Securities Markets ［J］ .Tsinghua Sicence and Technology，2003 （5）

［2］ 黄存捍，冯涛等.基于分形和支持向量机矿井涌水量的预测 ［J］.煤炭学报，2010 （5）

［3］ 李业学，刘建锋.基于R／S 分析法与分形理论的围岩变形特征研究 ［J］ .四川大学学报 （工程科学版），2010 （3）

［4］ 李远耀，殷坤龙，程温鸣.R／S 分析在滑坡变形趋势预测中的应用 ［J］.岩土工程学报，2010 （8）

［5］ J.Guan，N.B.Liu，Y.Huang，et al.Fractal characteristic in frequency domain for target detection within sea clutter［J］.IET Radar Sonar &Navigation，2012 （5）

［6］ V.Teverovsky，M.S.Taqqu，W.Willinger.A critical look at Lo's modified R／S statistic［J］.Journal of statistical Planning and Inference，1999 （1）

［7］ 李昊，李杰等.矿井瓦斯涌出量预测的GM （1，1）模型研究［J］.中国煤炭，2012 （9）

［8］ 王文静 .金融时间序列的长记忆特性及预测研究［D］.天津大学，2009

［9］ 黄诒蓉，罗奕.基于经典R／S分析方法的H 指数估计有效性评价 ［J］.统计与信息论坛，2009 （8）

［10］ 乔美英，马小平等.工作面瓦斯涌出量时间序列分形特性研究 ［J］.煤炭科学技术，2012 （12）

［11］ 刘彦伟，潘辉等.鹤壁矿区煤与瓦斯突出规律及其控制因素分析 ［J］.煤矿安全，2007 （12）

［12］ 林振山.长期预报的相空间分量组合法 ［J］.气象学报，1993 （3）

［13］ 何利文，施式亮，刘影.回采工作面瓦斯涌出的复杂性及其度量 ［J］.煤炭学报，2008 （5）