王昕宇，楊 鵬，2，戴健非

(1.北京聯(lián)合大學(xué)北京市信息服務(wù)工程重點(diǎn)實(shí)驗(yàn)室，北京 100101；2.北京科技大學(xué)土木與資源工程學(xué)院，北京 100083)

0 引言

尾礦庫是選擇有利地形筑壩攔截谷口或圍地形成的具有一定容積、用以貯藏尾礦和尾礦澄清水的專用場地，是一種高勢能的人造泥石流危險(xiǎn)源[1].尾礦壩是使尾礦庫形成一定容積，便于尾礦礦漿能堆存其中的一種由尾礦堆積碾壓而成的壩體，尾礦壩一旦潰壩，會(huì)嚴(yán)重威脅下游居民的生命安全，污染下游區(qū)域環(huán)境，造成大量經(jīng)濟(jì)財(cái)產(chǎn)損失.其安全穩(wěn)定對(duì)礦山生產(chǎn)、安全、環(huán)境保護(hù)至關(guān)重要[2].目前，自動(dòng)化系統(tǒng)雖然能夠?qū)ξ驳V壩主要指標(biāo)進(jìn)行監(jiān)測，但安全監(jiān)測僅是初步，能夠根據(jù)以往監(jiān)測數(shù)據(jù)建立預(yù)測模型，了解其形變規(guī)律才是最重要的.因此，將尾礦壩壩體位移數(shù)據(jù)與數(shù)學(xué)模型相結(jié)合實(shí)現(xiàn)預(yù)測預(yù)警，對(duì)尾礦庫的安全穩(wěn)定管理具有重要意義[3].

我國對(duì)于尾礦壩壩體位移預(yù)測起步相對(duì)較晚，借鑒了大量的大壩(水壩、堆石壩等)位移預(yù)測體系方法.目前國內(nèi)外研究者圍繞壩體位移預(yù)測模型展開大量的研究，其中黃定川等[4]根據(jù)尾礦庫初期壩沉降數(shù)據(jù)為例構(gòu)建BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型，重點(diǎn)研究其拓?fù)浜蛯W(xué)習(xí)算法；胡軍等[5]提出使用基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化的SVM算法進(jìn)行預(yù)測學(xué)習(xí)，并通過種族魚群優(yōu)化算法優(yōu)化參數(shù)；李豐旭等[6]針對(duì)形變模型發(fā)生早熟、預(yù)測不穩(wěn)定的現(xiàn)象，引入基于混沌搜索的多種群遺傳算法，改進(jìn)遺傳參數(shù)的計(jì)算方式；甄明秀[7]以白巖尾礦壩為媒介，比較回歸模型、時(shí)序模型、灰色模型在形變預(yù)測中的適應(yīng)能力；賈龍等[8]針對(duì)verhulst模型預(yù)測精度不足的問題，使用三類弱化緩沖算子處理數(shù)據(jù)，重新構(gòu)建verhulst模型；鐘登華等[9]提出了基于自適應(yīng)網(wǎng)絡(luò)模糊推理系統(tǒng)優(yōu)化的灰色模型；文獻(xiàn)[10]提出了一種新的分布式時(shí)間序列演化預(yù)測方法，將變形預(yù)測的思想與結(jié)構(gòu)化方法相結(jié)合，提高了變形預(yù)測的可靠性；文獻(xiàn)[11]針對(duì)回歸模型中擬合殘差序列的混沌特性，采用多尺度小波分析對(duì)殘差序列進(jìn)行分解和重構(gòu)，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型.

上述研究對(duì)壩體的安全穩(wěn)定管理有著重要意義，但并未充分涉及外部因素對(duì)形變的影響，僅考慮位移變化的時(shí)效性.合理的特征選擇能夠提升模型精度，降低時(shí)間復(fù)雜度.因此，本文提出建立特征遞歸消除(Recursive Feature Elimination，RFE)和隨機(jī)森林(Random Forest，RF)、極限梯度提升(Extreme Gradient Boosting，XGBoost)相融合的預(yù)測模型，并與BP神經(jīng)網(wǎng)絡(luò)、LSTM神經(jīng)網(wǎng)絡(luò)等模型進(jìn)行了對(duì)比.本文對(duì)防止?jié)问鹿屎偷V山安全管理具有一定參考價(jià)值.

