孫 遠， 楊 峰， 鄭 晶， 徐茂軒， 裴爍瑾

(1.中國礦業(yè)大學(北京) 機電與信息工程學院，北京 100083；2. 中國礦業(yè)大學(北京) 煤炭資源與安全開采國家重點實驗室，北京 100083)

隨著信號處理技術(shù)的發(fā)展，盲源分離(Blind Source Separation, BSS)又稱盲信號分離，得到了越來越廣泛的關(guān)注與研究。盲源分離是在無法獲取源信號和傳遞函數(shù)先驗知識的前提下，將混合后的各源信號分離出來的過程。它主要廣泛應(yīng)用于醫(yī)學信號處理、經(jīng)濟指標、語音信號識別，圖像處理、地震信號勘測以及通信信號處理等方面[1]。盲源分離中較為常用的假設(shè)是源信號的統(tǒng)計獨立性，當各個分量零均值且方差為1即相互獨立時，就稱為獨立成分分析(Independent Component Analysis, ICA)[2]。

ICA的目的是從混合后的多路觀測信號中，將混合前的源信號分離出來。例如在嘈雜的環(huán)境中提取出關(guān)心的某些聲音、多個傳感器記錄的腦電波、機械振動信號的分析與降噪等。ICA實質(zhì)是一種在某一判據(jù)條件下的尋優(yōu)計算，因此，ICA在分離信號的過程中主要進行優(yōu)化判據(jù)和尋優(yōu)算法兩方面工作。在ICA尋優(yōu)算法的選取上，較為傳統(tǒng)的方法是梯度法，梯度法容易受步長以及初值的影響，導致收斂速度較慢。針對這一問題，近些年有關(guān)群智能優(yōu)化算法被廣泛的應(yīng)用與ICA的尋優(yōu)算法中。文獻[3]將粒子群優(yōu)化算法(Particle Swarm Optimization, PSO)應(yīng)用到ICA中，較好的識別出滾動軸承故障的類型。然而粒子群優(yōu)化算法存在早熟收斂的缺點，粒子在全局搜索過程中遇到最優(yōu)值后容易陷入局部最優(yōu)值。文獻[4]將人工蜂群算法(Artificial Bee Colony, ABC)應(yīng)用到BSS中，成功的將混合信號分離開來，加快了收斂速度。但該算法含有較多參數(shù)，參數(shù)選取不當容易導致算法陷入局部極值點，降低分離效果。文獻[5]將模擬退火算法(Simulated Annealing Algorithm, SAA)與ICA相結(jié)合，提高了ICA的穩(wěn)態(tài)性能。但模擬退火算法的初始溫度參數(shù)設(shè)置較難控制，當初始溫度設(shè)置較大時，ICA算法需要消耗大量時間，當初始溫度設(shè)置較小時，則會影響全局搜索能力。文獻[6]將遺傳算法(Genetic Algorithm, GA)作為ICA的優(yōu)化算法，有效的分離了非高斯源信號。但遺傳算法局部搜索能力較弱，并且實現(xiàn)效果取決于問題的多種參數(shù)。

為了解決算法的早熟收斂、參數(shù)復(fù)雜、分離精度不高的問題，提出一種基于膜計算與粒子群算法的盲源分離方法，可以保證PSO算法局部搜索的精度，滿足全局搜索的多樣性，并且能避免PSO算法的早熟收斂問題，提高收斂速度和分離性能。

1 獨立成分分析原理

線性瞬時混疊ICA模型的原理圖，如圖1所示。

圖1 線性瞬時混疊ICA模型的原理圖Fig.1 Principle of linear instantaneous mixing ICA model

圖1中，s1(t)，s2(t)，…，sn(t)表示n個源信號，x1(t)，x2(t)，…，xm(t)表示n個源信號經(jīng)過線性混合后由m個傳感器所接收到的觀測信號在t時間的幅值，線性方程表達式如下

(1)

式中：aij(i=1,2,…,m,j=1,2,…,n)是理想狀態(tài)下xi(t)與si(t)之間的混合參數(shù)。

ICA的簡單模型如圖2所示。

圖2 ICA簡單模型

Fig.2 Simple model of ICA

