廣西盤龍鉛鋅礦涌水量時間序列變化特征分析及ARIMA預(yù)測

吳衛(wèi)忠，鄧 忠，陳余道，盧丹美，陳 盟，夏 源，鄒志坤，陸仁騫

(1.桂林理工大學(xué) 巖溶地區(qū)水污染控制與用水安全保障協(xié)同創(chuàng)新中心，廣西 桂林 541006；2.廣西壯族自治區(qū)水文地質(zhì)工程地質(zhì)隊，廣西 柳州 545006；3.自然資源部南方石山地區(qū)礦山地質(zhì)環(huán)境修復(fù)工程技術(shù)創(chuàng)新中心，南寧 530031)

0 引 言

礦山涌水的發(fā)生和演化受多重因素影響, 具有典型的非線性動力特征, 能夠?qū)е峦话l(fā)性涌水事件[1-2]。礦山突水是礦井事故的重要誘因[3], 開展礦山涌水量時間序列分析及預(yù)測研究, 對防治礦山突水、 保障安全生產(chǎn)具有重要意義[4-6]。

近年來, 時間序列的確定性分析法和非線性分析法在多領(lǐng)域的推廣應(yīng)用, 為發(fā)展礦山涌水量預(yù)測方法提供了新的視角[4, 7]。除了統(tǒng)計學(xué)方法外, 重標(biāo)極差法(R/S分析)和自回歸移動平均模型(ARIMA)較為典型。前者能夠用來分析時間序列的長程依賴性(LRD), 刻畫出序列變化的主要控制因素, 能較好地解釋非線性參數(shù)特征[8-9]、 降水量序列特征[10-11]和徑流特征[12]等; 后者作為時間序列預(yù)測的一種模型, 能較準(zhǔn)確地開展短期趨勢預(yù)測[13-15]。這兩種方法在礦山涌水量序列分析方面均有應(yīng)用的案例, 比如山東鄆城煤礦山工作面涌水量預(yù)測等[14]、 青龍煤礦月最大涌水量預(yù)測[16]、 東歡坨煤礦礦井涌水量預(yù)測[17]、 趙家寨煤礦礦井涌水量預(yù)測[18]等。這些應(yīng)用多集中在我國北方地區(qū)且以煤礦為主, 取得了很好的預(yù)測效果。對于南方降水充沛、 地表水豐富且與地下水聯(lián)系密切的巖溶區(qū)有色金屬礦山, 吳松明[19]采用水均衡法、 大井法和水文地質(zhì)比擬法等3種方法, 預(yù)測出鋁土礦區(qū)的礦坑涌水量, 但應(yīng)用R/S分析和ARIMA模型分析涌水量時間序列的案例未見報道。為此, 對廣西大水量巖溶礦山——盤龍鉛鋅礦開展了涌水量時間序列變化特征分析和ARIMA預(yù)測研究。

盤龍鉛鋅礦位于廣西武宣縣黔江河畔。自2008年開采以來, 礦業(yè)安全一直面臨地表水和地下水的侵害影響。根據(jù)礦山規(guī)劃, 未來將向深部開采, 開采標(biāo)高-440～-1 100 m; 加上2019年大藤峽水利樞紐黔江河段蓄水, 礦坑涌水也將成為影響安全生產(chǎn)的重要因素之一。

本文根據(jù)2010—2021年盤龍鉛鋅礦礦坑涌水量監(jiān)測數(shù)據(jù), 結(jié)合礦區(qū)水文地質(zhì)條件, 在分析礦坑涌水量時間序列結(jié)構(gòu)特征基礎(chǔ)上, 采用R/S方法分析涌水量序列LRD特征, 并應(yīng)用ARIMA模型進(jìn)行涌水量分段預(yù)測, 以期為礦山水安全監(jiān)測和水災(zāi)害防御提供參考。

1 礦區(qū)水文地質(zhì)概況

盤龍鉛鋅礦所在區(qū)域?qū)賮啛釒貪駳夂騾^(qū), 年平均氣溫21.1 ℃, 多年平均降水量1 370 mm, 降水多集中在5～9月份, 占全年降水量的69%。黔江是礦區(qū)東側(cè)毗鄰的河流, 多年平均流量4 240 m3/s, 常年水位38～45 m。地貌類型主要為構(gòu)造侵蝕、 溶蝕峰叢谷地及殘丘洼地(平原)。

