楊麗榮 江川 黎嘉駿 曹沖 周俊

(1 江西理工大學(xué)機(jī)電工程學(xué)院 贛州 341000)

(2 江西省礦冶機(jī)電工程研究中心 贛州 341000)

0 引言

在外界各種因素作用下，巖石內(nèi)部會(huì)產(chǎn)生應(yīng)力集中現(xiàn)象，但是這種應(yīng)力集中的狀態(tài)是不穩(wěn)定的，因此應(yīng)力會(huì)向穩(wěn)定的低能狀態(tài)進(jìn)行過(guò)度，在過(guò)程中應(yīng)變能會(huì)采用彈性波的形式向外擴(kuò)散，即聲發(fā)射現(xiàn)象[1-2]。聲發(fā)射技術(shù)[3]可以預(yù)測(cè)被測(cè)對(duì)象的內(nèi)部變化，是當(dāng)前最佳的礦井安全檢測(cè)技術(shù)[4]，且可以適用于巖石失穩(wěn)的預(yù)測(cè)，提前發(fā)現(xiàn)危險(xiǎn)狀況。在實(shí)際監(jiān)測(cè)過(guò)程中，能夠準(zhǔn)確獲知挖礦過(guò)程中巖石的破裂狀態(tài)，就能更快更多減少礦井安全事故發(fā)生[5-6]。巖石破裂狀態(tài)的獲取，需要對(duì)巖石破裂過(guò)程中的歷史數(shù)據(jù)進(jìn)行準(zhǔn)確的特征向量提取。為了保留特征信號(hào)中的有用信息并對(duì)大量的數(shù)據(jù)特征進(jìn)行簡(jiǎn)約，需要對(duì)特征數(shù)據(jù)進(jìn)行降維處理。因此融合特征下數(shù)據(jù)降維算法的好壞對(duì)于巖石破裂狀態(tài)分類及失穩(wěn)預(yù)測(cè)具有重要意義。

針對(duì)原始數(shù)據(jù)特征存在大量的冗余信息、數(shù)據(jù)的數(shù)維偏高及區(qū)分性不強(qiáng)等問(wèn)題，國(guó)內(nèi)外學(xué)者對(duì)此提出了不同的解決方法。Ibrhim等[7]抽取聲發(fā)射信號(hào)與力信號(hào)的各種統(tǒng)計(jì)參數(shù)，然后利用主元分析法(Principal component analysis,PCA)剔除冗余或與刀具磨損相關(guān)性小的信息，發(fā)現(xiàn)聲發(fā)射信號(hào)的均值、標(biāo)準(zhǔn)偏差和力信號(hào)在三個(gè)方向分量的最大值與刀具磨損相關(guān)性最大，達(dá)到降維的目的。劉開(kāi)南等[8]提出了一種用于圖像分類的Mod-局部線性嵌入(Locally linear embedding,LLE)算法，通過(guò)對(duì)圖像數(shù)據(jù)集進(jìn)行測(cè)試，結(jié)果表明圖像分類的精度較高且降維性能良好。Cheng等[9]提出了一種增量LLE算法，通過(guò)對(duì)衛(wèi)星的遙測(cè)數(shù)據(jù)進(jìn)行特征降維來(lái)實(shí)現(xiàn)故障檢測(cè)，其仿真的結(jié)果表明該方案的有效性。劉杏芳等[10]通過(guò)對(duì)3D 地震數(shù)據(jù)的層問(wèn)屬性進(jìn)行特征提取，比較LLE 和PCA 兩種降維方法的聚類效果，實(shí)例應(yīng)用表明LLE 算法較好地保持了原始數(shù)據(jù)的結(jié)構(gòu)特征，具有更高的特征提取能力。王玥[11]使用LLE、PCA、多維縮放(Multiple dimensional scaling,MDS)幾種降維算法對(duì)股票的收盤(pán)價(jià)數(shù)據(jù)進(jìn)行特征簡(jiǎn)約再預(yù)測(cè)分析，并和幾種直接回歸預(yù)測(cè)方法進(jìn)行對(duì)比，結(jié)果表明經(jīng)過(guò)LLE 算法降維后的預(yù)測(cè)精度提高顯著，更加精確。

通過(guò)以上分析，比較各種特征數(shù)據(jù)降維方法的優(yōu)劣，發(fā)現(xiàn)LLE 和PCA 算法作為主流的數(shù)據(jù)處理方法已經(jīng)廣泛應(yīng)用于科研的方方面面。通過(guò)對(duì)高維的數(shù)據(jù)進(jìn)行特征簡(jiǎn)約，保留了原始數(shù)據(jù)的重要特征，使得簡(jiǎn)單模型在小數(shù)據(jù)集上有更強(qiáng)的魯棒性。本文通過(guò)PCA與LLE算法的降維融合特征實(shí)驗(yàn)對(duì)比，分析4 個(gè)狀態(tài)下融合特征的聚類效果，選取聚類效果更集中的算法作為巖石破裂信號(hào)特征數(shù)據(jù)的降維方法，為分析巖體失穩(wěn)狀態(tài)提供一定的理論依據(jù)和基礎(chǔ)。