圖2中，s(t)=[s1(t),s2(t),…,sn(t)]T代表沒有先驗知識的源信號所組成的N維向量；A是M×N階混合矩陣；x(t)=[x1(t),x2(t),…,xm(t)]T即為M維觀測量；B是分離矩陣，y(t)=[y1(t),y2(t),…,yn(t)]T是傳感器接收到的觀測信號經(jīng)過B分離后得到的輸出量。ICA要解決的問題就是，在s(t)和A不具備先驗知識的情況下，通過對分離矩陣B的求解與更新，得到輸出信號y(t)，盡量保證y(t)中每個分量相互獨立[7]。此時，可以將輸出信號y(t)作為s(t)的有效估計值。

2 粒子群優(yōu)化算法

粒子群優(yōu)化算法是一種進化計算技術(shù)(Evolutionary Computing)，主要依靠各粒子之間的相互協(xié)作來完成對最優(yōu)解的尋找。在尋找最優(yōu)解的過程中，每個粒子根據(jù)與目標函數(shù)相關(guān)的適應(yīng)值(Fitness Value)對自己當前位置進行評價，找到自己在尋優(yōu)過程中經(jīng)歷的最好位置pbest，作為粒子的局部最優(yōu)值。另外，還需要根據(jù)計算適應(yīng)值找到所有粒子在尋優(yōu)過程中的群體最優(yōu)位置，即所有粒子發(fā)現(xiàn)的最好位置gbest(gbest是pbest中最好的值)。將此群體最優(yōu)值作為所有粒子的經(jīng)驗，根據(jù)式(2)、(3)來更新粒子的速度和位置。

vij(t+1)=wvij(t)+c1rand1(pij-xij)+

c2rand2(gij-xij)

(2)

xij(t+1)=xij(t)+vij(t+1)

(3)

式(2)、(3)中：i=1,2,…,N是粒子群的總數(shù)；j=1,2,…,D是種群搜索空間的維數(shù)；xij表示第i個粒子第j維分量的位置；vij表示第i個粒子第j維分量的速度；w為慣性權(quán)重；c1和c2是學習因子；rand1和rand2是介于(0,1)之間的隨機數(shù)；pij表示第i個粒子第j維分量當前最好的位置；gij表示所有粒子第j維分量當前最好的位置。在粒子尋優(yōu)過程中，粒子的速度vi的最大值為Vmax(Vmax>0)，如果vi>Vmax，則vi=Vmax。

3 基于MCPSO的ICA算法

從信息論角度講，信息極大化、互信息極小或極大似然估計均可以作為信號獨立性的判據(jù)。根據(jù)中心極限定理可知，當隨機變量X可以由若干相互獨立的隨機變量Si(i=1,2,…,N)之和表示時，對任意分布的Si，如果Si滿足有限的均值和方差這一條件，則Si與X相比，具有更強的非高斯性。因此，在ICA中，可以用非高斯性來衡量輸出信號y(t)中各分量的相互獨立性，非高斯性越強，表明分離的效果越好。在信息論中，負熵是度量非高斯性的一個較為穩(wěn)健的判據(jù)[8]，如式(4)所示

J(y)=HG(y)-H(y)

(4)

根據(jù)文獻[9]證明，多變量負熵可以近似表示為

(5)

式中：k3表示信號的三階累積量；k4表示信號的四階累積量。當信號的概率密度分布對稱時，k3=0，則

(6)

(7)

對某一分離矩陣B，J[y]越大，表明對源信號s(t)=[s1(t),s2(t),…,sn(t)]T的分離效果越好。由式(7)可知，尋找J[y]的最大值，對應(yīng)于PSO算法中尋找fiti的最小值。當用負熵來衡量信號的非高斯程度時，需要滿足信號為零均值以及E(yyT)=I(I為單位矩陣)的約束條件。因此，首先需要對信號進行中心化和白化處理，以滿足約束條件[11]。

在ICA優(yōu)化處理過程中，當用梯度法對分離矩陣B進行更新時，基于信息極大化、互信息極小化或極大似然估計目標函數(shù)，可以使用如下步長函數(shù)[12]

B(k+1)=B(k)+

μk[I-Ψ(y(k))yT(k)]B(k)

