弓長嶺露天礦采空區頂板位移超前預測算法

李相熙，朱萬成，任 敏

(東北大學 資源與土木工程學院，遼寧 沈陽 110819)

露天礦邊坡下層采空區不可避免的會因為爆破、機械鑿巖、地表行車、地下水侵蝕等外部擾動發生頂板崩落、整體失穩，最終影響露天礦山的生產作業安全[1-2]。為了保證露天礦山人員和設備安全，提高生產效率和降低開采難度，有必要對露天礦山邊坡下層采空區狀態進行超前預測。

隨著采空區穩定性監測技術的不斷提高，已基本實現對采空區位移場、應力場等的自動化實時監測[3-4]。由于受靜電、氣溫變化和機械振動等因素的干擾，導致采集的數據攜帶大量的誤差數據，給后期數據的分析和應用帶來巨大困擾。常見的時間序列數據超前預測方法對數據質量均有嚴格要求，直接使用時間序列原始數據使預測結果失真，進而無法提供長期、穩定、可信的預測結果。這就需要對所采集數據中的誤差數據進行科學剔除，以篩選出有效的數據集。目前剔除時間序列數據中誤差數據的方法主要基于數理統計[5]、距離計算、密度分析、聚類分析或小波分析等[6-8]中的一種算法來實現。其中，數理統計方法適用于對隨機誤差數據的剔除；基于距離計算、密度分析和聚類分析的方法是對“離群”較大的雜亂數據進行剔除。但是，露天礦山長期監測到的采空區頂板位移數據中隨機分布著大量雜亂信號和隨機誤差，單一的誤差數據識別剔除方法不足以對其充分剔除。

在數據超前預測算法方面，目前對時間序列數據長期超前預測技術的研究還處于初步探索階段。徐峰[9]等提出將三峽庫區八字門滑坡體的位移分解為趨勢項位移和周期項位移，分別采用灰色模型及AR自回歸模型進行位移預測；LIU Zhongqiang[10]等通過建立長短期記憶法、隨機森林法等多種神經網絡模型，實現對三峽庫區“階躍”型邊坡位移的預測，預測準確度較高；趙靜波[11]等利用灰色理論分階段預測了露天礦邊坡變形指標，結果表明預測值與實測值發展趨勢一致；陶慧[12]等基于遺傳算法對多個地壓監測變量進行重構，通過廣義回歸神經網絡( GRNN )實現對沖擊地壓的預測；蔡舒凌[13]等通過螢火蟲模型( FA模型 )與非線性自回歸模型( NAR模型 )的結合，提高單一自回歸模型預測的準確度，實現了對隧道圍巖變形的預測；曹文貴[14]等提出動態改變時間序列步長，再進行灰色預測的方法，提高了基樁峰值載荷預測的準確度；李恒凱[15]等建立由多種修改的灰色模型組成的預測模型庫，應用灰色評價模型以挑選最佳預測模型，并在盤古山鎢礦變形數據預測中取得了理想的效果。

針對露天礦邊坡下層采空區頂板位移超前預測的問題，本文以鞍鋼集團弓長嶺鐵礦獨木采區采空區頂板位移-時間序列為基礎，對K-means聚類算法進行改進，采用聚類-數理統計的組合模型實現對位移-時間序列數據的諸多誤差數據進行剔除，以保證監測結果的準確性，降低誤差數據對采空區頂板位移超前預測算法的干擾。在此基礎上，基于灰色系統理論實現對位移-時間序列數據進行超前預測，以提前預測采空區頂板位移量的變化，為礦山的安全生產和高效開采提供保障。

1 采空區頂板位移超前預測算法

首先應用層次迭代聚類算法實現雜亂數據的粗略檢測及剔除，使位移-時間序列數據能夠滿足正態分布；然后基于拉依達準則精細判斷隨機誤差數據，采用插值替補法替換隨機誤差數據；最后使用單一變量的一階灰色預測模型實現采空區頂板位移-時間序列數據的超前預測，并采用殘余誤差檢驗法、級比偏差檢驗法對預測模型進行檢驗，檢驗結果作為預測值的可信程度參考指標。技術路線如圖1所示。