1 LLE算法原理及特點(diǎn)分析

針對(duì)PCA忽略原始信號(hào)的局部信息的問(wèn)題，一種無(wú)監(jiān)督的非線性局部嵌入算法被提出，LLE 算法的主要思想是對(duì)局部鄰域線性化處理，通過(guò)將原始數(shù)據(jù)從高維映射到低維特征空間，從而使數(shù)據(jù)可視化[12]。該算法通過(guò)不同鄰域的不同數(shù)據(jù)點(diǎn)的線性重構(gòu)權(quán)值，得到局部鄰域中的每個(gè)數(shù)據(jù)點(diǎn)，再將其映射到低維空間，同時(shí)要保證局部鄰域與重構(gòu)權(quán)值的關(guān)系，這樣得到的低維嵌入空間的樣本點(diǎn)線性重構(gòu)權(quán)值和高維空間中表示的線性權(quán)值相同[13]。計(jì)算過(guò)程為以下3個(gè)步驟：

(1) 近鄰域選擇

若想要得到高維數(shù)據(jù)集X=[x1,x2,···,xn]∈RD×N各數(shù)據(jù)點(diǎn)的最近鄰點(diǎn)，可以利用各數(shù)據(jù)點(diǎn)之間的歐氏距離來(lái)尋找，其中，點(diǎn)xi的k個(gè)最近鄰點(diǎn)集合為R(xi)={xi1,xi2,···,xik}。

(2) 計(jì)算數(shù)據(jù)點(diǎn)局部最優(yōu)近鄰線性重構(gòu)權(quán)值

對(duì)數(shù)據(jù)集X=[x1,x2,···,xn]∈RD×N中的所有點(diǎn)，給予約束條件：若是xi的近鄰點(diǎn)，即xj ∈R(xi)，則線性重構(gòu)權(quán)值表示為wij；若xj不是xi的近鄰點(diǎn)，即xj/∈R(xi)，則wij=0。因此，可用最小二乘重構(gòu)誤差極小化原則來(lái)處理最佳線性重構(gòu)權(quán)值矩陣：

數(shù)據(jù)點(diǎn)xj對(duì)xi在重構(gòu)過(guò)程中的貢獻(xiàn)為權(quán)值wij。各數(shù)據(jù)點(diǎn)與其近鄰點(diǎn)之間的局部幾何關(guān)系為重構(gòu)權(quán)值矩陣W={wij}。由于數(shù)據(jù)點(diǎn)的近鄰線性重構(gòu)權(quán)值不會(huì)發(fā)生改變，因此當(dāng)數(shù)據(jù)集發(fā)生各種變化時(shí)，重構(gòu)權(quán)值矩陣都不會(huì)改變，因此，式(1)可以表達(dá)成如下形式：

(3) 求解原始信號(hào)在低維空間下的坐標(biāo)

低維嵌入空間下的坐標(biāo)及其近鄰點(diǎn)的計(jì)算，必須保持高維數(shù)據(jù)空間點(diǎn)之間的重構(gòu)權(quán)值固定不變，利用極小化誤差函數(shù)求低維嵌入坐標(biāo)。即：

式(3)中，yi為低維嵌入坐標(biāo)；Y={yi}為低維嵌入坐標(biāo)矩陣。低維坐標(biāo)矩陣在進(jìn)行旋轉(zhuǎn)、平移要想保持不變，必須滿足條件：

則計(jì)算低維嵌入坐標(biāo)的極小化誤差函數(shù)，式(3)可以用如下形式表示：

同理，將式(6)進(jìn)行推導(dǎo)，把低維嵌入坐標(biāo)Y求解轉(zhuǎn)化為最小特征值的求解問(wèn)題，即

式(6)～(7)中，M=(I -W)T(I -W)為n×n的矩陣；矩陣元素I表示單位矩陣；λ表示矩陣的特征值。M的最小d+1個(gè)特征值對(duì)應(yīng)的特征向量[γ2,γ3,···,γd+1]為低維嵌入坐標(biāo)Y，可表示為

根據(jù)分析LLE 算法計(jì)算步驟可知，由于低維嵌入坐標(biāo)的準(zhǔn)確性會(huì)隨著k值發(fā)生變化，因此該算法的前提是需要設(shè)定每個(gè)樣本點(diǎn)的近鄰參數(shù)k，同時(shí)，近鄰點(diǎn)數(shù)過(guò)多會(huì)導(dǎo)致原始數(shù)據(jù)中的非近鄰點(diǎn)也被納入近鄰區(qū)域；近鄰點(diǎn)個(gè)數(shù)過(guò)少會(huì)使高維原始數(shù)據(jù)集中的內(nèi)在結(jié)構(gòu)發(fā)生扭曲且近鄰域不連通。

2 砂巖單軸壓縮試驗(yàn)設(shè)計(jì)

2.1 巖樣選擇與制備

