基于R/S分析的礦井涌水量分形特性研究

喬美英++劉震震++程鵬飛

摘要：為研究受多重因素影響的礦井涌水量特性，利用重標(biāo)極差（R/S）分析法對(duì)千秋煤礦礦井涌水量時(shí)間序列進(jìn)行了分形特性研究。求得礦區(qū)及各采區(qū)涌水量Hurst指數(shù)為0.89～0.98，分形維數(shù)1.02～1.12，對(duì)應(yīng)時(shí)間序列的有效關(guān)聯(lián)長(zhǎng)度104.61～105.35。結(jié)果表明：千秋煤礦礦井涌水量序列具有非常顯著的分形特性，且涌水量上升趨勢(shì)明顯，礦區(qū)總涌水量和采區(qū)涌水量的內(nèi)在聯(lián)系在分形特性上也有一定體現(xiàn)，為研究受多重因素影響的礦井涌水量時(shí)間序列提供了新思路，同時(shí)也為礦井涌水的研究與治理增加了理論實(shí)證。

關(guān)鍵詞：R/S分析；礦井涌水量；分形特性；Hurst指數(shù)；時(shí)間序列

DOIDOI：10.11907/rjdk.171686

中圖分類號(hào)：TP319文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：16727800（2017）010014204

0引言

礦井水災(zāi)被稱為煤礦五大自然災(zāi)害之一，危害程度一定程度上取決于礦井涌水量大小[1]。礦井涌水量預(yù)測(cè)受多種不確定因素影響，伴隨著開(kāi)采深度增加和開(kāi)采面積擴(kuò)大，涌水量預(yù)測(cè)變得越來(lái)越復(fù)雜。但涌水量的變化并非無(wú)序和無(wú)規(guī)律性，其有著一定的分形特性[23]。分形特性是局部與整體具有某種相似性的特性。在分形特性研究中，R/S分析法是一種運(yùn)用十分廣泛的分析方法。該方法由英國(guó)科學(xué)家H·E·Hurst提出并用于水文方面研究，近年來(lái)在水文、氣象、股市等諸多領(lǐng)域的時(shí)間序列預(yù)測(cè)中皆有運(yùn)用，且效果很好[47]。R/S分析法在礦井涌水量預(yù)測(cè)方面也有較多成果，如將灰色理論與R/S分析相結(jié)合對(duì)礦井涌水量進(jìn)行預(yù)測(cè)[810]、基于分形和支持向量機(jī)的涌水量預(yù)測(cè)[11]、利用域重新標(biāo)度法對(duì)涌水量時(shí)間序列預(yù)測(cè)[12]等。

礦井涌水量時(shí)間序列是依據(jù)時(shí)間序列排列的礦井涌水集合。受到多種因素不同作用機(jī)理影響，涌水量變化相對(duì)復(fù)雜。通過(guò)對(duì)時(shí)間序列進(jìn)行分析，可以找出影響系統(tǒng)的特征量，將復(fù)雜系統(tǒng)研究轉(zhuǎn)換為簡(jiǎn)明的特征量研究[13]。本文將時(shí)間序列的R/S分析方法引用到礦井涌水量的分形特性分析，以河南省千秋煤礦礦井涌水量為例，分別求取兩組工作面和總礦井涌水量序列的Hurst值與相應(yīng)的有效關(guān)聯(lián)長(zhǎng)度，并進(jìn)行比較分析。此方法有助于防治礦井水害，對(duì)礦井涌水量變化規(guī)律的預(yù)測(cè)與研究具有指導(dǎo)意義。

1R/S分析基本原理

