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

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

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

0 引 言

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

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

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

本文根據2010—2021年盤龍鉛鋅礦礦坑涌水量監測數據, 結合礦區水文地質條件, 在分析礦坑涌水量時間序列結構特征基礎上, 采用R/S方法分析涌水量序列LRD特征, 并應用ARIMA模型進行涌水量分段預測, 以期為礦山水安全監測和水災害防御提供參考。

1 礦區水文地質概況

盤龍鉛鋅礦所在區域屬亞熱帶溫濕氣候區, 年平均氣溫21.1 ℃, 多年平均降水量1 370 mm, 降水多集中在5～9月份, 占全年降水量的69%。黔江是礦區東側毗鄰的河流, 多年平均流量4 240 m3/s, 常年水位38～45 m。地貌類型主要為構造侵蝕、 溶蝕峰叢谷地及殘丘洼地(平原)。

盤龍鉛鋅礦是廣西大瑤山西側的大型沉積巖型鉛鋅礦床, 大地構造上位于桂中凹陷帶與大瑤山隆起的結合部位, 經歷了多期構造作用[20]。礦區出露的地層有第四系、 石炭系、 泥盆系、寒武系(詳見圖1),主要由白云巖、 石灰巖、 泥灰巖、硅質巖等組成。其中, 上倫組上段(D1sl2)白云巖為主要賦礦層位,在礦區出露面積最大(厚度約970 m), 是構成大嶺礦段礦坑頂底板直接充水的含水層; 上倫組下段(D1sl1)為泥質灰巖與灰巖互層,富水性較弱。礦區發育多條斷層, 其中北東向逆斷層F1構成了礦區西側邊界; 近南北向平移斷層F2則將礦區分為崩山礦段和大嶺礦段(圖1), 前者已在2009年停止開采, 后者2008年以來處于開采狀態。

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

2 礦坑涌水來源分析

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

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

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

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

盤龍鉛鋅礦坑有-70、 -120、 -170、 -220、 -270、 -320、 -380和-440 m共計8個標高的開采中段, 各中段出水量成為礦坑涌水量主要來源, 2010年1月—2021年1月期間的監測結果如圖2所示。據礦區勘察, -120 m標高以上中段巖層平均滲透系數K=4.85 m/d。溶蝕裂隙、 溶洞發育, 鉆孔遇溶洞率48.24%, 巷道遇溶洞率74.07%, 巖溶中等發育。-120 m標高以下中段巖層的平均滲透系數為0.40 m/d, 鉆孔遇溶洞率9.52%, 巷道遇溶洞率18.52%, 巖溶弱發育。隨著開采中段的縱深開拓, 涌水點數量增加, 但涌水量主要來源于-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月異常增大, 達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年月均涌水量監測結果Fig.2 Monitoring results of water inflow in each middle section of Daling ore block from 2010 to 2021

3 涌水量時間序列結構特征分析

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

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

3.1 涌水量時間序列結構參數分析

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

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

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

3.2.1 R/S分析方法 設不同時間t1，t2，…，tn所對應的礦區涌水量分別為x1，x2，…，xn，記為

時間序列。在τ=tn-t1時間段內, 其礦區涌水量平均值為

(1)

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

(2)

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

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

(3)

同一個τ值下, 與R(τ)區間相對應的礦區涌水量標準偏差為

(4)

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

(5)

時間序列時段τ內, 極差和標準偏差的比值與τ/2的關系

(6)

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

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

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

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

4 涌水量序列預測

4.1 時間序列ARIMA模型建立

對于非平穩的涌水量時間序列, 經差分轉換為平穩時間序列, 可利用ARIMA模型對該序列進行預測, 其模型結構為[13, 25]

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

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

(7)

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

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

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

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

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

1.371et-1-0.371et-2。

(8)

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

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

0.867et-1-0.161et-2。

(9)

4.2 預測結果與討論

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

表3 -120 m標高上、下段涌水量分段預測結果

表4 -120 m標高上、下段涌水量預測評價結果

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

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

5 結 論

(1)2010—2021年盤龍鉛鋅礦礦坑涌水量年內變化受當地季節性降雨量變化的影響, -70和-320 m標高中段是主要的涌水段; 涌水量時間序列總體上呈現右偏長尾、 尖峰-平峰態和非平穩特征。

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

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