試驗(yàn)中的紅砂巖采自兩地礦山以確保實(shí)驗(yàn)結(jié)果的泛化能力。紅砂巖G采用贛州某地質(zhì)脆性粗粒紅砂巖，紅砂巖R 取自廣西某錫礦。為達(dá)到本次試驗(yàn)的預(yù)期目標(biāo)，將取回的巖體由鉆孔取樣機(jī)得到柱狀巖芯，經(jīng)巖體自動(dòng)切割機(jī)進(jìn)行分割，最后得到標(biāo)準(zhǔn)巖樣50 mm×100 mm為兩個(gè)端面被打磨后的巖樣，具體如表1和圖1所示。

圖1 不同紅砂巖試驗(yàn)巖樣Fig.1 Different red sandstone test samples

表1 試驗(yàn)巖樣規(guī)格Tabel 1 Test rock sample specifications

2.2 試驗(yàn)裝置與參數(shù)設(shè)置

試驗(yàn)采用RMT-150C 型加載系統(tǒng)，該加載系統(tǒng)在剛度、靜態(tài)特性和動(dòng)態(tài)特性方面具有優(yōu)勢(shì)。伺服控制系統(tǒng)由液壓動(dòng)力源提供穩(wěn)定輸出，可通過(guò)行程、位移、荷載以及不同組合形式、不同加載速率的控制方式，實(shí)現(xiàn)單軸、三軸、間接拉伸及剪切試驗(yàn)等多種試驗(yàn)?？刂颇Ｊ讲捎梦灰萍虞d，加載速率0.002 mm/s，力終點(diǎn)200 kN，力極限250 kN，位移終點(diǎn)2 mm，位移極限2.5 mm。本次實(shí)驗(yàn)加載如圖2所示。

圖2 RMT-150C 型巖石力學(xué)試驗(yàn)系統(tǒng)Fig.2 RMT-150C rock mechanics test system

多傳感器集合成的應(yīng)力應(yīng)變采集系統(tǒng)包含測(cè)量垂直液壓缸活塞桿的行程傳感器、兩種測(cè)量垂直液壓缸輸出力的力傳感器、測(cè)量試件軸向變形軸向位移傳感器。聲發(fā)射傳感器放置在試件中部軸線位置，試件和傳感器接觸部位采用真空脂作為耦合劑，最后采用透明膠布進(jìn)行固定，傳感器布局如圖3 所示。1000 kN力傳感器置于液壓缸內(nèi)；軸向位移傳感器放置在傳感器夾持器上，與上壓頭表面相互接觸。各傳感器量程范圍赫爾精度如表2所示。

圖3 聲發(fā)射傳感器布局Fig.3 Acoustic emission sensor layout

表2 應(yīng)力應(yīng)變采集傳感器性能參數(shù)Tabel 2 Stress and strain acquisition sensor performance parameters

2.3 紅砂巖破裂全過(guò)程的應(yīng)力-應(yīng)變特征分析

紅砂巖失穩(wěn)過(guò)程中垂直力和垂直形變隨時(shí)間變化的數(shù)據(jù)保存在RMT-150C 型巖石力學(xué)試驗(yàn)系統(tǒng)中，通過(guò)式(9)可計(jì)算紅砂巖破裂失穩(wěn)過(guò)程的應(yīng)力值和應(yīng)變值隨時(shí)間變化的關(guān)系：

式(9)中，σ為應(yīng)力，MPa；F為垂直力，kN；S為橫截面積，mm2。

根據(jù)式(9)繪制單軸壓縮試驗(yàn)下紅砂巖R10 和G15的應(yīng)力-應(yīng)變曲線如圖4所示。

圖4 單軸壓縮下不同紅砂巖應(yīng)力-應(yīng)變曲線Fig.4 Stress-strain curves of different red sandstones under uniaxial compression

由圖4(a)可以看出，R10 巖樣在應(yīng)力前期，應(yīng)變變化率減小，巖石內(nèi)部相對(duì)變形率在應(yīng)力達(dá)到12.1 MPa時(shí)逐漸保持不變，試件被破壞產(chǎn)生大量新的微裂紋匯合貫通的應(yīng)力為45.7 MPa，軸向應(yīng)變率增大。

從圖4(b)可知，在預(yù)加載時(shí)，G15 巖樣在未發(fā)生應(yīng)變時(shí)就有一定的軸應(yīng)力，此后先向上凹，斜率由小變大，在應(yīng)變?yōu)?.0014、斜率為直線時(shí)應(yīng)力為5.9 MPa，試件所能承受的極限為28.2 MPa，破壞后蠕變過(guò)程產(chǎn)生可能是加載應(yīng)力和加載的持續(xù)時(shí)間導(dǎo)致的。

綜上由R10 和G15 的應(yīng)力-應(yīng)變曲線可知，單軸破壞過(guò)程表現(xiàn)為4 個(gè)階段。壓密接階段(OA)：巖體內(nèi)部存在微小裂紋、空隙，等初始損傷受壓閉合，曲線呈現(xiàn)上凹形狀，斜率逐漸增大。彈性變形階段(AB)：巖體由不連續(xù)狀態(tài)進(jìn)入連續(xù)狀態(tài)，曲線呈現(xiàn)線性增長(zhǎng)趨勢(shì)，斜率保持不變。塑性變形階段(BC)：曲線斜率逐漸減小到零。失穩(wěn)破壞階段(CD)：巖體強(qiáng)度大大降低，宏觀破裂基本形成，但破裂面存在一定摩擦，且具備一定承載能力。