圖1 技術路線 Fig. 1 Technology roadmap

1.1 誤差數據剔除組合模型

首先將誤差數據分為雜亂數據和隨機誤差數據兩大類，在此基礎上提出層次迭代聚類算法，實現誤差數據粗略檢測及剔除，然后基于拉依達準則實現隨機誤差數據的精細判別，并使用線性內插法替換判別出的隨機誤差數據。

1.1.1 層次迭代聚類算法

在對采空區頂板位移的實時監測中，不可避免地采集到諸多雜亂數據，雜亂數據一般具備以下2點較為明顯的特征：① 雜亂數據一般偏差較大；

通過對弓長嶺露天邊坡下層采空區頂板多點位移計NA2測點的已有歷史數據進行統計分析，將雜亂數據按照值大小粗略分為3級，見表1。

表1 雜亂數據分級標準 Table 1 Clutter data grading criteria

依據表1所示分級標準，對2019年12月31日采集的位移-時間序列進行分級，結果如圖2所示。

圖2 原始數據分級現象示意 Fig. 2 Schematic diagram of original data classification phenomenon

從圖2可以看出位移原始數據大部分集中于( 0，5 )的正常值范圍內，少量雜亂數據在表1所劃分的3級區間內均有分布，故可以采用聚類的方法分步剔除這些雜亂數據。由于雜亂數據中少量第3級數據“離群”過遠，傳統基于距離的聚類方法，如Kmeans聚類算法，會使部分“離群”不太遠的雜亂數據被計入正常值，導致雜亂數據剔除效果有限。為改善位移-時間序列數據中對雜亂數據的聚類效果，本文對K-means聚類算法進行改進，用于檢測并剔除雜亂數據，改進后的算法稱為層次迭代聚類算法。該算法的基本流程如下：

( 1 ) 對位移原始數據序列D0取每條數據的絕對值，得到絕對值序列D(0)。部分原始數據中雜亂數據既“離群”很遠，又有正有負，直接進行聚類會出現3個聚類中心，為了簡化聚類算法，不直接基于原始時間序列數據聚類，而是基于原始時間序列的絕對值序列進行聚類，這樣聚類算法始終都只有2個聚類中心。

( 2 ) 獲取絕對值序列D(0)的最大和最小值dmax，dmin，分別作為絕對值序列D(0)的2個簇中心，2個簇分別記為簇D1，簇D2。

( 3 ) 計算絕對值序列D(0)中每條數據d(i0)到2個簇的歐氏距離分別為Δ1，Δ2，并比較Δ1，Δ2的大小。Δ1計算公式(Δ2計算方法相同 )為

式中，Δ1為數據到簇D1的歐氏距離，mm；di(0)為第i個數據點的值，mm；dmin為簇D1的簇中心，mm。

( 4 ) 遍歷完絕對值序列D(0)并將所有數據進行聚類后，將簇D1中的成員記為數據序列D(1)。

( 5 ) 計算數據序列D(1)所有成員的均值，以及數據序列D(0)的均值，如果有，則結束聚類迭代，數據序列D(1)為剔除雜亂數據后的初步清洗序列，否則令D(0)＝D(1)，D(1)清空后循環執行步驟( 1 )～( 4 )。

1.1.2 拉依達準則

采空區頂板位移監測過程中，因電壓不穩等偶然因素產生的隨機誤差數據波動較小，層次分級不明顯，導致聚類時隨機誤差數據和正常值始終在同一個簇。經層次迭代聚類得到的數據序列D(1)中不含雜亂數據，滿足正態分布，故可以使用數理統計方法進行隨機誤差數據的檢測。檢測出的隨機誤差數據采用線性內插法進行插值替換。

應用拉依達準則( 也稱做3σ準則 )對數據序列D(1)中可能含有的隨機誤差數據進行精細判定。拉依達準則認為：在一個數據序列中，如果有數據點的殘余誤差大于數據序列標準差的3倍，則認為該數據點含有偶然誤差[16]。