盤龍鉛鋅礦是廣西大瑤山西側(cè)的大型沉積巖型鉛鋅礦床, 大地構(gòu)造上位于桂中凹陷帶與大瑤山隆起的結(jié)合部位, 經(jīng)歷了多期構(gòu)造作用[20]。礦區(qū)出露的地層有第四系、 石炭系、 泥盆系、寒武系(詳見圖1),主要由白云巖、 石灰?guī)r、 泥灰?guī)r、硅質(zhì)巖等組成。其中, 上倫組上段(D1sl2)白云巖為主要賦礦層位,在礦區(qū)出露面積最大(厚度約970 m), 是構(gòu)成大嶺礦段礦坑頂?shù)装逯苯映渌暮畬? 上倫組下段(D1sl1)為泥質(zhì)灰?guī)r與灰?guī)r互層,富水性較弱。礦區(qū)發(fā)育多條斷層, 其中北東向逆斷層F1構(gòu)成了礦區(qū)西側(cè)邊界; 近南北向平移斷層F2則將礦區(qū)分為崩山礦段和大嶺礦段(圖1), 前者已在2009年停止開采, 后者2008年以來處于開采狀態(tài)。

礦區(qū)水文地質(zhì)條件復(fù)雜, 以溶蝕裂隙發(fā)育為主, 富水性由強到弱不等, 可分為4個地下水子系統(tǒng): 大嶺礦地下水子系統(tǒng)(Ⅱ-1)、 大坪嶺地下水子系統(tǒng)(Ⅱ-2)、 崩山礦地下水子系統(tǒng)(Ⅱ-3)與六沙地下水子系統(tǒng)(Ⅱ-4)。各子系統(tǒng)相互之間有一定的水力聯(lián)系。大嶺礦地下水子系統(tǒng)(Ⅱ-1), 南側(cè)主要由下泥盆統(tǒng)上倫組下段、 郁江組(D1y)、 那高嶺組(D1n)、 蓮花山組(D1l)泥灰?guī)r、淺變質(zhì)砂巖和泥頁巖等構(gòu)成,地層傾角約50°～70°;北側(cè)由下泥盆統(tǒng)二塘組(D1e)泥灰?guī)r構(gòu)成;東側(cè)毗鄰黔江河流;西側(cè)與F2斷層相接(圖1)。該子系統(tǒng)邊界條件可概括為:南側(cè)和北側(cè)為隔水邊界,東部為水頭邊界,西部為弱透水邊界[21]。

2 礦坑涌水來源分析

根據(jù)礦區(qū)地下水系統(tǒng)劃分(圖1)及開拓中段涌水量監(jiān)測, 大嶺礦段坑道涌水主要有垂向入滲和東、 西側(cè)補給3個方向的水源:

(1)大氣降水、 低洼地表水入滲補給: 礦區(qū)近5年(2017—2021年)年均降水量1 403.24 mm, 積水區(qū)主要位于重晶石采坑、 水塘、 積水坑和溝渠等, 多為常年積水, 加大了地表水滲入補給地下水的水量。

(2)東部黔江河水側(cè)向反補給: 大嶺礦段距離黔江河岸約500 m, 近河岸的地下水水位常年低于黔江(平水期水位約38 m), 礦坑疏干水位為-115.00 m, 地下水水力坡度達(dá)到34%。礦坑地下水水位隨著黔江水位升降而升降, 水力聯(lián)系較密切。黔江河床底部為裸露的碳酸鹽巖, 河水可以通過溶洞、 溶蝕裂隙帶補給礦坑疏干區(qū)。2019年大藤峽水利樞紐蓄水后黔江水位提高20～30 m, 礦山在東部建成了擋水帷幕(圖1), 帷幕底部高程-125 m, 削弱了黔江對礦區(qū)的側(cè)向補給。

(3)西部崩山礦地下水子系統(tǒng)的側(cè)向補給: 西部崩山礦地下水子系統(tǒng)與大嶺礦地下水子系統(tǒng)被弱透水F2斷層分隔。大嶺礦段疏干排水時, 崩山礦地下水子系統(tǒng)的地下水能穿越F2斷層滲流補給大嶺礦坑, 補給強度受裂隙空間大小、 連通性、 水頭高度等因素影響, 補給量有限[22]。