R/S分析法是分形特性研究中廣泛應(yīng)用的方法[14]，通過(guò)時(shí)間序列的尺度改變，探究在不同尺度范圍的統(tǒng)計(jì)規(guī)律，從而在不同尺度之間進(jìn)行轉(zhuǎn)換[15]。根據(jù)隨機(jī)時(shí)間序列Hurst指數(shù)的計(jì)算結(jié)果對(duì)時(shí)間序列加以區(qū)分，并進(jìn)行時(shí)間序列信號(hào)的分形特性分析和相關(guān)性判定，由此確定時(shí)間序列的過(guò)去與未來(lái)變化趨勢(shì)。通過(guò)Hurst指數(shù)可以判定礦井涌水量時(shí)間序列的分形結(jié)構(gòu)及狀態(tài)趨勢(shì)；礦井涌水量系統(tǒng)的記憶時(shí)間長(zhǎng)度也可由平均循環(huán)長(zhǎng)度計(jì)算得到。

R/S分析法求解[16]過(guò)程：

設(shè)定長(zhǎng)度為N的涌水量時(shí)間序列{x（k）}，將其等分為長(zhǎng)度為n的M個(gè)連續(xù)子序列，即Mn=N（n選取為不小于3的整數(shù)），其中，各連續(xù)子序列元素記為Dm（m=1，2，…M），對(duì)應(yīng)元素為xk，m，求取Dm的均值m。

構(gòu)建n個(gè)區(qū)間的Dm累積離差：

Xk，m=∑Nk=1（xk，m-m）（1）

計(jì)算式（1）中累積離差最大和最小值的差值Rm：

Rm=max1≤k≤N{Xk，m}-min1≤k≤N{Xk，m}（2）

標(biāo)準(zhǔn)差序列為：

Sm=（1K∑Nk=1（x（i）-（M））2）12（3）

根據(jù)式（2）和式（3）結(jié)果計(jì)算Dm的重標(biāo)極差：

（R/S）m=Rm/Sm（4）

對(duì)每個(gè)子序列進(jìn)行公式（1）～（4）的重標(biāo)極差運(yùn)算，得到對(duì)應(yīng)序列的重標(biāo)極差為：

（R/S）=1M∑Mm=1（R/S）m（5）

不斷增加n值，直至n=N/3，由此得出{x（k）}的重標(biāo)極差系列，再對(duì)其進(jìn)行取對(duì)數(shù)處理：

lg（R/S）n=lgC+Hlg n（6）

式（6）中，以lg n為自變量、等式左邊為因變量作散點(diǎn)圖，采用ploy fit擬合函數(shù)進(jìn)行擬合，得到直線的斜率即為Hurst值。

Hurst指數(shù)H不僅是R/S分析法中重要的統(tǒng)計(jì)量，也是一個(gè)現(xiàn)象表征量。研究發(fā)現(xiàn)：若Hurst值H接近0.5，則涌水量時(shí)間序列為完全不確定性的隨機(jī)序列；若H為0～1之間且不等于1的數(shù)值，則涌水量時(shí)間序列未來(lái)與過(guò)去具有相關(guān)的變化趨勢(shì)；當(dāng)0≤H<0.5時(shí)，存在狀態(tài)逆持續(xù)性，為一種反持久性序列，即涌水量序列在之前的時(shí)間段呈現(xiàn)向上（下）浮動(dòng)，序列在下一時(shí)間段會(huì)有向下（上）浮動(dòng)的趨勢(shì)，趨勢(shì)強(qiáng)度隨著H接近于0逐漸增加。若0.5

2算法實(shí)現(xiàn)步驟及程序驗(yàn)證

2.1Hurst指數(shù)求取算法步驟

實(shí)驗(yàn)準(zhǔn)備：MATLAB2014b軟件，實(shí)驗(yàn)程序，涌水量數(shù)據(jù)。

（1）將時(shí)間序列的涌水量數(shù)據(jù)導(dǎo)入軟件空間，并將序列作為子序列；計(jì)算數(shù)據(jù)平均值以及其方差序列。求方差的累加和向量程序采用cumsum函數(shù)，同時(shí)求出子序列中極差與標(biāo)準(zhǔn)差之比的平均值，得出平均重標(biāo)度極差值Rn/Sn。

（2）將涌水量序列分為2，3，4，…，n分別計(jì)算步驟（1），得出其重標(biāo)度極差值Rn/Sn，直到子序列的維數(shù)最后小于或等于給定的閾值，循環(huán)結(jié)束。