1 壩體位移模型構(gòu)建

1.1 隨機(jī)森林

隨機(jī)森林是Cutler與Breiman于2001年提出的一種統(tǒng)計(jì)學(xué)習(xí)方法[12].RF是根據(jù)一定準(zhǔn)則生成一定數(shù)量的決策樹并進(jìn)行組合，在訓(xùn)練決策樹的過程中隨機(jī)選擇特征以提高模型預(yù)測精度.可用于模擬多重非線性關(guān)系，該算法的實(shí)現(xiàn)步驟如下：

(1) 對(duì)含有m個(gè)樣本的數(shù)據(jù)集，使用自助法從數(shù)據(jù)集中隨機(jī)有放回的選取m個(gè)樣本.經(jīng)過N次選取，產(chǎn)生N個(gè)bootstrap樣本集，對(duì)N個(gè)樣本集訓(xùn)練構(gòu)成N個(gè)回歸樹，而每次選取后的剩余樣本全部放入N個(gè)袋外數(shù)據(jù)(OBB)作為測試樣本.

(3) 每個(gè)回歸樹自上往下遞歸分支，設(shè)定隨機(jī)森林中回歸樹的個(gè)數(shù)作為回歸樹生長的停止條件.

(4) 綜合N個(gè)回歸樹并取其平均值構(gòu)成隨機(jī)森林回歸模型，公式為

(1)

式中：f(x)表示隨機(jī)森林模型，Tj(x)是單個(gè)回歸樹模型，j=1，2，…，N，使用OBB中的樣本作為測試集并根據(jù)回歸分析的相關(guān)評(píng)估指標(biāo)來判斷模型的預(yù)測效果.

自助法是一種非參數(shù)統(tǒng)計(jì)方法，也被稱作bootstrap，其基本步驟如下：

(1) 從樣本中抽取一定量樣本，允許重復(fù)抽樣(抽取數(shù)量自擬).

(2) 統(tǒng)計(jì)抽取的樣本以計(jì)算待估計(jì)統(tǒng)計(jì)量M.

(3) 重復(fù)抽取N次，計(jì)算N個(gè)統(tǒng)計(jì)量M.

(4) 根據(jù)N個(gè)統(tǒng)計(jì)量M的方差估計(jì)統(tǒng)計(jì)量M的方差.

1.2 特征遞歸消除

特征遞歸消除是一種基于wrapper包裹型模式下的后向搜索算法，常用于特征選擇[13].相較于filter型算法，wrapper將特征選擇問題轉(zhuǎn)化為特征子集搜索問題，篩選多個(gè)特征子集，使用一個(gè)機(jī)器學(xué)習(xí)模型評(píng)估效果，重復(fù)該過程直至產(chǎn)生最優(yōu)的特征子集，RFE算法是在該方法基礎(chǔ)上的簡化，基本步驟如下：

(1) 選擇全部特征并構(gòu)造初始特征集.

(2) 在特征集上訓(xùn)練一個(gè)機(jī)器學(xué)習(xí)模型，本文采用隨機(jī)森林算法構(gòu)建回歸模型.

(3) 依據(jù)平均不純度的減少程度評(píng)估各特征的重要性，并由高到低排列.回歸樹模型以方差或最小二乘法作為不純度條件，公式為：

(2)

(3)

(4) 剔除重要性最低的特征，構(gòu)造新的特征集.

(5) 重復(fù)步驟2到步驟4，直到特征集沒有進(jìn)一步的特征可供保持，特征被剔除的順序是特征選擇的依據(jù)，越早被剔除的特征越應(yīng)舍棄，越晚被剔除的特征越應(yīng)保留.

1.3 極限梯度提升