盤龍鉛鋅礦坑有-70、 -120、 -170、 -220、 -270、 -320、 -380和-440 m共計8個標(biāo)高的開采中段, 各中段出水量成為礦坑涌水量主要來源, 2010年1月—2021年1月期間的監(jiān)測結(jié)果如圖2所示。據(jù)礦區(qū)勘察, -120 m標(biāo)高以上中段巖層平均滲透系數(shù)K=4.85 m/d。溶蝕裂隙、 溶洞發(fā)育, 鉆孔遇溶洞率48.24%, 巷道遇溶洞率74.07%, 巖溶中等發(fā)育。-120 m標(biāo)高以下中段巖層的平均滲透系數(shù)為0.40 m/d, 鉆孔遇溶洞率9.52%, 巷道遇溶洞率18.52%, 巖溶弱發(fā)育。隨著開采中段的縱深開拓, 涌水點數(shù)量增加, 但涌水量主要來源于-120 m以上中段。開采初期(2014年之前),-70、-120和-170 m中段平均涌水量分別為550.5、93.1和126.2 m3/h,以-70 m中段涌水為主; 2014年之后, -70 m中段涌水量平均為433.3 m3/h(在2021年9月異常增大, 達(dá)775.1 m3/h), -120～-270 m中段涌水量較小, -320 m中段涌水量平均為317.4 m3/h(最大值為2020年7月的597 m3/h), 2021年增加了-380和-440 m涌水點, 平均涌水量分別為104.1和79.5 m3/h。因此, 總體上歷年涌水量以-70 m和-320 m中段為主。

圖2 大嶺礦段各中段2010—2021年月均涌水量監(jiān)測結(jié)果Fig.2 Monitoring results of water inflow in each middle section of Daling ore block from 2010 to 2021

3 涌水量時間序列結(jié)構(gòu)特征分析

根據(jù)盤龍鉛鋅礦2010—2021年的月均涌水量時間序列(圖3), 年內(nèi)涌水量隨時間上下波動, 與降水特征相似, 受黔江水位影響, 具有季節(jié)性, 其中6—8月份涌水量較高, 具有中間高兩側(cè)低的分布特征。

圖3 盤龍鉛鋅礦2010—2021年月均礦坑涌水量時間序列Fig.3 Time series of mine water inflow of Panlong lead-zinc deposit from 2010 to 2021

3.1 涌水量時間序列結(jié)構(gòu)參數(shù)分析

對涌水量時間序列進(jìn)行統(tǒng)計, 得到平均值、 標(biāo)準(zhǔn)差、 偏度、 峰度和增強迪基-福勒(ADF)檢驗值(表1)。2010—2021年涌水量多年平均值為792.7 m3/h, 最小年均值為707.3 m3/h, 最大年均值為982.7 m3/h。2019、 2020和2021年年均涌水量有了明顯提高, 分別達(dá)865.8、 942.7和982.7 m3/h, 這不僅與礦坑縱深掘進(jìn)涌水量增加有關(guān), 而且可能與黔江蓄水增加側(cè)向補給有關(guān)[21]。

開采初期, 2010和2012年標(biāo)準(zhǔn)差較大, 與淺部含水層富水性差異大有關(guān), 導(dǎo)致涌水量波動幅度大; 隨著縱深開拓, 深部富水性差異逐漸變小, 地下涌水量序列趨于穩(wěn)定; 2019年1月大藤峽水利樞紐蓄水后又出現(xiàn)了涌水量大幅度增大?？傮w上, 2010—2021年涌水量序列的偏度為0.842, 具有右偏分布趨勢及其長尾特征, 說明后續(xù)隨著礦區(qū)的開采出現(xiàn)涌水量超常值可能性相對比較大。由于2021年礦區(qū)涌水量異常增長, 也使得2010—2021年涌水量時間序列呈現(xiàn)尖峰態(tài)特征; ADF檢驗統(tǒng)計量t值和概率P值分別為0.133和0.723, 涌水量序列為非平穩(wěn)時間序列[23]。

3.2 涌水量時間序列的長程依賴性分析

3.2.1 R/S分析方法 設(shè)不同時間t1，t2，…，tn所對應(yīng)的礦區(qū)涌水量分別為x1，x2，…，xn，記為

時間序列。在τ=tn-t1時間段內(nèi), 其礦區(qū)涌水量平均值為

(1)

在tj(tj-t1)這一時期的礦區(qū)涌水量對于平均值的積偏差為