（3）對(duì)各子序列數(shù)n與對(duì)應(yīng)得到的Rn/Sn分別取對(duì)數(shù)，再用最小二乘擬合函數(shù)ploy fit對(duì)這些雙對(duì)數(shù)坐標(biāo)點(diǎn)進(jìn)行擬合，通過(guò)計(jì)算擬合曲線的斜率值得到Hurst值。

（4）繪制擬合曲線，利用plot函數(shù)，以各子序列數(shù)n的對(duì)數(shù)為自變量，對(duì)應(yīng)的重標(biāo)度極差值Rn/Sn的對(duì)數(shù)為因變量，繪制出擬合曲線，曲線對(duì)應(yīng)的斜率為Hurst值，其中轉(zhuǎn)折點(diǎn)最大處所對(duì)應(yīng)的橫坐標(biāo)值即為突變點(diǎn)橫坐標(biāo)值。將此數(shù)值進(jìn)一步計(jì)算可得到有效關(guān)聯(lián)長(zhǎng)度，據(jù)此開(kāi)展實(shí)驗(yàn)數(shù)據(jù)的時(shí)間序列相關(guān)性研究。

2.2Hurst指數(shù)求取程序

白噪聲與有色噪音時(shí)間序列的Hurst數(shù)值分別對(duì)應(yīng)于接近0.5和略小于0.5，利用白噪聲與有色噪聲程序進(jìn)行驗(yàn)證。根據(jù)千秋煤礦礦井涌水量時(shí)間序列的長(zhǎng)短，采取相對(duì)應(yīng)長(zhǎng)短的序列進(jìn)行驗(yàn)證。由仿真軟件隨機(jī)生成序列長(zhǎng)度為60、120、240、1 200的有色噪聲，并分別將每組隨機(jī)數(shù)據(jù)在程序上運(yùn)行20次。求取有色噪音序列對(duì)應(yīng)Hurst值小于0.5的頻次，并以數(shù)據(jù)長(zhǎng)度60為例，繪出Hurst指數(shù)波動(dòng)圖，如圖1所示。從圖中可以看到，H值主要分布于0.4～0.6之間。同樣利用白噪聲序列數(shù)據(jù)進(jìn)行驗(yàn)證，隨機(jī)生成序列長(zhǎng)度60、120、240、1 200個(gè)白噪聲，求取白噪音序列對(duì)應(yīng)的Hurst均值。以數(shù)據(jù)長(zhǎng)度60為例，Hurst指數(shù)波動(dòng)如圖2所示。

白噪聲與有色噪聲驗(yàn)證結(jié)果匯總于表1。從表1可以看出本文所編程序的準(zhǔn)確性，且具有相對(duì)較強(qiáng)的穩(wěn)健性。

3實(shí)例研究及結(jié)果分析

3.1千秋礦區(qū)水文地質(zhì)條件

千秋煤礦井田南澗河旱季流量為2.26m3/s，雨季山洪暴發(fā)流量為1 446.5m3/s。煤礦附近的石河源于煤礦北部山區(qū)，河床水流受季節(jié)影響，平時(shí)水流量少，雨后水流增大，持續(xù)時(shí)間短，為典型的季節(jié)性河流。千秋煤礦礦山地質(zhì)構(gòu)造并不復(fù)雜，地表覆蓋著第四系粘土黃土，含水層之間水層穩(wěn)定，地表徑流中的大氣降水和地表水水量多由地表河流排出，向煤礦井下滲透較少。其中砂巖、礫巖含水層是礦井涌水主要因素。

3.2礦井涌水量數(shù)據(jù)分析

為了精準(zhǔn)反映礦井總涌水量與各采區(qū)情況的內(nèi)在聯(lián)系，通過(guò)收集整理一組礦區(qū)總涌水量數(shù)據(jù)和兩組采區(qū)工作面數(shù)據(jù)（數(shù)據(jù)均采用2011-2015年共60個(gè)觀測(cè)數(shù)據(jù)）進(jìn)行實(shí)例分析，見(jiàn)圖3。