(8)

式中：k為更新時的迭代次數(shù)；μk>0為時變的調(diào)節(jié)步長；I為單位矩陣；B(k)為第k次迭代的分離矩陣；Ψ(y(k))為y(k)的非線性奇函數(shù)，本文選擇為基于皮爾遜系統(tǒng)的分段非線性函數(shù)[13]。其表達式如下

(9)

式中：a，b0分別為皮爾遜系統(tǒng)的參數(shù)；α3，α4為信號的三階矩和四階矩，可以根據(jù)矩估計的方法求出。

針對傳統(tǒng)的梯度法易受步長影響，導致ICA算法的收斂速度較慢，文獻[14]采用粒子群優(yōu)化算法調(diào)整ICA算法的步長μk，根據(jù)式(8)更新分離矩陣B，可以減小ICA算法的穩(wěn)態(tài)誤差，提高收斂速度。但粒子群優(yōu)化算法存在著早熟問題，即在迭代更新過程中，粒子會根據(jù)當前其它粒子提供的最優(yōu)位置的經(jīng)驗，均飛向最優(yōu)位置，忽略其它位置的搜索，此時粒子的飛行速度會大大降低，甚至停止搜索，這導致粒子的全局搜索性能變差，過早的出現(xiàn)收斂。為解決早熟收斂的現(xiàn)象，常采用調(diào)節(jié)慣性權(quán)重w的方法來平衡PSO算法的全局搜索和局部搜索。文獻[15]提出線性遞減慣性權(quán)重策略(LDIW)，在粒子進行全局搜索時，w取較大值，使粒子可以在更廣的范圍進行搜索，在粒子進行局部搜索時，w取較小值，確保粒子能夠進行精細化的搜索。為增加迭代初期的全局搜索時間和迭代后期的局部搜索時間。文獻[14]提出一種新的非線性遞減慣性權(quán)重策略(NDIW)，解決了PSO的早熟收斂問題。但是，無論LDIW還是NDIW方法，都會增加計算的復(fù)雜度，降低ICA算法的收斂速度。考慮到膜計算具有并行處理的特點，可以簡化慣性權(quán)重的取值問題。因此，提出基于膜計算與粒子群算法(MCPSO)的盲源分離方法。

3.1 MCPSO算法

膜計算是一門基于計算機科學和自然科學的新的計算方法，主要受細胞的啟發(fā)，并從中提取出計算模型[16]。膜系統(tǒng)主要由膜結(jié)構(gòu)、對象和進化規(guī)則三部分構(gòu)成。膜計算的顯著特點在于并行性、分布式和非確定性。

MCPSO算法采用細胞型膜系統(tǒng)的單層膜結(jié)構(gòu)(One Level Membrane Structure, OLMS)，其簡單結(jié)構(gòu)示意圖如圖3所示。將粒子均勻地放入若干基本膜中，在每個基本膜內(nèi)，粒子進行精細化的局部搜索。圖4為MCPSO算法膜結(jié)構(gòu)示意圖。每個粒子的位置看作是膜區(qū)域內(nèi)的對象。將D維空間的n個粒子均勻的分配到m個基本膜中，表層膜為空。此膜系統(tǒng)的字符表示如下多元組

圖3 細胞型膜系統(tǒng)結(jié)構(gòu)示意圖Fig.3 Cell membrane system structure diagram Π=(V,O,H,μ,w0,…,wm,R0,…,Rm,io)

(10)

式中：V={xij(i=1,2,…,N;j=1,2,…,D)}；O?V是輸出對象的集合；H={0,1,…,m}；μ=[[]1[]2[]3…[]m]0；w0,…,wm是對應(yīng)的膜內(nèi)的對象多重集；R0,…,Rm表示各膜內(nèi)區(qū)域?qū)?yīng)的規(guī)則；io是膜系統(tǒng)輸出膜的標號。

圖4 MCPSO膜結(jié)構(gòu)示意圖Fig.4 Membrane system structure diagram of MCPSO

圖5 MCPSO算法過程示意圖Fig.5 Algorithm process diagram of MCPSO

