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

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

0 引言

尾礦庫(kù)是維持礦山企業(yè)正常運(yùn)營(yíng)的必要設(shè)施，由1個(gè)或多個(gè)尾礦壩堆筑攔截谷口或圍地所構(gòu)成的礦山生產(chǎn)設(shè)施，是具有高勢(shì)能的人造危險(xiǎn)源。其穩(wěn)定性與安全性對(duì)礦山企業(yè)的正常生產(chǎn)和區(qū)域安全至關(guān)重要[1-3]。因此，研究尾礦壩潰壩失效的早期預(yù)警是很有必要的。在Fund?o尾礦壩倒塌1 年后，超過4 300萬(wàn)m3的鐵礦石尾礦仍然造成環(huán)境破壞，污染了從Doce河流域到大西洋的668 km水道。污染物的數(shù)量和受影響的生態(tài)系統(tǒng)的范圍占據(jù)了前所未有的比例，涉及巴西大西洋森林河口，沿海和海洋環(huán)境[4]。對(duì)潰壩原因調(diào)查后發(fā)現(xiàn)[5]，尾礦壩由于一系列“施工缺陷”導(dǎo)致了砂土液化。其中，基礎(chǔ)排水溝的施工缺陷最為嚴(yán)重，排水系統(tǒng)有缺陷，導(dǎo)致庫(kù)內(nèi)尾礦飽和度過高，小型地震誘發(fā)潰壩。由此可見，排水能力是否正常在尾礦壩安全領(lǐng)域起著重要的指示作用。本文根據(jù)專家小組的數(shù)據(jù)支持預(yù)測(cè)排水?dāng)?shù)據(jù)flow即每秒排水量(L/s)。國(guó)內(nèi)外學(xué)者在尾礦庫(kù)安全預(yù)警方面做了大量研究。在尾礦壩數(shù)據(jù)預(yù)測(cè)方面， Dong等[6]通過考慮潛水線預(yù)測(cè)方程，極限平衡狀態(tài)參數(shù)，水庫(kù)水位和降雨量建立數(shù)值模擬模型，通過遠(yuǎn)程實(shí)時(shí)預(yù)警獲得尾礦壩穩(wěn)定或危險(xiǎn)的預(yù)警信號(hào)；王英博等[7]采用修正型果蠅算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)對(duì)灘頂高度、庫(kù)水位、浸潤(rùn)線、干灘長(zhǎng)度和安全超高5個(gè)指標(biāo)擬合,顯示出較高的預(yù)測(cè)精度；李娟等[8]使用SVM預(yù)測(cè)尾礦壩浸潤(rùn)線,實(shí)現(xiàn)了小樣本情況下的高精度預(yù)測(cè)； Tongle等[9]基于遺傳算法和BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)尾礦浸潤(rùn)線的未來(lái)狀態(tài)。在尾礦壩安全理論方面，何學(xué)秋等[10]經(jīng)過多次實(shí)驗(yàn)得出尾礦壩變形需經(jīng)歷衰減、穩(wěn)定、加速3個(gè)階段的結(jié)論，基于流變-突變理論,分別制定預(yù)警準(zhǔn)則；王昆等[11]綜合分析了國(guó)內(nèi)外尾礦庫(kù)潰壩災(zāi)害防控現(xiàn)狀及發(fā)展情況；張興凱等[12]采用物理模型試驗(yàn)方法建立尾礦庫(kù)漫頂潰壩演化模型，證明了尾礦庫(kù)洪水漫頂潰壩位移與壩體飽和程度有關(guān)，壩體浸潤(rùn)線越高,尾礦庫(kù)潰壩時(shí)滑動(dòng)位移越大，潰口破壞程度取決于溢流對(duì)壩體的沖刷侵蝕作用。在傳感器穩(wěn)定采集數(shù)據(jù)方面，陳凱等[13]研究了礦山監(jiān)測(cè)系統(tǒng)的防護(hù)、供電能力、網(wǎng)絡(luò)通信等技術(shù),來(lái)確保極端天氣環(huán)境下的數(shù)據(jù)獲取。