2.4 砂巖聲發(fā)射信號(hào)時(shí)域特征提取

對(duì)于時(shí)間序列模型，目前存在的主要問(wèn)題是僅對(duì)平穩(wěn)過(guò)程具有較好的分析效果，對(duì)非平穩(wěn)過(guò)程則表現(xiàn)得無(wú)能為力。而巖石破裂過(guò)程中產(chǎn)生的聲發(fā)射信號(hào)表現(xiàn)為較強(qiáng)的非平穩(wěn)過(guò)程，直接進(jìn)行AR 建模效果不太理想，而經(jīng)過(guò)改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble empirical mode decomposition,EEMD)正好解決了這一問(wèn)題。因此本文在建立AR模型之前先采用改進(jìn)EEMD 對(duì)巖石聲發(fā)射信號(hào)進(jìn)行預(yù)處理，將其分解為若干個(gè)有效的平穩(wěn)的本征模態(tài)函數(shù)(Intrinsic mode function,IMF)分量，然后對(duì)各IMF 分量采用AIC 準(zhǔn)則函數(shù)確定模型的最佳階次，建立各個(gè)IMF 分量的AR模型，提取AR模型的自回歸參數(shù)構(gòu)成特征向量。

采用改進(jìn)的EEMD 算法對(duì)第一破裂階段下連續(xù)聲發(fā)射時(shí)域波形進(jìn)行分解，取前8 個(gè)IMF 波形如圖5 所示?？梢钥闯觯琁MF5～I(xiàn)MF8分量是原始信號(hào)中分解出的時(shí)間尺度最長(zhǎng)、頻率較低的分量，代表信號(hào)中的低頻成分。

圖5 改進(jìn)的EEMD 算法分解后8 個(gè)IMF 波形Fig.5 Eight IMF waveforms decomposed by the improved EEMD algorithm

為了更加清晰地了解閾值選取的普遍性，采用改進(jìn)的EEMD 算法分別對(duì)10 組不同巖石破裂狀態(tài)下聲發(fā)射信號(hào)進(jìn)行分解，并求取分解出的各分量與原始信號(hào)的云相似值，最后取每組的前6 個(gè)分量的云相似值進(jìn)行曲線擬合，如圖6所示。

圖6 各樣本前6 個(gè)分量的云相似值擬合曲線圖Fig.6 The fitting curve of cloud similarity value of the first 6 components of each sample

從圖6 可知，前3 個(gè)分量(IMF1～MF3)的云相似值都很大。由云模型理論可知，云相似值越大則其含有原始信號(hào)的信息量越多，與原始信號(hào)的關(guān)聯(lián)性更強(qiáng)，可以完全反映原始信號(hào)的特征。而IMF4～I(xiàn)MF6 的值相對(duì)較小，說(shuō)明是與原始信號(hào)無(wú)關(guān)的噪聲。根據(jù)多組信號(hào)統(tǒng)計(jì)分析，將IMF1、IMF2、IMF3確定為有效分量。

為了消除其它因素對(duì)模型的影響，建模前采用式(10)對(duì)各IMF分量進(jìn)行能量歸一化處理：

式(10)中，ci(t)(i=1,2,3)為改進(jìn)的EEMD處理后得到的前3 個(gè)IMF 分量；?ci(t) 為能量歸一化后的IMF分量。

對(duì)砂巖同一破裂及不同破裂階段信號(hào)經(jīng)改進(jìn)的EEMD分解的各IMF 分量進(jìn)行AR建模，得到各分量的準(zhǔn)則函數(shù)AIC值隨模型階次的變化。對(duì)各階IMF分量，模型階次n ＞5 時(shí)，AIC值變化很小而且小于n為1～5的AIC值，因此AR模型的最佳階次為5，可建立5 階AR 模型，提取各IMF 分量模型系數(shù)，構(gòu)造15維特征向量，如式(11)：

式(11)中，φmn表示信號(hào)經(jīng)改進(jìn)的EEMD分解后第m個(gè)IMF分量的第n階模型系數(shù)。

2.5 砂巖聲發(fā)射信號(hào)頻域特征提取

巖石破碎過(guò)程中，一些不確定因素以及巖石材料密度不均勻，導(dǎo)致得到的聲發(fā)射信號(hào)具有較強(qiáng)的非線性特征。在實(shí)際工程應(yīng)用中，不僅需要信號(hào)的幅度信息，還需要信號(hào)的相位信息，雙譜不僅能夠檢測(cè)到信號(hào)的幅值，還可以得到信號(hào)的相位信息。均值為零的高斯過(guò)程，其三階累積量、雙譜均為零，因此雙譜分析法能很好地消除高斯噪聲。基于此，將砂巖不同破裂階段下的信號(hào)進(jìn)行雙譜分析，為了消除其他條件變化對(duì)采樣信號(hào)的影響，在應(yīng)用雙譜分析之前先對(duì)采樣數(shù)據(jù)進(jìn)行去均值及歸一化處理，如式(12)所示：