位于m個基本膜內(nèi)的粒子進行局部搜索，由于膜計算具有并行處理能力，所有基本膜內(nèi)對粒子最優(yōu)位置的搜索同時進行，在各基本膜進行局部搜索的同時，對于整個粒子群來說是在進行全局搜索。膜計算的引入，可以提高粒子局部搜索的精度，保證了全局搜索的多樣性，避免種群搜索過程中陷入局部最優(yōu)值。由于MCPSO算法粒子局部搜索和全局搜索同時進行，w只需要取較小值，來確保局部搜索的精度即可，膜計算并行處理的特點，簡化了慣性權(quán)重w的取值問題，降低了計算的復(fù)雜度，提高了PSO算法的收斂速度。

3.2 基于MCPSO的ICA算法流程

基于膜計算和粒子群優(yōu)化的ICA算法進行信號分離的流程為

步驟1對觀測信號進行中心化和白化處理；

步驟2初始化一個單層膜結(jié)構(gòu)，其中包括一個表層膜0和m個基本膜；初始化一個由n個粒子組成的D維空間里按照一定速度飛行的粒子種群，將所有粒子均勻的分配到m個基本膜中，表層膜為空；設(shè)定學習因子c1、c2、慣性權(quán)重w和最大迭代次數(shù)Tmax，隨機產(chǎn)生分離矩陣和粒子的飛行速度；

步驟3根據(jù)y(t)=Bx(t)計算輸出信號y(t)，并對y(t)進行中心化和白化處理，根據(jù)式(6)、(7)計算各基本膜內(nèi)粒子的適應(yīng)度；

步驟4將每個粒子的適應(yīng)度值與該粒子經(jīng)過的最好位置pbest做比較，選擇適應(yīng)值較小的粒子作為當前的個體最優(yōu)位置；將每個粒子的個體最優(yōu)位置與其所在基本膜內(nèi)的群體最優(yōu)位置gbest′做比較，選擇適應(yīng)值較小的粒子作為當前基本膜內(nèi)群體最優(yōu)位置；

步驟5將m個基本膜內(nèi)的群體最優(yōu)位置通過規(guī)則R送至表層膜內(nèi)，再通過對比此m個群體最優(yōu)位置，找到適應(yīng)值最小的一個最優(yōu)位置gbest，將gbest作為所有粒子的群體最優(yōu)位置，并將其由表層膜返回至各基本膜；

步驟6根據(jù)式(2)、(3)對所有粒子進行速度和位置的更新；

步驟7如果滿足終止條件(找到適應(yīng)度最小的粒子或者達到迭代次數(shù))，則停止迭代，輸出最優(yōu)解，否則返回步驟3；

步驟8將步驟7輸出的最優(yōu)解作為式(8)的步長μk，更新分離矩陣B；

步驟9根據(jù)y(t)=Bx(t)提取分離信號。

3.3 偽代碼

Begin

generation←1

Centralize and sphere observed signal x(t)//中心化、白化觀測信號

Initialize skin and elementary membrane structure //初始化表層膜和基本膜結(jié)構(gòu)

Generate n particles in each elementary membrane uniformly//初始化粒子群并均勻分配至基本膜

Initialize separable matrix randomly//隨機初始化分離矩陣

Compute output signal y(t) and then centralize and sphere it//計算輸出信號并中心化和白化

while(not termination condition) do//終止條件

—Evaluate every particle in each elementary membrane//計算粒子適應(yīng)度

—Find local best particle pbest //找到粒子個體最優(yōu)值

—Find local best particle gbest//找到各基本膜內(nèi)群體最優(yōu)值

—Execute communication rules//執(zhí)行膜內(nèi)規(guī)則

—Send global best from all elementary membranes to skin membrane//基本膜內(nèi)群體最優(yōu)值送到表層膜

—Chose best value of all received global best as global best among all particles//選擇表層膜內(nèi)最優(yōu)值作為群體最優(yōu)值

—Send best value back to each elementary membrane//表層膜內(nèi)群體最優(yōu)值送回各基本膜

—Update particle’s velocityV(t)//更新粒子速度

—Update particle’s positionX(t)//更新粒子位置

—generation←generation +1

end

Output best position//輸出最優(yōu)解

Update separable matrixB//更新分離矩陣

Extract separable signal//計算分離信號

end