拉依達準則精細判定隨機誤差數據的過程如下：

( 1 ) 計算數據序列D(1)的均值dabs：

式中，dabs為數據序列D(1)的均值，mm；di為數據序列D(1)中第i個數據點。

( 2 ) 計算數據序列D(1)的標準差σ：

( 3 ) 計算數據序列D(1)中每個數據點的殘差，并做3σ判定：

根據標準正態分布函數，滿足式( 4 )的數據點占數據序列D(1)的0.973；不滿足式( 4 )的數據點占數據序列D(1)的0.027[17]。對于不滿足拉依達準則的數據點，采用線性內插法進行替補，線性內插法公式為

1.1.3 誤差數據剔除效果檢驗

通過計算時間序列光滑度來檢驗誤差數據剔除效果。較為理想的誤差數據剔除效果是：上述2種隨機分布的誤差數據被剔除或替補后得到的時間序列中，任意相鄰2個數據點的值不會出現較大波動，時間序列曲線是平滑的。本文采用級比檢驗法檢驗位移-時間序列數據的光滑度。假設剔除誤差數據后的位移-時間序列為D0＝{d1，d2，…，dn}，則時間序列D0的級比檢驗公式為

式中，λi(0)為時間序列第i個數據的級比；di為第i個數據點，mm。

如果時間序列D0中任意相鄰2個數據點的級比滿足式( 7 )，那么可以判定位移-時間序列D0誤差數據剔除效果較好。

式中，n為時間序列長度。

1.2 位移-時間序列灰色預測模型

開采擾動、地質條件、地下水等諸多因素都會對采空區頂板位移產生作用，但是在對頂板位移量的預測研究中，不可能實現所有因素的量化，有些未知因素甚至無法進行量化描述，在這種信息量十分有限的情況下，可以將采空區頂板位移發展看作一個灰色系統。灰色系統是由我國學者鄧聚龍教授在20世紀80年代提出并不斷進行完善的一種慣性抽象系統[18]。灰色系統認為自然界有些事件在發展過程中，如果沒有受到外界環境突變影響，那么這一事件將會沿著一條不易改變的跡線產生慣性演變。灰色模型就是要基于有限的已知數據找出這條跡線。

本文對位移-時間序列使用單一變量的一階灰色模型進行超前預測，按照灰色理論，這種灰色模型被稱為GM( 1，1 )模型[19-20]。對位移-時間序列數據進行預測的基本原理：通過對剔除誤差數據后的位移-時間序列數據進行生成處理( 一般使用累加生成、指數生成等方法 )，使其形成呈指數發展趨勢的遞增序列。這種累加序列可實現原始離散數據的規律化，弱化灰色系統中的不確定性因素，在后續的遞減還原中又可還原原始數據中的特征數據。然后針對遞增序列建立一階差分方程，由一階差分方程變換得到參數未知的灰色預測模型。基于最小二乘原理近似計算出灰色模型計算參數，完成灰色模型建模。最后灰色預測模型通過遞減還原處理，得到模型值的重預測值及未來一段時間的預測值。

1.2.1 建立迭代灰色預測模型

位移-時間序列數據以( 時間，位移值 )鍵值對的形式存在，假設通過光滑度檢驗的位移絕對值時間序列為G(0)，做一次累加后得到的累加序列為G(1)，累加公式為

式中，gi(0)為第i個數據點，mm；gi(1)為第i個累加值，mm。

針對累加序列G(1)建立一階差分方程：

式中，a為灰色模型發展系數；b為灰色模型灰作用量。

基于最小二乘法基本原理，求解未知參數a，b。構造一個將a，b作為自變量的函數：

上式分別對a，b求偏導，求出使ω最小、擬合效果最優的a，b的值，代入式( 10 )整理得：

設通過最小二乘法得到的累加序列G(1)的擬合曲線為G(2)，則有灰色模型：