分別對(duì)巖石不同破裂階段下聲發(fā)射信號(hào)進(jìn)行雙譜分析，結(jié)果如圖7 所示。由圖7 可知，隨著軸應(yīng)力的增加，聲發(fā)射信號(hào)逐漸向低頻擴(kuò)散聚集，相比于第一破裂階段，第二階段在低頻逐漸聚集成塊，頻帶范圍也較寬。在第三破壞階段下，頻帶范圍最寬，低頻高頻都有，且高頻信息聚集；在巖樣破壞后第四階段，應(yīng)力值逐漸減小，但仍具有聲發(fā)射現(xiàn)象，但強(qiáng)度較弱，高頻信號(hào)消失，頻率主要又轉(zhuǎn)向低頻聚集。

圖7 巖石不同破裂階段下聲發(fā)射信號(hào)的雙譜圖Fig.7 Bispectral images of acoustic emission signals at different fracture stages of rocks

也可直觀地看出，信號(hào)經(jīng)去均值及歸一化處理后，砂巖聲發(fā)射信號(hào)的雙譜特征可以較好地區(qū)分巖石不同破裂階段，但最終是得到雙譜高階矩陣，數(shù)據(jù)量的維數(shù)較大，不適于機(jī)器識(shí)別，因此應(yīng)用奇異值分解理論對(duì)高階矩陣進(jìn)行處理，提取出多個(gè)奇異值作為最具代表性的特征向量，以便進(jìn)行智能識(shí)別。

在矩陣的奇異值中，根據(jù)奇異值的大小可判斷信號(hào)的種類，較大奇異值反映的是信號(hào)中的主要特征，剩余較小的部分為噪聲信號(hào)。由此定義奇異譜：

式(13)中，pi表示各個(gè)狀態(tài)變量在整個(gè)系統(tǒng)中所占能量的相對(duì)關(guān)系，也稱為由矩陣A經(jīng)奇異值分解得到的奇異譜。通常篩選前s個(gè)奇異譜累計(jì)貢獻(xiàn)率大于85%的分量作為特征向量并降維，即：

經(jīng)過(guò)大量數(shù)據(jù)處理發(fā)現(xiàn)，砂巖出于不同破裂狀態(tài)時(shí)，其累積貢獻(xiàn)率大于85%的奇異譜數(shù)量是不同的，但都分布在4～10 之間。因此，為了使構(gòu)造的特征向量維數(shù)一致，統(tǒng)一選取前10個(gè)奇異譜構(gòu)成特征向量Tho∈R1×10，

式(15)中，pi為由矩陣A經(jīng)奇異值分解得到的奇異譜。

為了說(shuō)明提取砂巖破裂聲發(fā)射信號(hào)奇異值特征的聚類效果，將奇異譜p1、p2和p1、p3的二維分布繪成散度圖，如圖8 和圖9 所示；結(jié)合兩圖分析可以看出，采用雙譜奇異值分解提取的砂巖破裂各階段特征的聚類效果明顯，可以較好地區(qū)分開(kāi)4 個(gè)階段，為后續(xù)分類與識(shí)別提供依據(jù)。

圖8 奇異譜p1和p2 的散度圖Fig.8 Divergence plots for singular spectra p1 and p2

圖9 奇異譜p1和p3 散度圖Fig.9 Divergence plots for singular spectra p1 and p3

3 融合特征向量提取

3.1 構(gòu)造聯(lián)合多特征向量

(1) 采集紅砂巖破裂全過(guò)程的實(shí)時(shí)數(shù)據(jù)，根據(jù)砂巖不同破裂階段的特點(diǎn)來(lái)劃分所采集到的聲發(fā)射信號(hào)，同時(shí)讀取該聲發(fā)射信號(hào)所對(duì)應(yīng)的應(yīng)力值，構(gòu)建巖石破裂狀態(tài)和聲發(fā)射信號(hào)之間相互對(duì)應(yīng)的關(guān)系。

(2) 將實(shí)驗(yàn)選取的各階段下的75組信號(hào)進(jìn)行特征提取。①采用EEMD 和云相似度結(jié)合的巖石聲發(fā)射信號(hào)預(yù)處理方法，由時(shí)域特征提取各IMF 分量的AR模型系數(shù)，構(gòu)造15維特征向量Tar。②將砂巖不同破裂階段下的信號(hào)進(jìn)行雙譜分析，并對(duì)雙譜高階矩陣進(jìn)行奇異值分解，將提取的多個(gè)奇異譜值作為砂巖聲發(fā)射信號(hào)頻域特征向量，構(gòu)造10維奇異譜特征向量Tho。

(3)將步驟(2)中 ①、②小步構(gòu)造的特征向量進(jìn)行首尾相連，構(gòu)造如式(16)所示的聯(lián)合多特征向量：

聯(lián)合多特征向量T為同一個(gè)樣本信號(hào)，采用不同的信號(hào)分析方法提取多個(gè)域的特征，組合為一個(gè)25維特征向量。