極限梯度提升是陳天奇博士于2016年提出的一種集成學(xué)習(xí)方法，以CART回歸樹和隨機(jī)森林為基礎(chǔ)的一種擴(kuò)展延伸，屬于Boosting方法中的一種提升算法[14].Boosting是將多個(gè)弱學(xué)習(xí)器提升為強(qiáng)學(xué)習(xí)器的一種方法，先訓(xùn)練第1個(gè)弱學(xué)習(xí)器，根據(jù)弱學(xué)習(xí)器的誤差表現(xiàn)改變訓(xùn)練樣本權(quán)重，使訓(xùn)練誤差大的樣本受到更多的關(guān)注.將改變權(quán)重的訓(xùn)練樣本導(dǎo)入第2個(gè)弱學(xué)習(xí)器，再根據(jù)訓(xùn)練誤差調(diào)整權(quán)重，如此反復(fù)，直到弱學(xué)習(xí)器的數(shù)量達(dá)到某一指定值n，將n個(gè)弱學(xué)習(xí)器整合為強(qiáng)學(xué)習(xí)器.XGBoost在Boosting的基礎(chǔ)上使用增量學(xué)習(xí)并增加L1和L2正則化項(xiàng)，提升泛化能力，減小過擬合，基本算法如下：

(4)

目標(biāo)函數(shù)為

(5)

1.4 RFE-RF-XGBoost預(yù)測模型構(gòu)建

RFE-RF-XGBoot預(yù)測模型構(gòu)建流程如下：

(1) 數(shù)據(jù)集整理，刪除缺失值，數(shù)據(jù)做最小-最大標(biāo)準(zhǔn)化(Min-Max Normalization)處理.

(2) 壩體垂直沉降作為輸出，其余特征作為輸入，導(dǎo)入RFE-RF模型進(jìn)行特征選擇.

(3) 按特征重要度重建數(shù)據(jù)集，采用留出法交叉驗(yàn)證(Holdout Cross Validation)，將數(shù)據(jù)依照時(shí)間順序劃分為訓(xùn)練集(60%)、驗(yàn)證集(20%)、測試集(20%).

(4) 設(shè)定XGBoost模型的初始參數(shù)對(duì)訓(xùn)練集擬合，利用驗(yàn)證集計(jì)算誤差指標(biāo)并配合網(wǎng)格搜索與模型歷史最優(yōu)指標(biāo)對(duì)比.

(5) 若模型預(yù)測精度提高，則當(dāng)前指標(biāo)為新的歷史最優(yōu)指標(biāo)，繼續(xù)調(diào)整參數(shù)，以此反復(fù).反之，若模型預(yù)測精度低于歷史最優(yōu)指標(biāo)，則說明模型出現(xiàn)過擬合現(xiàn)象，停止調(diào)整參數(shù).

(6) 測試樣本導(dǎo)入模型計(jì)算預(yù)測誤差，查看模型預(yù)測效果.基本流程如圖1所示.

1.5 模型評(píng)估指標(biāo)

本文使用平均絕對(duì)誤差(Mean Absolute Error，MAE)和相對(duì)誤差(Relative Error，RE)作為模型的評(píng)估指標(biāo)，MAE也被稱作平均絕對(duì)離差即絕對(duì)誤差的平均值，是參數(shù)估計(jì)值與參數(shù)真值之差絕對(duì)值的期望值.MAE可以評(píng)價(jià)模型預(yù)測的絕對(duì)準(zhǔn)確度，MAE的值越小，說明預(yù)測模型描述整體實(shí)驗(yàn)數(shù)據(jù)具有更好的精確度，模型更穩(wěn)定[16].其公式為

(6)

相對(duì)誤差RE是指被測量的真實(shí)值與預(yù)測值的絕對(duì)誤差與被測量真實(shí)值之比乘以100%所得的數(shù)值，以百分?jǐn)?shù)表示.其公式為

(7)

圖1 RFE-RF-XGBoost模型構(gòu)建流程

2 實(shí)例應(yīng)用