式( 12 )是累加序列G(1)的擬合曲線( 散點 )，對此擬合曲線做遞減還原，得到的時間序列即為原始時間序列G(0)的重預測值及超前預測值組成的預測時間序列G(3)：

式中，gi(3)為第i個預測值，mm。

通過以上過程可實現灰色模型建模，但灰色模型參數a，b是固定值，這就使得灰色模型不可應用于任意長度的時間序列。因為隨著時間序列步長的增加，參數a，b所確定的模型預測精度將會逐漸降低，直至模型失真。本文選擇時間序列的步長為10，即n=10，確保現在值與過去值的比例適中，再基于“先進先出”的原則，用現在值替補過去值，計算a，b的動態值，在迭代中及時修正灰色模型。迭代的灰色預測可以提高預測準確度，提供長期可信的穩定預測結果。

1.2.2 計算預測準確度指標

預測結果是否可信需要有直觀的評價指標作為參考。本文使用殘余誤差( 相對誤差 )檢驗法和級比偏差檢驗法對灰色預測模型的可信度進行評價。

殘余誤差檢驗法使用重預測值的平均相對誤差做為預測結果可信度評價指標。級比偏差可以判斷一組重預測值的光滑度。灰色模型能夠成立的假設是：在穩定環境下原始時間序列能夠按照一條平滑跡線發展。那么預測時間序列無論預測值是否準確，都應該首先滿足：預測時間序列也是一條光滑曲線。級比偏差檢驗法計算公式為

式中，λi(0)為原始時間序列的級比；λi(3)為重預測時間序列的級比；為平均級比偏差。

2 位移-時間序列超前預測微服務

基于Java語言完成算法程序的設計與開發，并創建Maven工程以方便項目管理。項目工程中使用Mybatis插件實現數據庫連接與數據讀取，使用Quartz插件實現算法程序的定時執行。最后使用Maven工具將項目打包成war包，直接運行于Tomcat服務器。

弓長嶺礦獨木采區因早期開采技術落后、盜采等原因形成的采空區一直未能得到有效治理。為確保采區生產安全，在已探明的采空區頂板布置多點位移計監測系統，對采空區頂板位移進行監測。每日采集的位移-時間序列數據最多達3 000條，其中隨機分布著諸多雜亂數據、隨機誤差。該位移監測系統自2018年6月開始運行，在運行初期階段，多點位移計運行狀況較好，監測值從0開始逐漸平穩增大，后在1～1.4 mm范圍內波動，監測數據質量較好。但由于長期受地表行車、降雨等影響，監測數據中逐漸出現“離群”較遠的誤差數據，位移監測值出現較大波動，導致監測數據序列產生整體偏移。多點位移計監測系統中NA2測點距離地表約13.61 m，距離采空區頂板約7.2 m，監測狀態穩定，能夠反映采空區頂板位移變化趨勢。故選取NA2測點的位移-時間序列數據，自2019年9月開始進行長期預測。

將存儲NA2測點位移原始數據的數據表接入到數據預測微服務程序中，借助Java Quartz定時任務插件，于每日0時讀取前10 d的位移監測值，通過誤差數據剔除程序及灰色預測程序得到新一天的采空區頂板位移預測值。

2.1 NA2測點剔除誤差數據

NA2測點平均每日可采集230條數據，若不進行誤差數據的識別與剔除，直接對原始數據求平均值，得到的結果如圖3所示。

圖3 NA2測點原始數據 Fig. 3 Original data of NA2 measuring point

由圖3可以看出，2019-09—2020-02期間，NA2測點處頂板位移在1.5～1.7 mm范圍內波動，由于受到誤差數據的干擾，監測數據中夾雜有較多“離群”較遠的散點，導致監測數據無法反映采空區頂板真實的位移變化。直接使用未經處理的原始數據進行預測，將會使預測結果失真。本文使用迭代聚類+數理統計的組合模型對誤差數據進行處理，其程序設計如圖4所示。

圖4 誤差剔除程序設計 Fig. 4 Eliminate error programming

在針對NA2測點位移-時間序列的微服務程序中，程序自動剔除NA2測點的誤差數據，得到NA2測點剔除誤差數據后的位移-時間序列，如圖5所示。