4 實驗結(jié)果及分析

4.1 仿真實驗

為驗證本文提出的MCPSO算法的有效性，實驗采用一組正弦波、一組方波、一組鋸齒波作為源信號，如圖6所示。隨機產(chǎn)生線性瞬時混合矩陣A如下，源信號經(jīng)線性瞬時混合矩陣A混合后的信號如圖7所示

圖6 源信號Fig.6 Source signals

仿真實驗在MATLAB 2012b下運行，用提出的MCPSO算法和文獻[3]中使用線性慣性權(quán)重的PSO算法(以下簡稱線性PSO算法)以及文獻[14]中使用非線性慣性權(quán)重的PSO算法(以下簡稱非線性PSO算法)分別對源信號進行分離。MCPSO算法的控制參數(shù)為：粒子群個數(shù)N=30,搜索空間D=5，基本膜個數(shù)m=5，粒子速度范圍v=[-1,1]，c1=c2=2，慣性權(quán)重w=0.4，最大迭代次數(shù)Tmax=100。粒子數(shù)與基本膜的個數(shù)是通過實驗對比分析，找到MCPSO算法分離效果最好時對應(yīng)的粒子數(shù)與基本膜個數(shù)。分離效果采用相似系數(shù)累衡量，其表達式如下

(11)

式中：cov(·)表示方差，si為源信號矢量s中的第i個源信號，yj為經(jīng)過盲源分離后的與si相對應(yīng)的分離信號。相似系數(shù)越接近1，源信號與分離信號越接近，分離效果越好。

圖7 混合信號Fig.7 Mixed signals

對于通常問題，30個粒子已足夠，對于較復(fù)雜的問題，粒子數(shù)可取為50，當種群粒子數(shù)小于50時，可發(fā)現(xiàn)種群大小對PSO性能有顯著的影響，當種群粒子數(shù)大于50時，對PSO性能影響較小[17]。因此在對比實驗分析中，粒子個數(shù)分別選為10、20、30、40、50、60、70；基本膜個數(shù)選擇為從2到8的整數(shù)。如圖8、9所示，分別為不同粒子數(shù)和基本膜個數(shù)的分離效果曲線。從圖8、9可以看出當粒子數(shù)為30，基本膜個數(shù)為5時，相似系數(shù)的值最大，表明MCPSO算法的分離效果最好。

圖8 基本膜個數(shù)對分離效果的影響Fig.8 The influence of membrane number to separation effect

圖9 粒子數(shù)對分離效果的影響Fig.9 The influence of particle number to separation effect

圖10 MCPSO分離信號Fig.10 Separated signals of MCPSO

圖11 線性PSO分離信號Fig.11 Separated signals of linear PSO

從圖10、11、12分離結(jié)果圖可以看出，線性PSO算法分離后的正弦信號、方波信號、鋸齒波信號均有變形情況，分離效果不夠理想；非線性PSO算法分離后的方波信號有變形情況，但正弦信號和鋸齒波信號分離效果較好；提出的MCPSO算法分離出的3組信號均與源信號有較好的吻合，分離效果較為理想，提高了分離精度。

圖12 非線性PSO分離信號Fig.12 Separated signals of nonlinear PSO

算法的分離效果的好壞可以選用性能指數(shù)PI(Performance Index)來衡量，其定義如下

(12)

式中：g(i,j)為矩陣G的第(i,j)個元素，而G為全局矩陣，即分離矩陣與混合矩陣的乘積，PI值越小表明分離的效果越好。圖13為MCPSO算法和線性PSO算法以及非線性PSO算法性能指數(shù)PI的收斂曲線比較結(jié)果。從圖13可以看出，MCPSO算法在迭代28步時即可收斂，這是由于膜計算并行處理的優(yōu)點，提高了算法的收斂速度，降低了計算復(fù)雜度；線性PSO算法在迭代68步后才開始收斂；非線性PSO算法在迭代62步開始收斂，雖然非線性PSO算法基本可以將源信號分離出來，但是在收斂速度上與MCPSO算法相比較低。因此，MCPSO算法在收斂速度上要優(yōu)于線性PSO算法和非線性PSO算法。

圖13 PI收斂曲線Fig.13 Convergence curves of performance index