2.1 數(shù)據(jù)來源

本文采用某尾礦庫安全監(jiān)測系統(tǒng)采集的2017年11月至2019年8月共7 000條實(shí)時(shí)數(shù)據(jù)進(jìn)行預(yù)測實(shí)驗(yàn).該尾礦庫規(guī)劃用地面積約1.76 km2.尾礦庫區(qū)匯水面積3.91 km2，屬于二級(jí)尾礦庫，尾礦庫初期壩屬透水堆石壩，初期壩設(shè)計(jì)標(biāo)高為450.00 m，高度約35 m，最終堆積標(biāo)高為647.00 m.尾礦堆積壩采用上游式堆筑法，即平行初期壩軸方向外尾礦庫內(nèi)逐級(jí)填筑堆積壩，增加庫容排放尾礦.尾礦堆積壩每級(jí)高度2 m，內(nèi)外坡比1∶2∶5，頂寬2 m，堆積壩總高度為197 m.特征包括3個(gè)監(jiān)測點(diǎn)處的浸潤線埋深(JRX1-1、JRX1-2、JRX1-3)、庫水位高度(KSW)、降雨量(RAIN)、壩體表面平行壩軸線位移(BMWY-X)、垂直壩軸線位移(BMWY-Y)、垂直沉降(BMWY-Z)，部分樣本特征數(shù)據(jù)如表1所示.

表1 部分實(shí)驗(yàn)樣本特征

2.2 特征選擇

憑借PyCharm搭配Anaconda開源包、環(huán)境管理器，利用python3語言編寫.以3點(diǎn)浸潤線埋深、庫水位高度、降雨量、壩體表面平行壩軸線位移、垂直壩軸線位移作為輸入，壩體垂直沉降作為輸出.首先調(diào)用preprocessing模塊的scale函數(shù)對(duì)輸入輸出特征進(jìn)行快速M(fèi)in-Max標(biāo)準(zhǔn)化，然后導(dǎo)入ensemble模塊中的RandomForestRegressor函數(shù)，默認(rèn)決策樹數(shù)值為100并擬合數(shù)據(jù)，計(jì)算輸入特征重要性.結(jié)果表明庫水位重要性占比為0.399，其次為兩點(diǎn)浸潤線埋深，分別占0.184和0.176.壩體表面平行壩軸線位移、垂直壩軸線位移分別占0.065和0.092.降雨量占比最低，為0.020，第一次篩選結(jié)果如圖2所示.

圖2 第一次特征篩選

待剔除降雨量特征后，數(shù)據(jù)重新導(dǎo)入模型并計(jì)算重要性，重要性占比發(fā)生變化：浸潤線埋深1-2和1-3占比降低，分別為0.063和0.173，第二次篩選結(jié)果如圖3所示.

經(jīng)多次篩選，重要性排名從高到低依次是庫水位、浸潤線埋深1-1、浸潤線埋深1-3、垂直壩軸線位移、平行壩軸線位移、浸潤線埋深1-2、降雨量.其中庫水位重要性最高，通過mic函數(shù)計(jì)算相關(guān)性為0.530，表明庫水位與垂直沉降具有強(qiáng)非線性相關(guān)性.本文選擇浸潤線埋深1-1、浸潤線埋深1-3、庫水位、垂直沉降、降雨量為特征重建新數(shù)據(jù)集，為預(yù)測t時(shí)刻的垂直沉降，以t-dt時(shí)刻的浸潤線埋深、庫水位、垂直沉降、降雨量為輸入，其中dt為1，2，3 d，新數(shù)據(jù)集部分樣本如表2所示，通過特征選擇可有效減少預(yù)測過程中的特征冗余，避免特征間的負(fù)面影響.

圖3 第二次特征篩選

表2 新數(shù)據(jù)集部分樣本特征

2.3 RFE-RF-XGBoost模型建立