(4) 選取相同或不同階段下的樣本數(shù)據(jù)，重復(fù)進(jìn)行步驟(2)、步驟(3)的處理，提取出對(duì)應(yīng)應(yīng)力此時(shí)聲發(fā)射信號(hào)的聯(lián)合多特征向量。

3.2 融合特征提取對(duì)比實(shí)驗(yàn)驗(yàn)證

3.2.1 PCA算法的降維特征融合實(shí)驗(yàn)

根據(jù)紅砂巖破裂狀態(tài)分類(表3)，在每個(gè)破裂階段下，間隔均勻地選取75 組不同的數(shù)據(jù)，例如第一階段，其應(yīng)力值范圍為：0～12.1 MPa，提取的樣本數(shù)據(jù)對(duì)應(yīng)的應(yīng)力值應(yīng)該盡量均勻地分布在0～12.1 MPa 下，那么將總共選取的300 組樣本采用上述方法提取聯(lián)合多特征向量，構(gòu)成一個(gè)300×25的特征矩陣M，每個(gè)樣本的聯(lián)合多特征表示為矩陣的每一行。

表3 R10 紅砂巖破裂臨界點(diǎn)應(yīng)力值Tabel 3 R10 red sandstone fracture critical point stress value

采用PCA 對(duì)矩陣M進(jìn)行特征降維，得到主元系數(shù)矩陣Coeff，即

在Coeff 矩陣中，第i列表示第i個(gè)主元成分的系數(shù)向量。此外，PCA 同時(shí)還獲得了樣本協(xié)方差矩陣latent，它是由15 個(gè)按降序排列的特征值構(gòu)成的列向量，它們代表了每個(gè)主元對(duì)整體的貢獻(xiàn)。所選取的主元個(gè)數(shù)按照累積貢獻(xiàn)率大于85%的準(zhǔn)則來(lái)確定。圖10 為各主元的貢獻(xiàn)率及累積貢獻(xiàn)率的變化趨勢(shì)。

圖10 各主元貢獻(xiàn)率及累積貢獻(xiàn)率Fig.10 Contribution rate and cumulative contribution rate of each main element

圖10 中紅色柱狀圖代表貢獻(xiàn)率，藍(lán)色的曲線代表累積貢獻(xiàn)率。從圖10 中可以看出：前3 個(gè)特征值相對(duì)較大，最大為0.126；后4 個(gè)特征值趨于平穩(wěn)，特征值為0.037；并且前7 個(gè)主元的累積貢獻(xiàn)率達(dá)到88%，超過(guò)85%。因此，選取前7個(gè)主元，這樣主元成分系數(shù)就變?yōu)槭?18)：

對(duì)于一個(gè)新的樣本，同樣按照聯(lián)合特征向量的方法，構(gòu)造出15 維聯(lián)合特征向量，然后向新的主元投影，得到降維后的特征向量Tnew∈R1×7。圖11分別為提取的融合特征向量Tnew中不同主元的二維和三維分布圖。

分析圖11(a)和圖11(b)可清晰地看出，PCA1和PCA2 主元區(qū)分的4 種狀態(tài)分布較開(kāi)，4 個(gè)狀態(tài)有較多的混疊，且在第一狀態(tài)和第四狀態(tài)下會(huì)出現(xiàn)一個(gè)樣本判別不了，而PCA1 和PCA3 主元區(qū)分的4 種狀態(tài)，在第一和第二狀態(tài)分布較為緊湊，第一、第三及第四狀態(tài)都出現(xiàn)了部分樣本判別不了。結(jié)合圖11(c)進(jìn)一步分析可知，通過(guò)使用PCA 降維后，4 種狀態(tài)的三維特征分布較開(kāi)，各狀態(tài)有交叉混疊，然而在第一狀態(tài)下，數(shù)據(jù)樣本錯(cuò)判及漏判嚴(yán)重。

3.2.2 LLE算法的降維特征融合實(shí)驗(yàn)

采用LLE 算法對(duì)矩陣M進(jìn)行特征降維，由于LLE 算法受低維子空間維數(shù)d和近鄰參數(shù)k選取的影響，本文采用極大似然估計(jì)法確定砂巖破裂各狀態(tài)特征數(shù)據(jù)集的低維子空間維數(shù)d。通過(guò)構(gòu)建近鄰間距離的似然函數(shù)，可得到全局?jǐn)?shù)據(jù)結(jié)構(gòu)的低維特征。

設(shè)xi為Rm獨(dú)立同步的觀察樣本，yi為Rd中呈光滑密度f(wàn)嵌入流形，則存在xi=g(yi)。假定Sx(t)是以x為球心、t為半徑的小球體，構(gòu)造二項(xiàng)式非其次過(guò)程，即

式(19)中，I{xi ∈Sx(r)}為示性函數(shù)，用泊松分布近似該過(guò)程可得

式(20)中，V(d)為d維空間下球的體積函數(shù)。令θ=lg(f(x))，則N(t)的似然函數(shù)可表示為

式(23)中，Tk(xi)是樣本點(diǎn)xi與第k個(gè)近鄰點(diǎn)的歐式距離。分別求解每個(gè)樣本點(diǎn)對(duì)應(yīng)的特征維數(shù)，并計(jì)算其平均值，將其當(dāng)作全局特征維數(shù)d，即

