基于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）回歸圖中，這一時間臨界特性表現為擬合直線斜率的變化?？捎筛餍甭食霈F明顯轉折點處的值可推算出長程相關的最大時間尺度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）