(2)

不同的τ對應(yīng)不同的x(t,τ)序列, 將每一個τ值對應(yīng)的x(t,τ)序列中的極差定義為區(qū)間, 記為R(τ)

R(τ)=maxx(t,τ)-minx(t,τ)。

(3)

同一個τ值下, 與R(τ)區(qū)間相對應(yīng)的礦區(qū)涌水量標(biāo)準(zhǔn)偏差為

(4)

英國科學(xué)家赫斯特(Hurst)等認(rèn)為, 時間序列時段τ內(nèi)的極差R(τ)和標(biāo)準(zhǔn)偏差S(τ)的比值與τ/2之間屬于一種冪律關(guān)系[24]。引入無量綱比值R(τ)/S(τ), 對R(τ)重新進(jìn)行標(biāo)度, 有

(5)

時間序列時段τ內(nèi), 極差和標(biāo)準(zhǔn)偏差的比值與τ/2的關(guān)系

(6)

式中:τ/2的指數(shù)H被稱為Hurst指數(shù)。由不同時間序列得出的R/S值建立雙對數(shù)坐標(biāo)系ln(τ/2)-ln(R/S), 并應(yīng)用最小二乘法進(jìn)行線性擬合, 得出的斜率即為Hurst指數(shù)(H)。Hurst指數(shù)是反映動力系統(tǒng)標(biāo)度不變性定律的重要指標(biāo), 可用來反映時間序列變化趨勢的持久性和非持久性, 即長程依賴性, 通常在0～1, 無量綱。若0

3.2.2 涌水量序列的長程依賴性 根據(jù)2010—2021年盤龍鉛鋅礦礦坑涌水量時間序列所作的R/S分析結(jié)果如表2所示。各年涌水量序列的Hurst指數(shù)范圍為0.780 1～0.994 9(平均0.918 6),表明礦區(qū)涌水量在年內(nèi)的變化趨勢十分相似, 持續(xù)性強,其原因與年內(nèi)降水量和黔江水位呈季節(jié)性變化的規(guī)律有關(guān)。然而,2020和2021年Hurst指數(shù)分別低至0.839 6和0.780 1,相對于之前的年份序列持續(xù)性下降明顯。結(jié)合表1中的序列峰度, 可以認(rèn)為: 2020和2021年這兩年涌水量異常升高、 呈現(xiàn)尖峰態(tài)特征是Hurst指數(shù)下降的主要原因。

表2 涌水量時間序列R/S分析

對2010—2021年際間涌水量序列開展R/S分析, 可得到Hurst指數(shù)為0.700 5(范圍0.5～1), 說明其具有較弱的長程依賴性, 未來序列變化趨勢可能會相對過去趨勢出現(xiàn)偏移。推測其原因, 可能與水利工程樞紐蓄水(2019年1月開始蓄水)導(dǎo)致地下水系統(tǒng)水動力場異常變化有密切關(guān)系。

4 涌水量序列預(yù)測

4.1 時間序列ARIMA模型建立

對于非平穩(wěn)的涌水量時間序列, 經(jīng)差分轉(zhuǎn)換為平穩(wěn)時間序列, 可利用ARIMA模型對該序列進(jìn)行預(yù)測, 其模型結(jié)構(gòu)為[13, 25]

yt=c+φ1Xt-1+φ2Xt-2+…+φpXt-p+et-

θ1et-1-θ2et-2-…-θqet-q,

(7)

式中:yt為t時序數(shù)據(jù);c為常數(shù)項;et為白噪聲序列;p、q為模型的階數(shù);φp、θq為自回歸和移動平均系數(shù)。

根據(jù)礦區(qū)坑道涌水點分布及2010—2021年涌水量組成(圖2), 按-120 m標(biāo)高以上段(-70和-120 m中段)和以下段(-170～-400 m中段), 利用2010—2019年月涌水量分別開展短期2020年和長期2021年各月涌水量預(yù)測, 通過對比預(yù)測值和實測值反映ARIMA模型的適用性。