為實(shí)現(xiàn)尾礦庫(kù)排水流量數(shù)據(jù)的實(shí)時(shí)預(yù)測(cè)，本文建立基于BP神經(jīng)網(wǎng)絡(luò)(Back Propagation Neural Network)和支持向量回歸(Support Vector Regression， SVR)的巴西Fund?o尾礦壩流量數(shù)據(jù)的預(yù)測(cè)回歸模型。分析庫(kù)水位、降雨量和干灘長(zhǎng)度對(duì)排水?dāng)?shù)據(jù)的影響。本文所述模型可為Fund?o壩的潰壩預(yù)警提供決策支持，對(duì)防止?jié)问鹿屎捅Ｗo(hù)尾礦安全具有參考價(jià)值。

1 數(shù)據(jù)來(lái)源與傳感器

本文所用數(shù)據(jù)均來(lái)自于文獻(xiàn)[5]。

巴西Fund?o尾礦壩的流量監(jiān)測(cè)是采用位于海拔826 m處的流量傳感器采集得到。在826 m處安裝有27個(gè)排洪隧洞進(jìn)行排水，材料為混凝土材質(zhì)。2011年的ITRB(Independent Tailings Review Board)報(bào)告指出，排水系統(tǒng)朝向右側(cè)和左側(cè)基臺(tái)，但是，只能從右側(cè)基臺(tái)附近的水槽進(jìn)行測(cè)量[5]。

排水?dāng)?shù)據(jù)監(jiān)測(cè)有每小時(shí)排水立方米和每秒排水升數(shù)2種計(jì)量模式。本文僅采用每秒排水升數(shù)作為模型的輸出特征，并記為flow。庫(kù)水位和降雨量來(lái)自于小組審查的Samarco現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)。安裝的監(jiān)測(cè)儀器類型是振弦式壓力計(jì)和水位指示器(立管)。主要數(shù)據(jù)來(lái)源是Samarco的月度儀表報(bào)告和Samarco電子數(shù)據(jù)表格。干灘長(zhǎng)度數(shù)據(jù)每月更新1次。數(shù)據(jù)來(lái)源共有4種：月度尾礦排放報(bào)告、使用無(wú)人機(jī)地形測(cè)量估算干灘長(zhǎng)度及上升率電子表格、每月巖土監(jiān)測(cè)報(bào)告和地形測(cè)量數(shù)據(jù)。選取的潰壩災(zāi)害數(shù)據(jù)包含3個(gè)特征屬性：庫(kù)水位(m)、降雨量(mm)、干灘長(zhǎng)度(m)。本文提取了2011年7月至2015年11月的巴西Fund?o壩的上述3種潰壩災(zāi)害特征數(shù)據(jù)和排水?dāng)?shù)據(jù)共62組數(shù)據(jù)作為數(shù)據(jù)源。其中，2011年7月至 2015年3月的數(shù)據(jù)作為訓(xùn)練集，部分歷史數(shù)據(jù)如表1所示。

表1 Fund?o壩部分歷史數(shù)據(jù)Table 1 Partial historical datas of the Fund?o

2 基于BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型

人工神經(jīng)網(wǎng)絡(luò)技術(shù)是對(duì)大腦功能的模擬、抽象和簡(jiǎn)化，具有學(xué)習(xí)、聯(lián)想和記憶功能，以及高度的非線性預(yù)測(cè)能力[14]。本文采用基于Levenberg-Marquardt算法的BP雙隱層神經(jīng)網(wǎng)絡(luò)對(duì)巴西Fund?o壩的排水?dāng)?shù)據(jù)進(jìn)行預(yù)測(cè)建模。該雙隱層網(wǎng)絡(luò)包含1個(gè)輸入層，1個(gè)輸出層，2個(gè)隱藏層。其中輸入層有3個(gè)神經(jīng)元節(jié)點(diǎn)，輸出層包含1個(gè)神經(jīng)元節(jié)點(diǎn)。為了全面分析巴西Fund?o壩的排水流量規(guī)律，本文提取了2011年7月至2015年3月的庫(kù)水位、降雨量和干灘長(zhǎng)度數(shù)據(jù)做訓(xùn)練集的輸入，而對(duì)應(yīng)的排水流量則為輸出特征。在BP神經(jīng)網(wǎng)絡(luò)模型中，回歸系數(shù)R用于評(píng)估模型的有效性，R越接近1，模型的預(yù)測(cè)能力越精準(zhǔn)。

2.1 Levenberg-Marquardt算法

Levenberg-Marquardt算法結(jié)合了高斯牛頓法和梯度下降法的優(yōu)點(diǎn)，是目前使用廣泛的優(yōu)化算法。設(shè)x(i)為當(dāng)前狀態(tài)的權(quán)重向量，x(i+1)為BP神經(jīng)網(wǎng)絡(luò)中更新后的權(quán)重向量，Δx為權(quán)重更新值，有：