圖5 NA2測點位移-時間序列誤差剔除效果 Fig. 5 Error elimination effect of NA2 measuring point displacement-time series

在剔除誤差數據后，NA2測點位移監測值在1.4～1.6 mm范圍內波動，整體呈現增大趨勢。對該曲線進行光滑度檢驗，根據式( 6 )計算得到剔除誤差數據之后位移-時間序列的級比如圖6所示，符合由式( 7 )計算得到的級比區間( 0.833 8，1.199 4 )，說明誤差數據剔除效果較好。

圖6 光滑度檢驗結果 Fig. 6 Inspection result of smoothness

同時該曲線符合監測點處于采空區頂板中較密實位置，受外界環境影響相對較小，采空區頂板應力處于相對平衡狀態，故在很長一段監測時間內，測點處累積位移量僅有輕微增加。說明本文通過迭代聚類+數理統計的組合模型，能有效剔除位移-時間序列中的雜亂數據、隨機誤差數據，剔除誤差數據后的位移-時間序列可以較準確地反映采空區頂板實際監測位移量。綜合考量以上定量計算及定性分析結果，可以認為誤差數據剔除算法合理有效。

2.2 NA2測點灰色預測

使用剔除誤差數據后的位移-時間序列執行灰色預測程序進行超前預測。預測程序首先得到位移-時間序列的累加序列。引入累加序列可以很好地將對波動較大的原始位移-時間序列曲線的擬合轉化為對累加序列曲線的擬合，即將弱規律或無規律的復雜曲線擬合問題轉化為單調遞增的簡單曲線擬合，擬合難度大大降低。然后基于最小二乘原理得到累加序列的有參擬合函數，節選部分位移-累加序列如圖7所示。

由圖7可以看出，原本復雜、波動的位移-時間序列得到了一個可以使用簡單函數進行精準擬合的累加序列。基于此累加序列的擬合結果進行遞減還原，得到采空區頂板歷史位移數據的重預測值及未來預測值。預測程序在2019-09-01—2020-02-10期間運行得到的預測結果如圖8所示。

圖7 節選部分位移累加序列 Fig. 7 Extract part of the displacement accumulation sequence

從圖8可以看出剔除誤差數據后，灰色預測模型預測效果較好，對采空區頂板NA2測點位移自2019年9月—2020年2月的預測中，預測值與實測值基本吻合。得到的預測結果能夠反映采空區頂板位移量的變化。分別采用殘余誤差檢驗法和級比偏差檢驗法對預測結果的可信度進行檢驗，NA2測點位移-時間序列超前預測微服務自運行以來，預測值的殘余誤差率均值為0.22%，最大不超過0.96%，預測值的級比偏差均值為0.21%，最大不超過0.89%。

圖8 位移-時間序列預測結果 Fig. 8 Predictive results of displacement-time series

3 結 論

( 1 ) 結合位移-時間序列數據特征，在時間序列數據中含有誤差數據的情況下，提出層次迭代聚類+拉依達準則組合算法，實現誤差數據的自動檢測及剔除。通過改進K-means聚類算法得到的層次迭代聚類算法可以有效剔除位移-時間序列數據中“離群”較遠的雜亂數據，為確保超前預測程序的長期穩定運行提供了保障。

( 2 ) 基于灰色系統理論，對采空區頂板位移-時間序列數據進行長期超前預測分析。采空區頂板位移的歷史預測結果與對應的實測值基本吻合，預測值的殘余誤差和級比偏差均在1%以內，能夠正確預測弓長嶺露天邊坡下層采空區頂板位移量的變化趨勢，預測結果與實際采集值具有對比意義。

( 3 ) 從近6個月的位移預測結果及實測值可以看出，弓長嶺露天邊坡下層采空區頂板位移不明顯，位移量呈現出極緩慢的增加趨勢( 6個月內采空區頂板僅產生0.2 mm左右的位移 )，短期內沒有發生失穩破壞的風險。