為了定量的說明信號的盲源分離效果，需要多種評價指標來衡量，來從不同角度反映分離信號與源信號之間的誤差度量[18]。本文采用相似系數(shù)ρij、二次殘差VQM以及性能指數(shù)PI三種評價指標。

二次殘差采用帶幅值修正因子的計算公式[19]，表達式如下

(13)

表1為MCPSO算法、線性PSO算法以及非線性PSO算法的盲源分離評價指標的比較。從表1可知，MCPSO算法與線性PSO算法和非線性PSO算法相比，相似系數(shù)更接近于1，二次殘差VQM更小，PI性能指數(shù)更小。表明MCPSO算法的盲源分離效果最好。

表1 盲源分離評價指標對比Tab.1 Evaluation indexes comparison of blind source separation

4.2 應(yīng)用實例

本節(jié)利用所提出的算法處理人員腳步和車輛實測振動信號。采集裝置使用Geopen地震儀，如圖14所示。采樣頻率設(shè)置為1 kHz，采樣點數(shù)為16×103，采樣率設(shè)置為1 ms。實驗主要分為三組：第一組實驗設(shè)計為單獨采集單人行走腳步振動信號；第二組實驗設(shè)計為單獨采集車輛(漢蘭達2.0T)振動信號；第三組實驗設(shè)計為同時采集第一組試驗中的人員和第二組實驗中的車輛振動混合信號。前兩組實驗的設(shè)計是為了與分離出的信號作對比，證明MCPSO算法的有效性和準確性。如圖15所示，分別為人員腳步振動信號和車輛振動信號的時域波形。圖16所示為傳感器的觀測信號，即采集到的單人行走腳步振動和車輛振動的混合信號。圖17為MCPSO算法對觀測信號進行分離后的信號。圖18為人員腳步振動源信號與分離后信號的頻譜對比。圖19為車輛振動源信號與分離后信號的頻譜對比。

圖14 采集實驗Fig.14 The experiment of vibration signal collection

圖15 人員腳步和車輛振動信號時域波形Fig.15 The time domain waveform of step and vehicle vibration signals

圖16 傳感器觀測信號Fig.16 The sensor signals

圖17 分離信號Fig.17 The separated signals

圖18 人員腳步振動源信號與分離信號的頻譜對比Fig.18 Spectrum comparison between the step source signals and the separated signals

圖19 車輛振動源信號與分離信號的頻譜對比Fig.19 Spectrum comparison between the vehicle source signals and the separated signals

為使結(jié)果更可靠，對觀測信號多次分離，計算時頻域評價指標，取多次分離實驗的平均值作為最終的評價指標。如表2所示，三種時域評價指標均能反應(yīng)出MCPSO算法的分離效果較好。對人員腳步振動源信號和分離信號以及車輛振動源信號和分離信號做頻譜分析，兩種分離信號的頻率與兩種源信號的頻率基本相同，腳步振動源信號和分離信號頻率主要在50 Hz以下，車輛振動源信號和分離信號頻率主要在200 Hz以下。

表2 MCPSO的盲源分離時頻域評價指標Tab.2 The time-frequency domain evaluation indexes of MCPSO blind source separation

5 結(jié) 論

本文將膜計算引入到粒子群優(yōu)化算法中來進行信號的盲源分離，將粒子放入基本膜中，通過計算適應(yīng)度，將各基本膜中粒子最優(yōu)位置通過膜規(guī)則輸送到表層膜中，再將表層膜中粒子最優(yōu)值作為群體最優(yōu)值返回給各基本膜來對粒子的位置和速度進行更新，在保證粒子局部搜索精度的同時，也滿足了粒子全局搜索的多樣性。膜計算的引入簡化了慣性權(quán)重的取值問題，只需要考慮局部搜索時慣性權(quán)重的取值，降低了計算復(fù)雜度，提高了算法的收斂速度，有效地解決了PSO全局搜索和局部搜索之間的矛盾。MCPSO算法可以提高信號的分離精度，且具有較快的收斂速度，有效的避免了標準PSO算法粒子容易陷入早熟的問題。通過仿真實驗對比和實例應(yīng)用驗證了MCPSO算法的有效性。