從圖中可以看出，一水平采區(qū)涌水量變化較快，二一采區(qū)相對(duì)走勢(shì)平坦；礦井總涌水量峰值可達(dá)537.40m3/h，平均值375.68m3/h；一水平采區(qū)涌水量平均值157.13m3/h，二一采區(qū)涌水量平均值136.87m3/h。整體來(lái)看涌水量大小波動(dòng)明顯，蘊(yùn)含一定的分形特性。這種波動(dòng)看似毫無(wú)規(guī)律，但除卻隨機(jī)因素，包含地下水演變等在內(nèi)的規(guī)律性因素影響卻隱藏其中。

4討論

礦井涌水量擬合曲線如圖6所示，其中橫縱坐標(biāo)分別代表不等長(zhǎng)度n的對(duì)數(shù)值和重標(biāo)度極差Rn/Sn的對(duì)數(shù)大小。擬合曲線的率即Hurst指數(shù)值。

最后得到礦區(qū)Hurst指數(shù)為0.98，一水平采區(qū)Hurst指數(shù)0.89，二一采區(qū)Hurst指數(shù)0.97，H值均接近1，表現(xiàn)出非常明顯的持續(xù)效應(yīng)。

礦區(qū)和兩個(gè)采區(qū)R/S分析結(jié)果匯總見(jiàn)表2。

由表2可知，3個(gè)時(shí)間序列數(shù)據(jù)服從分形布朗運(yùn)動(dòng)，且時(shí)間序列具有長(zhǎng)期的正相關(guān)特征。千秋礦區(qū)礦井涌水量變化既非純粹的隨機(jī)性，亦非完全確定，而是隨機(jī)性和確定性相摻雜的復(fù)雜變化。

二一采區(qū)時(shí)間序列的Hurst指數(shù)值最大，其正趨勢(shì)性最強(qiáng)，而且突變點(diǎn)的橫坐標(biāo)值也比一水平采區(qū)大，表明二一采區(qū)的有效關(guān)聯(lián)長(zhǎng)度相對(duì)更長(zhǎng)，即三組實(shí)驗(yàn)時(shí)間序中最長(zhǎng)。三組時(shí)間序列的Hurst指數(shù)H分布在0.89～0.98區(qū)間內(nèi)，大小均接近于1，說(shuō)明涌水量是趨勢(shì)增強(qiáng)的序列，上升趨勢(shì)非常明顯；二一采區(qū)和礦區(qū)涌水量的Hurst指數(shù)H分別為0.97和0.98，說(shuō)明其涌水上升趨勢(shì)尤其明顯。

三組時(shí)間序列的分形維數(shù)都大于1，其值分布在1.02～1.12之間，表明礦區(qū)涌水量變化的復(fù)雜性，特別是一水平采區(qū)的復(fù)雜性明顯要強(qiáng)于礦區(qū)總涌水量和二一采區(qū)涌水量，這是因?yàn)橐凰讲蓞^(qū)涌水時(shí)間序列的分形維數(shù)最大，表明其防控系統(tǒng)更復(fù)雜。圖4中一水平采區(qū)涌水呈現(xiàn)出較大波動(dòng)也印證了這點(diǎn)。因此，在今后的礦區(qū)涌水量檢測(cè)與防治中要給予更多的關(guān)注。

從表1可以看到，三組時(shí)間序列的突出點(diǎn)橫坐標(biāo)較接近，說(shuō)明礦區(qū)和采區(qū)間涌水量有一定的內(nèi)在聯(lián)系。礦井、二一采區(qū)及一水平采區(qū)有效關(guān)聯(lián)長(zhǎng)度105.35>104.61>104.83，在這種趨勢(shì)下，可預(yù)測(cè)時(shí)間也是礦井涌水時(shí)間序列大于二一采區(qū)和一水平采區(qū)。

5結(jié)語(yǔ)