根據(jù)RFE-RF特征選擇排列重建新數(shù)據(jù)集后，數(shù)據(jù)集采用留出法交叉驗(yàn)證，依照時(shí)間順序劃分為訓(xùn)練集(60%)、驗(yàn)證集(20%)、測試集(20%)，其中訓(xùn)練集導(dǎo)入XGBoost框架的XGBRegressor函數(shù)中訓(xùn)練擬合，驗(yàn)證集為模型泛化能力驗(yàn)證指標(biāo)，測試集用于檢驗(yàn)?zāi)Ｐ皖A(yù)測效果.XGBoost模型參數(shù)主要包括：回歸樹數(shù)量(n_estimators)、回歸樹最大深度(max_depth)、回歸樹最小葉子節(jié)點(diǎn)樣本權(quán)重和(min_child_weight)、損失函數(shù)減小閾值(gamma)、回歸樹隨機(jī)采樣率(subsample)、回歸樹特征采樣率(colsample_bytree)、L1正則化權(quán)重(reg_alpha)、L2正則化權(quán)重(reg_lambda)、學(xué)習(xí)率(learning_rate).本文根據(jù)驗(yàn)證集預(yù)測指標(biāo)，配合網(wǎng)格搜索方法逐步調(diào)整模型參數(shù)，過程如下：

(1) 調(diào)整回歸樹數(shù)量，n_estimators從[10，300]之間取值計(jì)算驗(yàn)證集誤差指標(biāo)，結(jié)果如圖4所示.[10，40]區(qū)間內(nèi)誤差曲線呈下降趨勢，當(dāng)n_estimators=40時(shí)，誤差最小，MAE=0.072，[50，100]區(qū)間內(nèi)誤差曲線逐漸增高，泛化能力降低.

圖4 回歸樹數(shù)量對(duì)MAE的影響

(2) 調(diào)整回歸樹最大深度和最小葉子節(jié)點(diǎn)樣本權(quán)重和，max_depth調(diào)整區(qū)間為[0，10]，min_child_weight調(diào)整區(qū)間為[0，20]，結(jié)果如圖5所示.max_depth=6，min_child_weight=19時(shí)，誤差最小，MAE=0.068.

(3) 調(diào)整損失函數(shù)減小閾值，gamma調(diào)整區(qū)間為[0，1]，gamma=0.05時(shí)，誤差最小，MAE=0.067.

(4) 調(diào)整回歸樹隨機(jī)采樣率和回歸樹特征采樣率，subsample和colsample_bytree調(diào)整區(qū)間為[0，1]，subsample=0.5，colsample_bytree=0.9時(shí)，誤差最小，MAE=0.065.

圖5 回歸樹最大深度與最小葉子節(jié)點(diǎn)

(5) 調(diào)整L1正則化權(quán)重和L2正則化權(quán)重，reg_alpha和reg_lambda調(diào)整區(qū)間為[0，7]，reg_alpha=0.06，reg_lambda=0.1時(shí)，誤差最小，MAE=0.062.

(6) 調(diào)整學(xué)習(xí)率，learning_rate調(diào)整區(qū)間為[0，1]，learning_rate=0.2時(shí)，誤差最小，MAE=0.059.

經(jīng)多輪參數(shù)調(diào)整，最小平均絕對(duì)誤差MAE=0.059.表3為驗(yàn)證集的10組樣本預(yù)測結(jié)果和相對(duì)誤差，可知垂直沉降預(yù)測最大相對(duì)誤差為7.21%，最小相對(duì)誤差為0.17%.

表3 RFE-RF-XGBoost對(duì)驗(yàn)證集樣本預(yù)測結(jié)果

2.4 不同方法預(yù)測結(jié)果比較