首先對盤龍鉛鋅礦2010—2019年-120 m標(biāo)高上段涌水量序列進(jìn)行ADF檢驗, 統(tǒng)計量t值為-0.951(大于1%水平臨界值),概率P為0.303(大于10%概率),可判別為非平穩(wěn)序列;對數(shù)據(jù)進(jìn)行一階差分,再進(jìn)行ADF檢驗,t值為-7.209(小于10%水平臨界值),概率P為0.000(小于1%概率),表明涌水量序列在一階差分后變成了平穩(wěn)序列;進(jìn)一步對平穩(wěn)序列進(jìn)行自相關(guān)性(AC)和偏自相關(guān)性(PAC)檢驗, 表明可應(yīng)用ARIMA模型[26]進(jìn)行分析預(yù)測。

通過赤池信息準(zhǔn)則(AIC)與施瓦茨準(zhǔn)則(SC)確定ARIMA模型的階數(shù)p和q, 當(dāng)p和q均等于2時, ARIMA模型的AIC與SC數(shù)值均最小, 說明該模型與原始礦坑涌水量序列擬合效果最佳。由此得到ARIMA模型表達(dá)式

yt=-2.648+1.521yt-1-0.722yt-2+et+

1.371et-1-0.371et-2。

(8)

對ARIMA模型進(jìn)行誤差檢驗, 表明該ARIMA模型符合要求可進(jìn)行涌水量預(yù)測。同理, 對盤龍鉛鋅礦2011—2019年-120 m標(biāo)高下段涌水量序列建立ARIMA模型, 可得到模型表達(dá)式

yt=4.419+1.070yt-1-0.649yt-2+et+

0.867et-1-0.161et-2。

(9)

4.2 預(yù)測結(jié)果與討論

利用上述ARIMA模型對2020和2021年盤龍鉛鋅礦-120 m標(biāo)高上段、 下段涌水量序列進(jìn)行擬合預(yù)測, 結(jié)果見表3, 模型相關(guān)的平均絕對誤差(MAE)、 均方根誤差(RMSE)、 平均絕對百分比誤差(MAPE)與全年礦坑總涌水量誤差見表4。

表3 -120 m標(biāo)高上、下段涌水量分段預(yù)測結(jié)果

表4 -120 m標(biāo)高上、下段涌水量預(yù)測評價結(jié)果

2020年-120 m標(biāo)高上段和下段礦坑涌水量預(yù)測值分別為4 323.7和6 274.0 m3/h, 與實測值誤差分別為2.49%和3.14%; 2021年預(yù)測值分別為 000.4和6 837.0 m3/h, 與實測值誤差分別為26.38%和9.64%。

2020年分段預(yù)測的MAPE均小于10%, 說明ARIMA模型預(yù)測準(zhǔn)確, 并且-120 m標(biāo)高上段預(yù)測值的MAE、RMSE與全年礦坑涌水量總誤差表現(xiàn)優(yōu)于-120 m標(biāo)高下段, 均能夠滿足實際工程短期預(yù)報的需求。模型具有可靠性, 可用于盤龍鉛鋅礦礦坑突水預(yù)防工作。對于2021年涌水量預(yù)測結(jié)果, 預(yù)測時間較遠(yuǎn), 且礦區(qū)開拓層不斷加深, 其精度比2020年結(jié)果差, 主要原因是模型模擬未能考慮2020和2021年的年度異常涌水。因此, 用ARIMA模型作短期預(yù)報更為合適。

5 結(jié) 論

(1)2010—2021年盤龍鉛鋅礦礦坑涌水量年內(nèi)變化受當(dāng)?shù)丶竟?jié)性降雨量變化的影響, -70和-320 m標(biāo)高中段是主要的涌水段; 涌水量時間序列總體上呈現(xiàn)右偏長尾、 尖峰-平峰態(tài)和非平穩(wěn)特征。

(2)年內(nèi)涌水量序列具有明顯的長程依賴性, 未來發(fā)展趨勢與歷年趨勢之間呈現(xiàn)出較強的持續(xù)性, 與氣象水文要素的季節(jié)性變化有關(guān)。年際涌水量序列(2010—2021年)具有較弱的長程依賴性, 過去的趨勢存在偏移, 可能與礦山開采深度加大、 水文地質(zhì)條件變化和黔江蓄水活動有關(guān)。

(3)ARIMA模型可用于盤龍鉛鋅礦礦坑涌水量的短期預(yù)測, 按-120 m標(biāo)高上下分段預(yù)測較為合理, 能滿足礦坑突水防御的需要。該模型預(yù)測需及時更新監(jiān)測數(shù)據(jù), 校正模型, 以提高模型精度。