通過(guò)上述極大似然估計(jì)法計(jì)算砂巖聲發(fā)射信號(hào)特征數(shù)據(jù)集的低維子空間維數(shù)d=4，因此，利用LLE算法對(duì)樣本數(shù)據(jù)降維后得到4維的融合特征向量。

由于不同近鄰參數(shù)取值會(huì)對(duì)所得融合特征產(chǎn)生影響，為分析其不同取值的影響效果，先選擇固定本征維數(shù)d=4，然后改變近鄰參數(shù)k的取值，利用支持向量機(jī)(Support vector machines,SVM)對(duì)融合特征向量進(jìn)行特征識(shí)別，評(píng)價(jià)標(biāo)準(zhǔn)為各個(gè)階段的識(shí)別準(zhǔn)確率。圖12為不同近鄰參數(shù)k下各階段識(shí)別率曲線關(guān)系。

圖12 不同近鄰參數(shù)下各融合特征的敏感度之和Fig.12 The sum of the sensitivities of each fusion feature under different neighbor parameters

由圖12 中可以看出：隨著k值的變化，降維后所得到的融合特征經(jīng)SVM 識(shí)別后，識(shí)別率發(fā)生較大波動(dòng)，第二階段下識(shí)別率明顯偏高，且當(dāng)k=5時(shí)，4 個(gè)階段的識(shí)別率都達(dá)到最大，分別為0.84、0.88、0.84、0.8。因此確定近鄰參數(shù)k=5，此時(shí)LLE 算法降維后融合特征的聚類效果如圖13所示。

圖13 LLE 算法降維后的融合特征的聚類效果Fig.13 Clustering effect of fusion features after dimension reduction of LLE algorithm

從圖13(a)和圖13(b)可以看出：dim1 和dim2特征區(qū)分得第一階段與后三個(gè)階段分布相隔較遠(yuǎn)，而dim1 和dim3 特征區(qū)分的4 種狀態(tài)分布都靠的比較近，呈一條水平線趨勢(shì)，且各狀態(tài)交叉混疊數(shù)目較少，只有在第四狀態(tài)下存在一個(gè)樣本判別不了。結(jié)合圖13(c)可以看出，第一狀態(tài)沒(méi)有一個(gè)樣本錯(cuò)判，且4 個(gè)狀態(tài)相比于PCA降維后的聚類效果更集中。

4 兩種算法融合特征值的敏感度對(duì)比分析

為了更好地比較兩種降維方法的優(yōu)劣，利用散步矩陣法計(jì)算PCA 及LLE 算法在d=4、k=5參數(shù)下對(duì)矩陣M降維后的融合特征的敏感度之和。散布矩陣包含類間和類內(nèi)兩種散布矩陣，特征的散布值可由這兩個(gè)矩陣分析得到。設(shè)有M類破裂狀態(tài)，每種狀態(tài)的樣本數(shù)為N。原始特征向量x=(x1,x2,···,xD)，原始特征向量的維數(shù)用D表示。

類間散布矩陣：定義一個(gè)Sb，表達(dá)式如式(25)所示。特征的類間散布值越大，其類間辨識(shí)度則越好。

式(25)中，ui表示第i類特征值的均值；uo表示總體樣本的全局均值向量。

類內(nèi)散布矩陣：同理定義一個(gè)Sw，表達(dá)式如式(26)所示。特征的類內(nèi)散布值的大小，代表特征的類內(nèi)聚集度好壞，越小代表聚集度越好。

根據(jù)這兩種矩陣的分布特性，定義評(píng)價(jià)特征優(yōu)劣的敏感度算法，如式(27)所示：

式(27)中，tr(Sb)和tr(Sw)分別表示Sb和Sw的跡。當(dāng)Sb越大或Sw越小時(shí)，φ也越大，表明該特征對(duì)砂巖破裂狀態(tài)分類性能越強(qiáng)，反之則越弱。通過(guò)式(27)計(jì)算兩種算法下各融合特征的敏感度，結(jié)果如表4 所示。

表4 兩種算法降維后融合特征的敏感度之和Tabel 4 The sum of the sensitivities of the fusion features after dimensionality reduction of the two algorithms

由表4 可以看出：經(jīng)過(guò)LLE 算法降維得到的融合特征值的敏感度比PCA算法近乎兩倍大。前3個(gè)特征值的敏感度逐漸減小，第四個(gè)特征的敏感度增加，與PCA算法得到的融合特征值的敏感度的變化趨勢(shì)幾乎一樣，呈現(xiàn)先減小后增大，且融合特征敏感度之和遠(yuǎn)遠(yuǎn)大于PCA 算法，說(shuō)明經(jīng)過(guò)LLE 算法降維后得到的融合特征更多地表征了原始信號(hào)包含的局部信息。圖11(c)和13(c)也證明了LLE算法相比PCA算法具有更好的聚類效果，為后續(xù)砂巖破裂狀態(tài)識(shí)別分類的準(zhǔn)確性提供了有利的驗(yàn)證。