本文使用tesorflow神經(jīng)網(wǎng)絡(luò)框架對(duì)BP神經(jīng)網(wǎng)絡(luò)編寫及測試，網(wǎng)絡(luò)結(jié)構(gòu)包括輸入層、隱藏層、輸出層.其中隱藏層激活函數(shù)為tanh函數(shù)，輸出層激活函數(shù)為線性函數(shù).輸入層神經(jīng)單元數(shù)為24，壩體表面垂直沉降為輸出，輸出單元個(gè)數(shù)為1.使用SGD優(yōu)化器實(shí)現(xiàn)權(quán)重、閾值更新.學(xué)習(xí)率為0.01，Dropout為0.2，迭代3 000次.BP神經(jīng)網(wǎng)絡(luò)的隱藏層神經(jīng)元個(gè)數(shù)一般是根據(jù)以往經(jīng)驗(yàn)和反復(fù)的試驗(yàn)選取，為尋找預(yù)測效果最好的網(wǎng)絡(luò)結(jié)構(gòu)，本文通過Kolmogorov定理計(jì)算隱含層神經(jīng)元個(gè)數(shù)，公式為

h=2m+1.

(8)

其中：m為輸入層神經(jīng)元個(gè)數(shù)，h為隱含層神經(jīng)元個(gè)數(shù)，網(wǎng)絡(luò)結(jié)構(gòu)為24-49-1時(shí)，預(yù)測誤差MAE=0.133.

作為傳統(tǒng)深度學(xué)習(xí)模型代表之一，LSTM神經(jīng)網(wǎng)絡(luò)模型通過keras深度學(xué)習(xí)框架搭建，網(wǎng)絡(luò)結(jié)構(gòu)包括輸入層、LSTM層、全連接層、輸出層[18].輸入層神經(jīng)元24個(gè)，LSTM層神經(jīng)元60個(gè)，輸出層神經(jīng)元1個(gè).LSTM層激活函數(shù)為tanh，輸出層激活函數(shù)為線性函數(shù).優(yōu)化器是Adam，學(xué)習(xí)率為0.01.Dropout=0.2，batch size=64，epochs=46時(shí)，預(yù)測誤差MAE=0.078.SVR預(yù)測模型通過sklearn模塊的SVR函數(shù)構(gòu)建，算法使用高斯核函數(shù)，正則化系數(shù)為0.8時(shí)，MAE=0.080[19].圖6與表4為RFE-RF-XGBoost、XGBoost、LSTM神經(jīng)網(wǎng)絡(luò)、BP神經(jīng)網(wǎng)絡(luò)、SVR預(yù)測結(jié)果對(duì)比.

表4 多方法預(yù)測結(jié)果

圖6 多方法預(yù)測結(jié)果對(duì)比

表4為各模型樣本預(yù)測結(jié)果和相對(duì)誤差，由表4可知：RFE-RF-XGBoost模型最小相對(duì)誤差為1.01%，最大相對(duì)誤差為6.79%.XGBoost模型最小相對(duì)誤差為3.02%，最大相對(duì)誤差為12.84%.各模型MAE依次為0.060，0.062，0.080，0.133，0.078，平均相對(duì)誤差分別為3.60%，7.53%，11.20%，16.76%，9.25%，其中RFE-RF-XGBoost預(yù)測效果最好，平均相對(duì)誤差低于XGBoost3.93%.

3 結(jié)論

本文以某尾礦庫2017年11月至2019年8月的歷史數(shù)據(jù)為實(shí)驗(yàn)樣本，提出使用隨機(jī)森林配合特征遞歸消除與極限梯度提升模型相融合.實(shí)驗(yàn)結(jié)果表明：庫水位與垂直沉降重要性最高，非線性相關(guān)性為0.530.RFE-RF-XGBoost對(duì)比XGBoost模型平均絕對(duì)誤差降低了0.2%，平均相對(duì)誤差降低了3.93%.所提出模型能夠提高壩體位移的預(yù)測精度，總體精度高于XGBoost模型、SVR模型、BP神經(jīng)網(wǎng)絡(luò)、LSTM神經(jīng)網(wǎng)絡(luò).借助RFE-RF-XGBoost預(yù)測模型，礦山管理部門可依據(jù)預(yù)測結(jié)果及早實(shí)施相應(yīng)安全措施，對(duì)于礦山施工決策、安全管理，減少潰壩事故具有參考意義.