本文通過(guò)R/S分析法對(duì)河南千秋煤礦涌水量時(shí)間序列進(jìn)行分形分析，得到千秋煤礦涌水量時(shí)間序列Hurst數(shù)值大小接近且略小于1，說(shuō)明三組涌水量序列具有很強(qiáng)的趨勢(shì)性和持續(xù)性，表明該研究方法可對(duì)礦井涌水量進(jìn)行預(yù)測(cè)。通過(guò)對(duì)不同采區(qū)及總礦井涌水量的時(shí)間序列R/S分析，以及相關(guān)性分析，獲得了較好的驗(yàn)證，這為研究受多重因素影響的礦井涌水量時(shí)間序列提供了新的思路，也為礦井涌水研究與治理增加了理論實(shí)證。

參考文獻(xiàn)參考文獻(xiàn)：

[1]田洪勝.淺析礦井涌水量預(yù)測(cè)中的常用方法[J].地下水，2013，35（3）：2527.

[2]武強(qiáng)，趙蘇啟，董書寧，等.煤礦防治水手冊(cè)[M].北京：煤炭工業(yè)出版社，2013.

[3]陳玉華，楊永國(guó)，彭高輝.礦井涌水量混沌時(shí)間序列分析與預(yù)測(cè)[J].煤田地質(zhì)勘探，2008，36（4）：3436.

[4]李潤(rùn)求，施式亮，伍愛(ài)友，等.煤礦瓦斯災(zāi)害事故的分形特征[J].中國(guó)安全生產(chǎn)科學(xué)技術(shù)，2014，10（9）：2529.

[5]施式亮，李潤(rùn)求，何利文，等.基于分形學(xué)的瓦斯爆炸事故時(shí)序數(shù)據(jù)分析模型及應(yīng)用[J].中國(guó)安全科學(xué)學(xué)報(bào)，2011，21（10）：1015.

[6]施式亮，李潤(rùn)求，念其峰.煤礦安全狀況關(guān)鍵指標(biāo)變化特征的R/S分析[J].中國(guó)安全科學(xué)學(xué)報(bào)，2012，22（9）：7984.

[7]喬美英，馬小平，喬鐵柱，等.工作面瓦斯涌出量時(shí)間序列分行分形特性研究[J].煤炭科學(xué)技術(shù)，2011，39（12）：5155.

[8]李寶玲，李建華，昝明君，等.河流年徑流量的R/S灰色預(yù)測(cè)[J].水文，2015，35（2）：4448.

[9]王皓.基于灰色模型的煤礦涌水量預(yù)測(cè)研究[J].中國(guó)煤炭地質(zhì)，2014，26（3）：3942.

[10]李建林，李志強(qiáng)，王心義，等.基于R/S分析的礦井涌水量灰色預(yù)測(cè)[J].安全與環(huán)境學(xué)報(bào)，2015，15（5）：610.

[11]黃存捍，馮濤，王衛(wèi)軍，等.基于分形和支持向量機(jī)礦井涌水量的預(yù)測(cè)[J].煤炭學(xué)報(bào)，2010，35（5）：806810.

[12]唐依民，肖江.礦井涌水量時(shí)間序列分析的域重新標(biāo)度法[J].煤炭學(xué)報(bào)，2000，25（1）：6670.

[13]張殷欽，胡偉，劉俊民.基于R/S分析法的地下水位動(dòng)態(tài)化趨勢(shì)分析[J].中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，43（12）：49124916.

[14]張濟(jì)忠.分形[M].北京：清華大學(xué)出版社，2011.

[15]ALVAREZRAMIREZ J，ECHEVERRIA J C，RODRIGUEZ E.Performance of a highdimensional R/S method for Hurst exponent estimation[J].Physica A：Statistical Mechanics and Its Applications，2008，387（26）：64526462.

[16]CHEN Y，SUN R，ZHOU A.An Improved hurst parameter estimator based on fractional fourier transform[J].Telecommunication Systems，2010，43（3/4）：197206.

責(zé)任編輯（責(zé)任編輯：杜能鋼）endprint