馮先丁,魏鏡弢,吳張永,錢 杰,浦友尚
(昆明理工大學(xué) 機(jī)電工程學(xué)院,云南 昆明 650500)
磨礦作業(yè)是選礦流程中非常關(guān)鍵的環(huán)節(jié)[1],其負(fù)責(zé)最大限度地解離礦石中的有用礦物和脈石礦物。磨礦產(chǎn)品的質(zhì)量直接影響整個(gè)選礦廠的技術(shù)指標(biāo)和經(jīng)濟(jì)效益。磨礦作業(yè)是選礦廠中能源消耗最大的環(huán)節(jié),其生產(chǎn)成本占整個(gè)選礦廠生產(chǎn)總成本的45%~60%,因此準(zhǔn)確有效地預(yù)測(cè)球磨機(jī)內(nèi)部的運(yùn)行狀態(tài)及優(yōu)化磨礦控制技術(shù)是選礦行業(yè)提高選礦效率、提升企業(yè)效益的關(guān)鍵環(huán)節(jié)[2]。
早期的工人根據(jù)球磨機(jī)的聲音和個(gè)人經(jīng)驗(yàn)[3]來判斷球磨機(jī)內(nèi)部的負(fù)荷狀態(tài)。由于球磨機(jī)工作時(shí)影響因素眾多,內(nèi)部環(huán)境復(fù)雜,故這些方法誤差較大。文獻(xiàn)[4~6]中,研究人員通過采集球磨機(jī)外部信號(hào)來預(yù)測(cè)球磨機(jī)負(fù)荷。近年來,磨音信號(hào)成為了研究的熱點(diǎn)。文獻(xiàn)[7]通過采集球磨機(jī)磨音信號(hào),對(duì)起磨過程磨音和停磨過程磨音進(jìn)行了功率譜估計(jì),得到磨音的有效頻段。然后將有效頻段分段后的能量作為特征,采用徑向基(Radial Basis Function,RBF)神經(jīng)網(wǎng)絡(luò)建立預(yù)測(cè)模型來預(yù)測(cè)球磨機(jī)的負(fù)荷。文獻(xiàn)[8]采用主元分析法(Principal Component Analysis,PCA)對(duì)振動(dòng)譜能量值進(jìn)行降維,利用支持向量機(jī)(Support Vector Machine,SVM)建立磨機(jī)負(fù)荷參數(shù)分類模型,實(shí)現(xiàn)了磨機(jī)負(fù)荷的預(yù)測(cè)。該研究證明,球磨機(jī)磨音信號(hào)與球磨機(jī)負(fù)荷有緊密的聯(lián)系。上述方法在對(duì)球磨機(jī)信號(hào)進(jìn)行頻譜特征提取時(shí),由于受到人為分析不精確和考慮不足等因素的影響,導(dǎo)致產(chǎn)生頻段被硬性劃分或無法準(zhǔn)確找到反映球磨機(jī)負(fù)荷特征等問題。此外,上述方法在建立預(yù)測(cè)模型時(shí)沒有考慮 SVM 的最佳參數(shù),導(dǎo)致球磨機(jī)預(yù)測(cè)精度低,模型性能不穩(wěn)定。因此,本文提出一種基于PCA-PSO-SVM的球磨機(jī)負(fù)荷預(yù)測(cè)模型。
球磨機(jī)產(chǎn)生的聲音主要包括磨音和各種干擾噪音[9]。磨音主要分為兩部分:(1)由磨礦介質(zhì)與球磨機(jī)內(nèi)襯板碰撞產(chǎn)生的聲音。這部分聲音主要受到球磨機(jī)的轉(zhuǎn)速和磨礦介質(zhì)的填充量的影響。一般在工程中使用的球磨機(jī)轉(zhuǎn)速較小,所以這部分磨音主要集中在中低頻[10];(2)由筒體內(nèi)磨礦介質(zhì)與磨礦介質(zhì)、磨礦介質(zhì)與物料之間相互撞擊產(chǎn)生磨音。由于磨礦介質(zhì)數(shù)量多、體積小,碰撞機(jī)會(huì)大,所以產(chǎn)生的聲音以中高頻為主。
在磨音中,磨礦介質(zhì)和磨礦介質(zhì)之間、磨礦介質(zhì)和物料之間撞擊產(chǎn)生的聲音主要受筒體內(nèi)的負(fù)荷影響。當(dāng)負(fù)荷較小時(shí),所發(fā)出的磨音主要是磨礦介質(zhì)之間的沖擊和碰撞聲。此時(shí)磨音的頻率高,聲響大,聽起來比較清脆。當(dāng)球磨機(jī)內(nèi)的負(fù)荷較大時(shí),由于筒體內(nèi)磨礦介質(zhì)之間的空隙逐漸被物料填充,磨礦介質(zhì)以及物料間的碰撞和沖擊次數(shù)減少,此時(shí)的磨音頻率低,聲響小,聽起來比較低沉。本文根據(jù)磨音信號(hào)與負(fù)荷之間的聯(lián)系,通過對(duì)磨音信號(hào)進(jìn)行采集、處理和分析,可有效預(yù)測(cè)球磨機(jī)負(fù)荷對(duì)。
Welch算法是對(duì)周期圖法的改進(jìn)算法[11]。將樣本數(shù)據(jù)進(jìn)行分段加窗處理,求得各段功率譜后再求平均功率譜[12]。將信號(hào)x(n)的N個(gè)數(shù)據(jù)點(diǎn)分為L(zhǎng)段,每段長(zhǎng)度為M。
如果對(duì)每一段的數(shù)據(jù)點(diǎn)加矩形窗,第i段的數(shù)據(jù)變成
(1)
由此可以得到平均周期圖法
(2)
這種方法也稱為Bartlett法。
采用Welch法對(duì)信號(hào)數(shù)據(jù)點(diǎn)分段時(shí),對(duì)每段數(shù)據(jù)之間做一個(gè)交疊,并且在對(duì)數(shù)據(jù)段加窗時(shí)可以選擇多種窗。按Bartlett法對(duì)數(shù)據(jù)段求取功率譜,對(duì)所求的結(jié)果進(jìn)行歸一化處理,得到
(3)
式中,U為歸一化因子;d2(n)為數(shù)據(jù)加窗函數(shù)。
相比于周期法,Welch法可有效改善譜曲線的光滑性,并提高譜估計(jì)的分辨能力[13]。
主元分析法是特征提取和數(shù)據(jù)壓縮中的經(jīng)典方法[14-15]。PCA用較少的主元變量反映大部分原始數(shù)據(jù)的信息,可以有效降低特征維數(shù)。PCA的實(shí)現(xiàn)算法步驟如下[16]:
步驟1數(shù)據(jù)標(biāo)準(zhǔn)化處理。對(duì)原始數(shù)據(jù)進(jìn)行z-score標(biāo)準(zhǔn)化處理,使得不同指標(biāo)具有相同的數(shù)量級(jí)和量綱,便于后續(xù)計(jì)算和分析。標(biāo)準(zhǔn)化計(jì)算式為
(4)
式中,yij為矩陣中的第i個(gè)樣本的第j個(gè)觀測(cè)值;yj為矩陣Y中第j列均值;sj為第j列的標(biāo)準(zhǔn)差;
步驟2計(jì)算樣本矩陣的相關(guān)系數(shù)矩陣R。相關(guān)系數(shù)矩陣計(jì)算為
(5)
(6)
式中,m為樣本總數(shù)據(jù);
步驟3求協(xié)方差矩陣的特征根和特征向量。根據(jù)特征方程|λE-R|=0可求出特征根,并按照大小順尋排列為λ1≥λ2≥…≥λn≥0。與特征根相對(duì)應(yīng)的特征向量為α1,α2,…,αn。特征向量組成的矩陣特為特征向量矩陣,如式(7)所示。
(7)
式中,n為樣本數(shù)據(jù)的觀測(cè)量;
步驟4計(jì)算主元矩陣T。主元ti=Yαi表示了原始信息矩陣在某個(gè)主元下對(duì)應(yīng)向量上的投影,它的大小表示了Y在αi方向上的覆蓋程度。主元矩陣T為
T=Ym×nLm×n=t1,t2,…,tn
(8)
式中,Ym×n為樣本原始數(shù)據(jù)矩陣;
步驟5計(jì)算主元貢獻(xiàn)率和累計(jì)貢獻(xiàn)率。主元貢獻(xiàn)率為主元所攜帶的信息占所有信息的百分比。每個(gè)成分變量所攜帶信息可以由其方差來表示,所以每個(gè)主元的方差占總方差S的百分比就是該成分的貢獻(xiàn)率,如式(9)所示。
(9)
前r個(gè)主元累計(jì)貢獻(xiàn)率也就是前r個(gè)特征值的累計(jì)方差值貢獻(xiàn)率,如式(10)所示。
(10)
依據(jù)累計(jì)貢獻(xiàn)率來選取主元,當(dāng)累計(jì)貢獻(xiàn)率達(dá)到 90%或者以上時(shí),前r個(gè)主元就能反映出原有數(shù)據(jù)信息樣本變量包含的大部分信息,即用較少相關(guān)的特征代替了原有的特征。
粒子群優(yōu)化(Particle Swarm Optimization, PSO)算法通過追隨當(dāng)前最優(yōu)值來尋找全局最優(yōu)。PSO 算法具有并行性、分布式處理的特點(diǎn),得到了廣泛的研究和應(yīng)用[17-18]。
PSO算法中粒子群迭代更新的過程為:假設(shè)在一個(gè)D維空間中,一共有N個(gè)粒子,第i個(gè)粒子在D維空間中的飛行速度向量為
vi=(vi1,vi2,vi3,…,viD)
(11)
位置向量為
xi=(xi1,xi2,xi3,…,xiD)
(12)
每次迭代,粒子根據(jù)適應(yīng)度函數(shù)算出相對(duì)應(yīng)的適應(yīng)值,依據(jù)自己和群體的經(jīng)驗(yàn)來更新下一次迭代的位置和速度,將表述轉(zhuǎn)化為直觀的計(jì)算式
(13)
(14)
式中,k為迭代次數(shù);xik為第i個(gè)粒子在第k次迭代中的位置向量;vik為第i個(gè)粒子在第k次迭代中的速度向量;pi是第i個(gè)粒子個(gè)體經(jīng)歷過的最好位置;gi是迭代到目前種群所經(jīng)歷過的最好位置;r1和r2為隨機(jī)函數(shù);c1和c2為認(rèn)知因子和社會(huì)學(xué)習(xí)因子;ω為慣性權(quán)重。
支持向量機(jī)是建立在結(jié)構(gòu)風(fēng)險(xiǎn)最小原則和統(tǒng)計(jì)學(xué)理論基礎(chǔ)上的一種機(jī)器學(xué)習(xí)方法。SVM最早被用于解決分類問題,后又將其應(yīng)用于解決回歸分析和函數(shù)擬合問題[19]。最優(yōu)分類函數(shù)表示為
(15)
式中,λi為拉格朗日乘子;b為偏置;xi為數(shù)據(jù)點(diǎn);yi為類別。
將上式用于處理非線性樣本時(shí),其無法準(zhǔn)確劃分樣本,因此引入了懲罰函數(shù)C、不敏感損失參數(shù)ε和核參數(shù)g。懲罰因子主要表示模型對(duì)誤差的容忍度。不敏感損失參數(shù)主要影響學(xué)習(xí)機(jī)的擬合效果,其取值越小,模型的復(fù)雜程度越大,擬合的精度也就越高。核函數(shù)決定了數(shù)據(jù)從原始空間映射到更高維的特征空間的分布情況,從而影響模型的泛化性能。
搭建球磨機(jī)磨音采集系統(tǒng),采集10個(gè)負(fù)荷下球磨機(jī)的磨音信號(hào),每組實(shí)驗(yàn)得到40個(gè)樣本,共400個(gè)樣本。采用PCA進(jìn)行特征提取,得到每個(gè)樣本的21維特征,進(jìn)而得到400×21的樣本集。將樣本集隨機(jī)分為訓(xùn)練集(樣本個(gè)數(shù)為300)和測(cè)試集(樣本個(gè)數(shù)為100)。為了便于分析并提高計(jì)算效率,在進(jìn)行模型建立之前需對(duì)樣本數(shù)據(jù)進(jìn)行歸一化處理,將原始數(shù)據(jù)規(guī)范到[-1,1]之間,表達(dá)式為
(16)
式中,x為需要被歸一化的數(shù)據(jù);y為歸一化后的結(jié)果。部分?jǐn)?shù)據(jù)如表1所示。
表1 歸一化后的樣本集(部分)
對(duì) SVM 中的參數(shù)進(jìn)行初始化設(shè)置:懲罰因子C的取值范圍為(0.1,1 000),核參數(shù)g的取值范圍為(0.01,100),不敏感損失系數(shù)ε的取值范圍為(0,1)。PSO 中的參數(shù)設(shè)定為:種群粒子個(gè)數(shù)為 20,慣性權(quán)重w=1,學(xué)習(xí)因子c1=1.5,c2=1.7。將交叉驗(yàn)證法(K-CV)的均方誤差平均值(Mean Square Error,MSE)作為適應(yīng)度函數(shù),以最大迭代次數(shù) 200 作為終止條件,得到支持向量機(jī)的3個(gè)參數(shù)最優(yōu)取值。
本文提出的基于PCA-PSO-SVM的球磨機(jī)負(fù)荷預(yù)測(cè)模型整體步驟如下:
步驟1利用實(shí)驗(yàn)室搭建的磨音信號(hào)采集系統(tǒng)采集球磨機(jī)正常工作時(shí)不同負(fù)荷下的磨音信號(hào),并對(duì)時(shí)域信號(hào)進(jìn)行樣本分段;
步驟2利用Welch法計(jì)算磨音信號(hào)的功率譜,進(jìn)行頻譜分析,并進(jìn)行分段,得到基本特征;
步驟3利用PCA主元分析法對(duì)基本特征進(jìn)行特征提取,降低特征的維度和特征之間的冗余度;
步驟4利用PSO算法對(duì)SVM的3個(gè)參數(shù)進(jìn)行尋優(yōu),得到最優(yōu)參數(shù);
步驟5利用尋優(yōu)后得到的參數(shù)建立SVM預(yù)測(cè)模型,用訓(xùn)練集數(shù)據(jù)訓(xùn)練模型,最終對(duì)測(cè)試集進(jìn)行球磨機(jī)負(fù)荷預(yù)測(cè)。
采用主元分析法對(duì)球磨機(jī)磨音信號(hào)功率譜進(jìn)行特征提取,將特征值從大到小進(jìn)行排列,計(jì)算每個(gè)特征值對(duì)應(yīng)的主元貢獻(xiàn)率和累計(jì)貢獻(xiàn)率。主元累計(jì)貢獻(xiàn)率的變化如圖1所示。
為了確保特征提取后的特征能夠反映原始數(shù)據(jù)的絕大部分信息,本文將累計(jì)貢獻(xiàn)率96%作為確定主元個(gè)數(shù)的指標(biāo),因此本文將主元個(gè)數(shù)確定為 21。
圖1 主元累計(jì)貢獻(xiàn)率變化圖 Figure 1. Diagram of principal cumulative contribution rate change
通過 150次進(jìn)化迭代計(jì)算,得到支持向量機(jī)的最優(yōu)懲罰因子C=89.751 2,核參數(shù)g=0.109 9,不敏感損失系數(shù)ε=0.422 35。得到的最佳適應(yīng)度變化進(jìn)行曲線如圖2所示。
圖2 最佳適應(yīng)度進(jìn)化曲線Figure 2. Evolution curve of best fitness
由圖可知,適應(yīng)度隨進(jìn)化代數(shù)的增加而呈階梯式下降,說明尋優(yōu)算法沒有陷入局部最優(yōu)值。當(dāng)?shù)螖?shù)到 22 次時(shí)適應(yīng)度不再改變,得到最優(yōu)參數(shù)值。
本文利用優(yōu)化后的參數(shù)建立PCA-PSO-SVM預(yù)測(cè)模型。用訓(xùn)練集數(shù)據(jù)對(duì)模型進(jìn)行訓(xùn)練,之后分別對(duì)訓(xùn)練集和測(cè)試集進(jìn)行預(yù)測(cè),得到預(yù)測(cè)結(jié)果。測(cè)試集預(yù)測(cè)結(jié)果誤差如圖3所示。
圖3 PCA-PSO-SVM 測(cè)試集預(yù)測(cè)結(jié)果誤差圖
為了更好地證明本文提出PCA-PSO-SVM預(yù)測(cè)模型的精確性和穩(wěn)定性,將這個(gè)模型的預(yù)測(cè)結(jié)果分別與 SVM 和 PCA-SVM預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析。圖4、圖5分別為SVM 和 PCA-SVM的預(yù)測(cè)誤差圖。
為了更好地說明新預(yù)測(cè)模型的精確性和穩(wěn)定性,本文引入3個(gè)評(píng)價(jià)指標(biāo):均方根誤差(Root Mean Square Error,RMSE)、平均絕對(duì)誤差(Mean Absolute Error,MAE)以及平均百分比誤差(Mean Average Percentage Error,MAPE)。對(duì)比結(jié)果如表2所示。
圖4 SVM測(cè)試集預(yù)測(cè)結(jié)果誤差圖Figure 4. Error graph of SVM test set prediction result
圖5 PCA-SVM測(cè)試集預(yù)測(cè)結(jié)果誤差圖Figure 5. Error graph of PCA-SVM test set prediction result
表2 3種預(yù)測(cè)結(jié)果評(píng)價(jià)指標(biāo)
由表中信息可知,PCA-PSO-SVM負(fù)荷預(yù)測(cè)模型的各項(xiàng)指標(biāo)都優(yōu)于PCA-SVM和SVM,且PCA-SVM模型的評(píng)級(jí)指標(biāo)也高于SVM。PCA-PSO-SVM負(fù)荷預(yù)測(cè)模型不僅降低了預(yù)測(cè)的誤差,也提高了預(yù)測(cè)速度。
本文通過對(duì)樣本原始數(shù)據(jù)進(jìn)行Welch法功率譜分析,并利用PCA對(duì)數(shù)據(jù)進(jìn)行特征提取和降維處理,不僅降低了特征之間的冗余度和特征維度,還減少了樣本的干擾信息,簡(jiǎn)化了樣本數(shù)據(jù)。本文通過PSO對(duì)SVM的相關(guān)參數(shù)進(jìn)行尋優(yōu),得到建立該模型所需的最優(yōu)參數(shù),提高了預(yù)測(cè)模型的推廣能力,提升了預(yù)測(cè)模型的預(yù)測(cè)精度和穩(wěn)定性,降低了預(yù)測(cè)誤差。本文提出的PCA-PSO-SVM負(fù)荷預(yù)測(cè)模型通過球磨機(jī)磨音信號(hào)準(zhǔn)確地預(yù)測(cè)其內(nèi)部負(fù)載,誤差范圍小,算法運(yùn)算時(shí)間較短,計(jì)算效率有所提高。
由于本文的研究對(duì)象是實(shí)驗(yàn)室小型球磨機(jī),故實(shí)驗(yàn)采用間斷式給料磨礦,這與實(shí)際現(xiàn)場(chǎng)連續(xù)給料的球磨機(jī)采集的信號(hào)數(shù)據(jù)仍有差異。在今后的實(shí)際應(yīng)用中,需要考慮到這一點(diǎn)。