x(i+1)=x(i)+Δx

(1)

Δx=-[2E(x)]-1E(x)

(2)

(3)

式中：E(x)為損失函數(shù)，選用均方誤差；ei為排水?dāng)?shù)據(jù)真實(shí)值與預(yù)測(cè)值之差；2E(x)為損失函數(shù)的Hessian矩陣，即JT(x)J(x)；E(x)為其梯度；N為樣本個(gè)數(shù)。這時(shí)對(duì)于BP神經(jīng)網(wǎng)絡(luò)：

x(i+1)=x(i)-[JT(xi)J(xi)+μiI]-1JT(xi)ei

(4)

即

Δx=-[JT(xi)J(xi)+μiI]-1JT(xi)ei

(5)

式中：μ為大于0的比例系數(shù)；I為單位矩陣；J(x)為損失函數(shù)E(x)的雅可比矩陣，即誤差向量與更新權(quán)重向量的一階微分矩陣。

此時(shí)的Levenberg-Marquardt算法由μ來(lái)控制，μ越大，越接近梯度下降法，當(dāng)μ=0時(shí)，Levenberg-Marquardt變?yōu)楦咚古ｎD法。

2.2 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和模型訓(xùn)練

多層前向神經(jīng)網(wǎng)絡(luò)具有理論上可逼近任意非線性連續(xù)映射的能力，因而非常適合于非線性系統(tǒng)的建模及控制[15]。根據(jù)專家小組提供的數(shù)據(jù)和Fund?o大壩排洪隧洞遭到破壞的實(shí)際情況，影響尾礦排水能力的應(yīng)變因素主要包括庫(kù)水位(Res)、降雨量(Rai)和干灘長(zhǎng)度(Bre)3個(gè)因素[13]。將這3個(gè)因素作為BP網(wǎng)絡(luò)的輸入層神經(jīng)元，排水?dāng)?shù)據(jù)flow作為輸出層的神經(jīng)元。而隱含層神經(jīng)元個(gè)數(shù)則根據(jù)式(6)初步確定[16]：

(6)

式中：m為隱含層節(jié)點(diǎn)數(shù)；n為輸入層節(jié)點(diǎn)數(shù)；l為輸出層節(jié)點(diǎn)數(shù)；α為1～10之間的常數(shù)。網(wǎng)絡(luò)的初步結(jié)構(gòu)設(shè)定為3-3-3-1。

網(wǎng)絡(luò)隱含層的傳遞函數(shù)選用雙曲正切S型傳遞函數(shù)tansig，輸出層則使用線性傳遞函數(shù) purelin。模型的訓(xùn)練選用trainlm，即Levenberg-Marquardt優(yōu)化算法。憑此設(shè)置，使用MATLAB工具建立BP神經(jīng)網(wǎng)絡(luò)模型，經(jīng)

過多次調(diào)整確定網(wǎng)絡(luò)結(jié)構(gòu)為3-4-4-1，網(wǎng)絡(luò)結(jié)構(gòu)如圖1。此時(shí)的回歸系數(shù)是0.985 6，說(shuō)明網(wǎng)絡(luò)的訓(xùn)練效果良好，訓(xùn)練效果如圖2。其中，最大迭代次數(shù)為5 000次，網(wǎng)絡(luò)的學(xué)習(xí)率為0.01，目標(biāo)誤差為0.65×10-3。

圖1 BP雙隱層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Structure of BP double hidden layer neural network

圖2 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練效果Fig.2 BP neural network training effect

2.3 模型預(yù)測(cè)效果驗(yàn)證

為了驗(yàn)證模型的效果好壞，取10組訓(xùn)練樣本之外的數(shù)據(jù)即2015年3月-11月的數(shù)據(jù)輸入模型中進(jìn)行預(yù)測(cè)，其中，只保留庫(kù)水位、降雨量和干灘長(zhǎng)度。預(yù)測(cè)的結(jié)果與實(shí)際真值對(duì)比如圖3、表2所示。

圖3 BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)結(jié)果Fig.3 BP neural network prediction results