5 LLE特征融合下的砂巖破裂狀態(tài)分類實(shí)驗(yàn)驗(yàn)證

分別將得到的時(shí)域、頻域及融合特征向量分為兩組，用于訓(xùn)練和驗(yàn)證。訓(xùn)練樣本數(shù)為300 個(gè)，每類樣本為75 個(gè)；測(cè)試樣本數(shù)為100 個(gè)，每類樣本為25個(gè)。使用訓(xùn)練樣本的特征集訓(xùn)練基于PSO-LSSVM分類器，并優(yōu)化LSSVM 中參數(shù)σ和γ，設(shè)定PSO 算法中的局部搜索能力c1為1.5，全局搜索能力為c2為1.7，種群數(shù)量為50，迭代次數(shù)為200，迭代曲線如圖14所示。

圖14 各域特征迭代曲線圖Fig.14 Iterative graph of each domain feature

由圖14 可知，在時(shí)域特征集下搜索到的最佳懲罰因子γ=0.1，σ=207.267；在頻域特征集下搜索到的最佳懲罰因子γ=0.1，σ=105.438；在融合特征集下搜索到的最佳懲罰因子γ=3.672，σ=55.924。進(jìn)一步對(duì)比分析得到融合特征下迭代適應(yīng)度更集中，分布在90～98之間。

將測(cè)試樣本特征集放入訓(xùn)練好的PSO-LSSVM多模式分類器中進(jìn)行分類預(yù)測(cè)。由分類器的輸出值來(lái)確定巖石破裂狀態(tài)類型。設(shè)定1 表示壓密階段，2表示彈性階段，3 表示失穩(wěn)破壞階段，4 表示失穩(wěn)破壞后階段。圖15 為時(shí)域、頻域及融合特征向量集經(jīng)PSO-LSSVM算法分類的結(jié)果圖。

圖15 各域特優(yōu)化后的結(jié)果分類圖Fig.15 Result classification diagram after optimization of each domain

圖15 中藍(lán)色*型和紅色°型分別表示測(cè)試樣本預(yù)測(cè)與實(shí)際類別。由圖15 可知：各域特征集經(jīng)PSO-LSSVM 優(yōu)化分類，只出現(xiàn)少數(shù)樣本錯(cuò)分，而LLE 降維后得到的融合特征集經(jīng)PSO-LSSVM 優(yōu)化后識(shí)別效果最理想，計(jì)算各階段的分類識(shí)別率如表5所示。

表5 不同特征的識(shí)別率比較Tabel 5 Comparison of recognition rates of different features

對(duì)比時(shí)域特征、頻域特征及融合兩域特征的識(shí)別率，采用融合特征下的總識(shí)別率為90%，較單一的時(shí)域特征識(shí)別提高了6%，可知，使用LLE 算法對(duì)信號(hào)的聯(lián)合特征進(jìn)行數(shù)據(jù)降維用于巖石破裂狀態(tài)分類可以有效提高系統(tǒng)識(shí)別率。同時(shí)該方法可以將特征中次要成分去掉，降低了特征向量的維數(shù)，大大減少了后續(xù)訓(xùn)練及識(shí)別的時(shí)間，避免了“維數(shù)災(zāi)難”的發(fā)生。

6 結(jié)論

(1) 通過(guò)單軸壓縮試驗(yàn)，計(jì)算紅砂巖破裂失穩(wěn)過(guò)程的應(yīng)力值和應(yīng)變值隨時(shí)間變化的關(guān)系，得到不同砂巖失穩(wěn)過(guò)程的應(yīng)力應(yīng)變曲線，通過(guò)比較分析，將砂巖破裂失穩(wěn)劃分為4 個(gè)階段，并分析各階段下巖石內(nèi)部狀態(tài)變化，為預(yù)測(cè)砂巖破壞失穩(wěn)提供依據(jù)。

(2) 將構(gòu)造的時(shí)、頻域聯(lián)合多維特征向量采用PCA 方法和流形學(xué)習(xí)LLE 算法進(jìn)行特征約簡(jiǎn)，比較發(fā)現(xiàn)LLE 算法降維后融合特征向量的聚類效果更好，樣本錯(cuò)判漏判現(xiàn)象較少，且LLE 算法降維后的融合特征的敏感度之和遠(yuǎn)大于PCA方法，表明該融合特征更多地包含了原始信號(hào)的局部特征信息。

(3) 利用粒子群優(yōu)化最小二乘支持向量機(jī)算法進(jìn)行巖石破裂狀態(tài)分類，比較分析基于粒子群改進(jìn)的最小二乘支持向量機(jī)算法對(duì)單一特征與經(jīng)LLE算法降維融合特征下的分類識(shí)別率，可以得出采用融合特征后的分類識(shí)別率顯著提高且識(shí)別的效果明顯要優(yōu)于單一特征識(shí)別。表明，充分考慮多個(gè)域特征向量之間相互制約和影響的關(guān)系相較于單一特征向量，在巖石破裂狀態(tài)分類方面具有更為準(zhǔn)確的預(yù)測(cè)結(jié)果。