序號(hào)日期庫(kù)水位/m降雨量/mm干灘長(zhǎng)度/m樣本真值/(L·s-1)樣本預(yù)測(cè)值/(L·s-1)相對(duì)誤差/%12015-3-1884.988.66209.1370.5169.131.9622015-4-1885.4033.82194.7170.5169.151.9332015-5-1886.9011.67197.1164.6066.162.4142015-6-1887.542.82290.8656.2857.672.4752015-6-5887.9713.17304.0861.7559.653.4062015-6-10888.180.00295.8759.7857.184.3572015-7-1889.462.86168.2665.4764.681.2182015-9-1890.952.90290.8657.8156.622.0792015-10-1891.590.00216.3570.9569.192.48102015-11-1892.8716.26240.3869.8569.081.10

由圖3可知，神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果與實(shí)際結(jié)果之間相差不大，變換趨勢(shì)基本相同，表示該模型訓(xùn)練成功。由表2可知，排水?dāng)?shù)據(jù)flow(L/s)對(duì)應(yīng)的10組驗(yàn)證數(shù)據(jù)的相對(duì)誤差范圍在1.10%至4.35%之間，該BP神經(jīng)網(wǎng)絡(luò)的運(yùn)行時(shí)間為8 s?？梢?，BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型的誤差范圍相對(duì)較小，準(zhǔn)確度較高。

3 基于SVR算法的回歸預(yù)測(cè)模型

3.1 SVR回歸預(yù)測(cè)模型的建立

SVR回歸即找到1個(gè)高維的回歸平面，讓數(shù)據(jù)集內(nèi)的所有樣本到該平面距離最近。SVR也是小樣本數(shù)據(jù)集下性能可觀的預(yù)測(cè)模型。

已知數(shù)據(jù)樣本是三維向量，訓(xùn)練集可表示為{(xi,yi),…,(xk,yk)}，其中xi∈R3是輸入變量即庫(kù)水位、降雨量和干灘長(zhǎng)度，yi∈R是排水?dāng)?shù)據(jù)實(shí)際值，i=1,2,…,k是訓(xùn)練樣本的個(gè)數(shù)。此時(shí)構(gòu)建回歸超平面：

f(x)=ω·x+b

(7)

式中：f(x)為排水?dāng)?shù)據(jù)的預(yù)測(cè)值；ω為權(quán)值向量，ω·x做內(nèi)積運(yùn)算；b為偏移值。

(8)

式中：ε為損失界限，控制了模型對(duì)數(shù)據(jù)的誤差寬度。引入拉格朗日(Lagrange)乘子a,a*將凸二次優(yōu)化問題轉(zhuǎn)變?yōu)閷?duì)偶問題：

(9)

求解參數(shù)，符合決策變量的數(shù)據(jù)成為支持向量回歸，此時(shí)的目標(biāo)函數(shù)為：

(10)

因?yàn)檩斎胩卣髋c輸出特征屬于非線性關(guān)系，本模型采用徑向基核函數(shù)，函數(shù)如式(11)：

(11)

式中：y為核函數(shù)中心；σ為核函數(shù)的寬度參數(shù)， 控制了函數(shù)的徑向作用范圍；x為樣本點(diǎn)。高斯核函數(shù)將數(shù)據(jù)從其特征空間映射到無(wú)窮維空間，解決了非線性不可分的問題，擬合出1條高維的曲線達(dá)到預(yù)測(cè)效果[18]。最后模型的目標(biāo)函數(shù)為：

(12)

模擬過程如圖4所示。

圖4 SVR預(yù)測(cè)模型流程Fig.4 SVR prediction model flow chart

1)回歸預(yù)測(cè)模型是基于特征屬性搭建的，將2011年7月至2015年3月Fund?o尾礦大壩的庫(kù)水位(Res)、降雨量(Rai)和干灘長(zhǎng)度(Bre)作為特征輸入，排水?dāng)?shù)據(jù)作為輸出。

2)采用和BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)排水?dāng)?shù)據(jù)一樣的測(cè)試集。調(diào)用sklearn中的svm庫(kù)，核函數(shù)設(shè)置為rbf徑向基函數(shù)，擬合SVR回歸曲線，采用試湊法經(jīng)過多次試驗(yàn)找到最優(yōu)參數(shù)確定為gamma=26,C=1。

3)完成模型回歸，得到預(yù)測(cè)結(jié)果，與真值對(duì)比計(jì)算誤差。

3.2 結(jié)果驗(yàn)證

借助Anaconda+Pycharm集成開發(fā)環(huán)境，調(diào)用sklearn庫(kù)中的svm模塊，對(duì)數(shù)據(jù)集進(jìn)行支持向量回歸預(yù)測(cè)。目前，常用的最優(yōu)參數(shù)選取方法有K折交叉驗(yàn)證法、留一法和試湊法。留一法和交叉驗(yàn)證法不能固定訓(xùn)練樣本，無(wú)法與BP神經(jīng)網(wǎng)絡(luò)方法做結(jié)果對(duì)比，因此選用試湊法。試湊法的基本原理為：首先為函數(shù)固有的參數(shù)賦初始值，然后開始實(shí)驗(yàn)測(cè)試，根據(jù)測(cè)試精度重復(fù)調(diào)整參數(shù)值，直至得到滿意的測(cè)試精度為止[19]。選取默認(rèn)參數(shù)gamma=1，C=1，固定gamma選用不同C值確定最低均方根誤差(Root Mean Squard Error, RMSE)，此時(shí)gamma=1，C=5，RMSE=3.61。再固定C值，選用不同gamma確定最低RMSE值，此時(shí)gamma=25，C=5，RMSE=2.76。重復(fù)上述步驟直至RMSE值不再降低。經(jīng)過多次調(diào)試，最優(yōu)參數(shù)設(shè)置為gamma=26,C=1，此時(shí)RMSE=2.74。調(diào)試過程見圖5。表3為SVR模型預(yù)測(cè)結(jié)果，可知，預(yù)測(cè)排水?dāng)?shù)據(jù)相對(duì)誤差范圍在1.64%～9.21%之間，最大誤差小于10%，運(yùn)行時(shí)間為0.08 s，所以認(rèn)為該模型的預(yù)測(cè)效果良好。

圖5 SVR模型參數(shù)調(diào)試Fig.5 SVR model parameter modulation

序號(hào)日期樣本真值/(L·s-1)樣本預(yù)測(cè)值/(L·s-1)相對(duì)誤差/%12015-3-170.5171.671.6522015-4-170.5171.681.6532015-5-164.6067.405.0142015-6-156.2859.535.7852015-6-561.7562.921.8862015-6-1059.7865.299.2172015-7-165.4761.496.0882015-9-157.8158.972.0192015-10-170.9572.121.64102015-11-169.8571.011.67

綜合對(duì)比BP神經(jīng)網(wǎng)絡(luò)與SVR回歸預(yù)測(cè)模型，二者的預(yù)測(cè)機(jī)制不同，預(yù)測(cè)效果也不一樣。前者通過反向傳播誤差更新權(quán)值來(lái)預(yù)測(cè)數(shù)據(jù)，后者通過在高維空間擬合超平面使樣本點(diǎn)距離平面最近預(yù)測(cè)數(shù)據(jù)。2個(gè)模型對(duì)比結(jié)果如表4所示，SVR回歸預(yù)測(cè)模型的平均相對(duì)誤差為3.66%，而BP神經(jīng)網(wǎng)絡(luò)模型的平均相對(duì)誤差則只有2.34%。從時(shí)間上看，SVR預(yù)測(cè)模型要更快。

表4 模型對(duì)比結(jié)果Table 4 Results of model comparison

4 結(jié)論

1)以庫(kù)水位、降雨量和干灘長(zhǎng)度為輸入，運(yùn)用BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型和SVR回歸預(yù)測(cè)模型2種方法預(yù)測(cè)排水?dāng)?shù)據(jù)，2種模型的最大誤差分別為4.35%和9.21%，均小于10%，說(shuō)明模型的建立和預(yù)測(cè)結(jié)果是可行的。

2)BP神經(jīng)網(wǎng)絡(luò)和SVR預(yù)測(cè)模型克服了排水?dāng)?shù)據(jù)復(fù)雜且非線性的特點(diǎn)。前者的預(yù)測(cè)精度高，平均誤差為2.34%，后者的運(yùn)行速度快，運(yùn)行時(shí)間為0.08 s。其中，核函數(shù)為徑向基函數(shù)，參數(shù)gamma=26，C=1時(shí)，SVR模型預(yù)測(cè)結(jié)果最優(yōu)。

3)借助BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型和SVR回歸預(yù)測(cè)模型，礦山安全管理部門可根據(jù)模型預(yù)測(cè)結(jié)果預(yù)測(cè)排水?dāng)?shù)據(jù)，對(duì)于尾礦庫(kù)安全管理具有一定參考價(jià)值。