密度空間與密度峰值聚類的欠定混合矩陣估計

何選森 何 帆

（1.廣州商學(xué)院信息技術(shù)與工程學(xué)院，廣州，511363；2.湖南大學(xué)信息科學(xué)與工程學(xué)院，長沙，410082；3.北京理工大學(xué)管理與經(jīng)濟學(xué)院，北京，100081）

引 言

在信源和傳輸信道均未知的情況下,僅利用傳感器采集獲得的觀測信號來分離信源的過程稱為盲源分離(Blind source separation,BSS)[1]。在無噪情況下,BSS的時域模型為:x(t)=As(t),其中s(t)=[s1(t),s2(t),…,sM(t)]T為源信號向量,x(t)=[x1(t),x2(t),…,xN(t)]T為觀測信號向量,A∈RN×M為混合矩陣。當(dāng)M=N,稱BSS為適定的;若M＜N則稱BSS為超定的;若M＞N,則稱為欠定的盲源分離(Underdetermined BSS,UBSS)[1]。在UBSS問題中,由于混合矩陣的逆不存在,對它的估計是很困難的[2],且辨識混合矩陣與分離信源成為兩個截然不同的問題。因此,兩步法[3]成為UBSS的主要解決方案,其中第一步是估計混合矩陣,第二步是分離源信號。而至關(guān)重要的是對混合矩陣的估計,因為它的結(jié)果直接影響到BSS的性能。

近些年,利用群體智能算法解決UBSS問題成為國內(nèi)外學(xué)者的研究熱點,主要包括蟻群優(yōu)化算法[4],人工蜂群算法[5],粒子群算法[6]等;而對于單通道的UBSS,基于粒子濾波[7]以及粒子流濾波[8]的方法在改進(jìn)UBSS性能方面也取得了一定的進(jìn)展。另外,充分利用自然信號(如音頻,圖像等)本身具有的稀疏特性來解決UBSS問題已成為業(yè)界的共識。因此,稀疏表示[9]和稀疏分量分析[10]是處理UBSS問題最基礎(chǔ)和最有效的方法。由于稀疏信號具有線性聚類特性，且聚類形成的直線方向向量就是混合矩陣的列向量[11]，因此對觀測信號進(jìn)行聚類分析能實現(xiàn)混合矩陣的估計。常用的聚類方法有K-均值[12],霍夫變換[13],勢函數(shù)[3]等。另外,把兩種不同聚類方法組合能形成更有效的聚類分析。例如:算法KHough[14]利用霍夫變換對K-均值的聚類中心進(jìn)行修正,同時K-Hough還克服了霍夫變換在進(jìn)行峰值提取時的峰值簇?fù)韱栴};DBSCAN-Hough[15]采用具有噪聲的基于密度空間聚類(Density based spatial clustering of applications with noise,DBSCAN)算法[16]確定聚類的數(shù)目,通過霍夫變換對聚類中心進(jìn)一步修正。正是從這些組合的方法中得到啟發(fā),本文利用單源點(Single-source-point,SSP)檢測以增加信源的線性聚類特性,然后把DBSCAN算法和快速搜索與尋找密度峰值聚類(clustering by fast search and find of density peaks,CFSFDP)算法[17]相結(jié)合形成新的聚類分析方法。選擇DBSCAN是因為它能對觀測數(shù)據(jù)進(jìn)行聚類以自動確定信源數(shù)目,從而克服了K-均值算法需事先預(yù)知信源數(shù)目的缺陷。然而,DBSCAN算法對于高維數(shù)信號,或密度不均勻、聚類間距相差很大的信源,其聚類效果不理想。為此,在DBSCAN聚類分析的基礎(chǔ)上,利用CFSFDP對每一類數(shù)據(jù)分別計算出對應(yīng)的密度峰值,并把峰值點作為修正后的聚類中心,從而提高混合矩陣的估計精度。把DBSCAN和CFSFDP兩種聚類算法相結(jié)合的另一個優(yōu)勢是DBSCAN能彌補CFSFDP算法需要人為干預(yù)的不足。

1 單源點檢測

對于時域BSS模型x(t)=As(t),兩邊取短時傅里葉變換(Short-time Fourier transform,STFT),則

式中：X(t,k)=[X1(t,k),X2(t,k),…,XN(t,k)]T和S(t,k)=[S1(t,k),S2(t,k),…,SM(t,k)]T分別為x(t)∈RN和s(t)∈RM在時頻點(t,k)的STFT的系數(shù)；am為混合矩陣A的第m個列向量。與時域相比,時頻域中稀疏信號的直線聚類特性得到了更好的體現(xiàn),但仍存在有一些觀測數(shù)據(jù)不能聚集在直線上,而是處于直線之外,這就造成混合矩陣的估計性能下降。為此,本文采用SSP檢測[18]。所謂SSP是指在這些時頻點上只有一個主導(dǎo)信源的能量具有較大的值,而其余信源能量很小以至于可忽略。而不滿足SSP條件的時頻點稱為多源點(Multi-source-points,MSP)。顯然,SSP檢測后的數(shù)據(jù)具有顯著直線聚類的方向性。

對于任意一個時頻點(t,k),觀測信號X(t,k)的實部和虛部分別為

這二者之間的夾角θ為

如果信源各分量的實部與虛部之比是相等的,即滿足以下關(guān)系

則它們之間的夾角θ=0°或θ=π(180°)。這時,就稱該時頻點(t,k)是 SSP。

在實際應(yīng)用中,恒等式(4)成立的條件是很苛刻的。為此可采用另一個判斷條件[18]:若R[X(t,k)]和I[X(t,k)]的絕對方向相同,則對應(yīng)的時頻點(t,k)就稱為SSP。通常情況下,若R[X(t,k)]和I[X(t,k)]的絕對方向夾角小于某個閾值Δθ時,就認(rèn)為時頻點(t,k)是SSP,即

式中：|z|表示z的絕對值,而||Z||=(ZTZ)1/2。本文采用的閾值為Δθ=0.8°。

通過SSP檢測之后,由于剔除了MSP,則觀測信號的數(shù)據(jù)分布凸顯出了明確的方向性。由所有SSP的數(shù)據(jù)點組成的集合記為Xssp。

一般地,經(jīng)過原點的直線方向可以被兩個方向相反的向量來表示。例如,在三維空間中同一條直線可被方向向量(1,1,1)或(-1,-1,-1)來描述。為了用唯一的方向向量描述直線,采用鏡像映射[4]方式,將負(fù)方向的向量映射到對應(yīng)的正方向上。在信號處理領(lǐng)域,利用觀測信號的歸一化可實現(xiàn)鏡像映射的過程[11]

從式(6)可知,鏡像映射是把線性聚類轉(zhuǎn)換成致密聚類,便于利用基于密度聚類方法搜索到密集數(shù)據(jù)堆中的關(guān)鍵數(shù)據(jù)點(聚類中心)。該數(shù)據(jù)點的方向就是稀疏信號線性聚類的直線方向,即混合矩陣的列向量。因此,通過對數(shù)據(jù)X*(k)進(jìn)行聚類分析就可實現(xiàn)對欠定混合矩陣的估計。

2 聚類分析

在眾多聚類方法中,DBSCAN作為基于密度聚類的典型代表,能夠在有噪聲的數(shù)據(jù)中發(fā)現(xiàn)各種形狀和各種大小的數(shù)據(jù)簇。DBSCAN的核心思想是在數(shù)據(jù)堆中找到密度較高的數(shù)據(jù)點,再通過搜索鄰近的其他高密度數(shù)據(jù)點,逐步將高密度數(shù)據(jù)點連成一片,從而生成各種形狀的數(shù)據(jù)簇。

DBSCAN算法是利用參數(shù)(eps,MinPts)來描述鄰域的樣本分布緊密程度。首先,以每個數(shù)據(jù)點為圓心,以鄰域eps為半徑畫個圓圈,落在該圈內(nèi)的數(shù)據(jù)點數(shù)就是該點的密度值。然后,利用密度閾值MinPts判斷數(shù)據(jù)點的密度級別,若圓圈內(nèi)數(shù)據(jù)點數(shù)小于MinPts,則其圓心的數(shù)據(jù)點是低密度點,而點數(shù)大于或等于MinPts的圓心數(shù)據(jù)點為高密度點(核心點Core point)。若某個低密度數(shù)據(jù)點落在高密度點的圓圈內(nèi),則把低密度點連到最鄰近的高密度點上,并稱它為高密度點的邊界點。不在任何高密度點圈內(nèi)的低密度點為異常點(噪聲)。

若某個樣本點y在點x的eps鄰域內(nèi),且x為核心點,則稱y從x直接密度可達(dá)。假設(shè)給定一連串?dāng)?shù)據(jù)樣本點x1,x2,…,xn且x=x1,y=xn,若xi+1從xi(i∈[1,n])直接密度可達(dá),則稱y從x密度可達(dá)。對于樣本點xi和xj,若存在核心點xk,使xi和xj都可以由xk密度可達(dá),則稱樣本點xi和xj密度相連。由密度可達(dá)關(guān)系得到的最大密度相連的樣本集合,即為聚類分析最終得到的一個類別(簇)。

DBSCAN的聚類效果取決于參數(shù)eps和MinPts的選取。MinPts的選取原則是:MinPts≥D+1,其中D為待聚類數(shù)據(jù)的維度;一般地,MinPts必須選擇大于等于3的值。參數(shù)eps的選擇也要適中:若eps值太小,會造成大部分?jǐn)?shù)據(jù)不能聚類;若eps值過大,會使多個數(shù)據(jù)簇被合并到同一個簇中。eps選擇可通過繪制K-距離曲線來實現(xiàn),在曲線的明顯拐點位置對應(yīng)于合適的eps參數(shù)值。

與傳統(tǒng)的K-均值聚類相比,DBSCAN算法的優(yōu)勢主要體現(xiàn)在:(1)可以對任意形狀的稠密數(shù)據(jù)集進(jìn)行聚類,而K-均值僅適用于凸數(shù)據(jù)集;(2)可以自動獲得聚類的數(shù)量,而K-均值需事先給定聚類數(shù);(3)在聚類的同時還能找出異常(噪聲)點;(4)聚類的結(jié)果沒有偏移,而K-均值的初始值對聚類結(jié)果影響很大。

然而,由于使用全局的密度閾值參數(shù)MinPts,DBSCAN只能發(fā)現(xiàn)密度值不少于MinPts的數(shù)據(jù)點所組成的簇,而很難發(fā)現(xiàn)不同密度的數(shù)據(jù)簇。為此,本文采用可視化的CFSFDP算法對DBSCAN產(chǎn)生的初步聚類中心進(jìn)行修正,以提高關(guān)鍵數(shù)據(jù)的定位精度。CFSFDP的基本思想是：由于每個數(shù)據(jù)簇都有一個最大密度的數(shù)據(jù)點作為簇中心,在它的周圍都是密度比它低的數(shù)據(jù)點;使得不同的數(shù)據(jù)簇的中心相距較遠(yuǎn),從而可區(qū)分出不同密度的數(shù)據(jù)簇。顯然,CFSFDP彌補了DBSCAN算法僅適合于稠密數(shù)據(jù)集的缺陷。

假設(shè)觀測數(shù)據(jù)集為X={xi|i=1,2,…,n},數(shù)據(jù)點xi和xj之間距離為dij=dist(xi,xj),該距離可以采用歐氏距離、馬氏距離、漢明距離和曼哈頓距離等。對于數(shù)據(jù)集X中的任何點xi,定義它的兩個基本屬性:局部密度ρi以及它與最近的高密度數(shù)據(jù)點的距離δi。ρi的計算方式有兩種:Cut-off kernel和Gaussian kernel,其中Cut-off kernel方式的計算公式為[17]

而χ函數(shù)定義為

參數(shù)dc＞0為截斷距離,需要用戶提前設(shè)置,且CFSFDP算法對截斷距離不敏感。由定義式(7)可以看出,局部密度ρi是數(shù)據(jù)集X中與xi之間距離小于dc的數(shù)據(jù)點的個數(shù)。

局部密度ρi的Gaussian kernel方式的計算式為[17]

比較定義式(7)和式(9)可知,Cut-off kernel方式的計算結(jié)果為離散值,Gaussian kernel方式的計算結(jié)果為連續(xù)值。因此,在Gaussian kernel的結(jié)果中,發(fā)生不同數(shù)據(jù)點具有相同局部密度值的概率會更小。

設(shè)qi(i=1,2,…,n)是ρi(i=1,2,…,n)的一個降序排列的序列,即ρq1≥ρq2

≥ρq3

≥ … ≥ρqn,則點qi與高密度數(shù)據(jù)點的距離定義為[17]

顯然,δi表示數(shù)據(jù)集X中任一點xi和所有局部密度大于它密度的點之間的最小距離;當(dāng)xi的局部密度是最大時,δi是其他所有聚類中最大的一個。局部密度最大的點一定也是一個數(shù)據(jù)簇的聚類中心。

對數(shù)據(jù)集X中每個點xi,計算出二元對(ρi,δi),以ρi為橫坐標(biāo),δi為縱坐標(biāo)畫出(ρi,δi)的決策圖,在該圖中,同時具有較大ρi和δi值的數(shù)據(jù)點會脫穎而出,這些點就可以看作是聚類中心;而ρi值很小且δi值很大的數(shù)據(jù)點就是離群點。顯然,利用決策圖確定聚類中心屬于定性分析而非定量分析,需要人為干預(yù)。為了減少人為因素的影響,定義一個綜合考慮ρi與δi值的量

ξi的值越大,則所對應(yīng)的數(shù)據(jù)點xi越有可能是聚類中心。

CFSFDP算法中截斷距離dc的確定方法如下。分配給每個點的平均鄰居數(shù)量約為數(shù)據(jù)點總數(shù)的1%～2%,所謂鄰居就是指某點在dc距離范圍內(nèi)的數(shù)據(jù)點。對數(shù)據(jù)集X中每個點,它與其他的n-1個點都有一個距離,總共有n(n-1)個距離。因為每個點對應(yīng)的距離都被計算了兩次,因此這些距離有一半是重復(fù)的。為此,把距離dij(i＜j)按升序排列為d1≤d2≤…≤dM。若取dc=dk(k=1,2,…,M),則在所有n(n-1)個距離中,小于dc的距離所占的比例約為t=k/M,即大約有(k/M)n(n-1)個距離小于dc。比值t就是選取參數(shù)dc的重要指標(biāo)。

參數(shù)dc值的選取決定著CFSFDP的聚類效果。若dc值過大,會使每個數(shù)據(jù)點的ρi值都很大,導(dǎo)致聚類的區(qū)分度不高;在極端情況下dc=dM,則所有的數(shù)據(jù)點都?xì)w屬于一個簇。若dc值太小,同一聚類中就可能被拆分成多個簇;在極端的情況下dc＜d1,則每一個數(shù)據(jù)點都單獨成為一個簇。選取dc值是依賴于具體問題的。通過采取比例值t來確定參數(shù)dc的策略,可降低對具體問題的依賴性。

從以上分析可知,DBSCAN是單獨以密度為指標(biāo)選擇聚類的類別,而CFSFDP綜合考慮了密度和距離兩個指標(biāo)作為選擇聚類的類別。因此,CFSFDP能夠區(qū)分出兩個密度值很接近的不同的兩個聚類。在利用聚類分析估計混合矩陣過程中,本文首先把計算得到的ξi(i=1,2,…,n)值按降序排列,然后從大到小截取由算法DBSCAN得到聚類中心數(shù)量相對應(yīng)的ξi值的數(shù)據(jù)點作為聚類中心。這意味著,把DBSCAN獲得的聚類數(shù)量作為CFSFDP的輸入?yún)?shù),通過CFSFDP進(jìn)一步搜索密度峰值以修正DBSCAN的聚類中心位置;同時DBSCAN也彌補了CFSFDP需要人為干預(yù)的不足,使兩種算法揚長補短,得到最佳的組合。

本文方法的基本流程如圖1所示。

圖1 本文提出方法的基本流程圖Fig.1 Basic flow chart of the proposed method

在利用聚類方法獲得混合矩陣的估計之后,采用最短路徑方法[3]就可以容易地分離(恢復(fù))出信源。

3 仿真結(jié)果與分析

為了驗證本文方法的有效性,在混合矩陣估計和信源分離兩個方面進(jìn)行仿真測試。特別地,把K-means算法[11,12]、基于層次的聚類(Hierarchical clustering,HC)算法[19]與本文算法在混合矩陣估計的性能方面進(jìn)行比較。仿真測試中,為了獲得盡可能公平的結(jié)果,所有算法的仿真環(huán)境都是同樣的。

3.1 仿真環(huán)境與性能指標(biāo)

仿真的PC平臺:Intel(R)Celeron(R)1007U-1.5 GHz的CPU,4 GB內(nèi)存,操作系統(tǒng)Windows 10,所有仿真都是運行在MATLAB 9(R2016a)上。用于測試的信源為SixFlutes數(shù)據(jù)集[3]中的長笛演奏音樂信號,其采樣率為44.1 kHz,信號樣本長度為216=65536。源信號向量記作s(t)=[s1(t),s2(t),s3(t),s4(t),s5(t),s6(t)]T,6路信源的時域波形如圖2所示。

圖26路音樂源信號的時域波形Fig.2 Time domain waveforms of six music source signals

為了測試算法對混合矩陣的估計精度,采用角度偏差[11]作為技術(shù)指標(biāo)

式中：a為原始混合矩陣A的列向量,b為估計出混合矩陣B的列向量,〈a,b〉表示向量a與b的內(nèi)積。如果角度偏差d(a,b)的值越小,說明混合矩陣的估計精度越高。另外,采用均方誤差(Mean square error,MSE)作為測度混合矩陣估計的另一個技術(shù)指標(biāo)，表達(dá)式為

式中：ak為原始混合矩陣A的第k個列向量,bk為估計出混合矩陣B的第k個列向量,ak和bk都是歸一化的向量。ak和bk的方向越接近,其|akTbk|的值越接近于1,說明混合矩陣估計精度越高。

在混合矩陣被估計出來之后,利用最短路徑法即可分離(恢復(fù))源信號。為了度量原始的信源與估計的信源之間的相似性,采用相關(guān)系數(shù)[11]作為性能指標(biāo)，表達(dá)式為

式中：si(t)為時域中第i個信源,rj(t)為時域中第j個恢復(fù)的信源,T為時域中信號的樣本數(shù)。如果恢復(fù)的信源與原始的信源越相似,其相關(guān)系數(shù)ρii就越接近于1或-1,而ρij(i≠j)越接近于0。

在下面的仿真中,首先在時域中把6個信源隨機地混合成3路觀測信號;然后利用STFT把時域信號變換到時頻域中進(jìn)行混合矩陣的估計,在進(jìn)行STFT操作時,利用Hanning窗截斷來實現(xiàn)對信號的分幀,每一幀的長度為L=8192,兩個連續(xù)幀的重疊率為60%;最后將恢復(fù)出的源信號通過逆STFT再變換到時域中,在時域中進(jìn)行相關(guān)系數(shù)的計算。這里對STFT的參數(shù)選擇說明如下。對于音頻信號的處理,一般來說是需要進(jìn)行分幀的;為了避免信號所含信息的損失，要求連續(xù)兩幀之間的樣本要重疊。由于每個信號的樣本長度為65536,要將信號分成整數(shù)幀(一般為2的整數(shù)次冪,本文取23=8),則可得每幀長度為L=65536/8=8192。對于連續(xù)兩幀之間重疊的樣本數(shù),在實際應(yīng)用中,一般取d=round(0.15×L×4),其中round(x)為對數(shù)據(jù)x取整操作[11]。這樣就可得到重疊的樣本數(shù)為d=round(0.6×L)=4915,即連續(xù)兩幀的重疊率為60%。利用窗函數(shù)截斷來實現(xiàn)音頻信號的分幀則是信號處理最重用的方法,在文獻(xiàn)[11]中對Barlett,Blackman,Flot top,Hanning,Hamming,Kaiser,Tukey,Welh,Boxcar等常見的窗函數(shù)的特點及適用的信號類型進(jìn)行了詳細(xì)的分析。Hanning窗函數(shù)能夠在較好幅度精度的情況下,提供良好的頻率分辨率和泄露保護(hù)[11];另外,在對音頻信號處理中,Hanning窗還具有非常低的頻率混疊的優(yōu)勢[11]。因此,本文采用Hanning窗實現(xiàn)對音頻信號的分幀處理。

3.2 結(jié)果及分析

仿真中,混合矩陣是利用MATLAB命令A(yù)=rand(3,6)隨機產(chǎn)生，即

3路觀測信號x(t)=[x1(t),x2(t),x3(t)]T=As(t)是由6個信源s(t)混合而成。x(t)的時域波形如圖3所示。

圖33路觀測信號的時域波形Fig.3 Time domain waveforms of three observed signals

圖4給出了x(t)的時域散點圖。所謂散點圖是指信號的數(shù)據(jù)點在直角坐標(biāo)系上的分布圖,其表示的是因變量隨自變量而變化的大致趨勢。從圖4可知,由于時域的數(shù)據(jù)點密集地分布在一起,沒有顯示出稀疏信號的線性聚類特性,因而無法分辨信源的數(shù)量和方向。為此,需要對時域信號x(t)進(jìn)行STFT,得到時頻域中的觀測信號X(t,k)=[X1(t,k),X2(t,k),X3(t,k)]T。在時頻域中信號X(t,k)的散點圖如圖5所示。

從圖5可以看出,時頻域中稀疏信號的線性聚類特性得到了明顯的增強。但在幾條由數(shù)據(jù)點組成的直線之間仍然分布著很多數(shù)據(jù)

點,這就造成了線性聚類的直線數(shù)量和方向都具有一定的不確定性。為解決這個問題,本文在時頻域中對X(t,k)進(jìn)一步采用SSP檢測,得到數(shù)據(jù)集Xssp,其散點圖如圖6所示。從圖6可看出,Xssp的散點圖明確地給出了6條直線,即反映了源信號的數(shù)目為6個。

為了應(yīng)用密度基的聚類算法對觀測數(shù)據(jù)進(jìn)行分析,本文通過鏡像映射方式(歸一化處理)把線性聚類的數(shù)據(jù)Xssp轉(zhuǎn)變?yōu)橹旅芫垲惖臄?shù)據(jù)X*(k)。觀測信號X*(k)的散點圖如圖7所示。

從圖7可看出,X*(k)在上半個單位球上形成了密集的數(shù)據(jù)堆,信源的數(shù)量由數(shù)據(jù)堆的個數(shù)給定。盡管在密集的6堆數(shù)據(jù)之外也分布著某些數(shù)據(jù)點,但通過聚類算法的不斷搜索即可實現(xiàn)這些數(shù)據(jù)的歸類。

在利用數(shù)據(jù)集X*(k)對混合矩陣的估計過程中,各種算法的主要參數(shù)設(shè)置如下。對于K-means,聚類數(shù)量事先給定為K=6,初始聚類中心位置隨機地生成,算法最大迭代次數(shù)為100,停止規(guī)則為產(chǎn)生的分配不再變化;對于HC算法,初始的每個數(shù)據(jù)點為一類,任意兩點間采用歐氏距離,通過合并距離最小的兩個類來發(fā)現(xiàn)數(shù)據(jù)簇,停止規(guī)則為合并后剩余類數(shù)量為6個。對于本文方法,DBSCAN算法的類簇鄰域esp=0.04,密度閾值MinPts=10,把DBSCAN獲得聚類數(shù)量作為CFSFDP算法的輸入?yún)?shù),CFSFDP算法中數(shù)據(jù)點之間采用歐氏距離,利用CFSFDP搜索密度峰值對聚類中心位置進(jìn)行修正。這里對DBSCAN參數(shù)的選擇說明如下。通過繪制觀測數(shù)據(jù)X*(k)的K-距離曲線,找出該曲線的拐點位置,可得到對應(yīng)的esp=0.04。而參數(shù)MinPts的選取原則為MinPts≥D+1,這里D=6,則MinPts≥7。在數(shù)據(jù)X*(k)的三維散點圖中,密集數(shù)據(jù)堆外還分布著一些數(shù)據(jù)點,經(jīng)計算可得這些點與最鄰近數(shù)據(jù)堆的空間距離約為2個單位長度,即MinPts＞7+2,于是取MinPts=10。

圖43路觀測信號的時域散點圖Fig.4 Time domain scatter plot of three observed signals

圖53路觀測信號的時頻域散點圖Fig.5 Time-frequency domain scatter plot of three observed signals

經(jīng)過本文方法、K-means算法和HC算法估計出的混合矩陣(記作B1,B2,B3)分別為

把原始混合矩陣A和各種方法估計得到的混合矩陣Bk(k=1,2,3)中各列向量代入到角度偏差d(a,b)和均方誤差MSE的定義式中,計算得到結(jié)果如表1所示。

由表1可得出以下結(jié)論:(1)K-means算法對混合矩陣第3個列向量估計的角度偏差到達(dá)了71.6920,這是不能接受的結(jié)果;同樣地,K-means算法估計的均方誤差也比較大。(2)HC聚類算法的估計性能與K-均值算法基本相當(dāng),而它估計的均方誤差是3種算法中最大的。(3)本文方法不但能夠克服K-均值算法要求預(yù)先知道源信號數(shù)目的缺陷,而且估計的性能指標(biāo)(無論是角度偏差還是均方誤差)是3種方法中最好的,這說明本文方法在估計欠定混合矩陣方面是非常有效的。

圖6 經(jīng)SSP檢測后觀測 信號的散點圖Fig.6 Scatter plot of observed signals after SSP detection

本文方法在獲得混合矩陣之后,利用最短路徑法[3]對源信號進(jìn)行恢復(fù),并計算出原始信號與恢復(fù)信號的相關(guān)系數(shù)分別為:ρ11=0.9498,ρ22=0.9546,ρ33=0.9333,ρ44=0.9888,ρ55=0.9681,ρ66=0.965 2。各信號的相關(guān)系數(shù)均大于0.93,且6個信號的平均相關(guān)系數(shù)為0.96,這說明本文方法的估計量具有很好的一致性。

在以上的仿真中,僅僅做了一次實驗。為了進(jìn)一步測試各種算法對欠定混合矩陣的平均估計性能,在下面的仿真中,對類似的實驗反復(fù)進(jìn)行了15次。

在每次實驗中,首先利用MATLAB函數(shù)rand(3,6)隨機地生成一個3×6的混合矩陣A從而得到時域觀測信號x(t)=As(t);然后通過STFT把時域信號轉(zhuǎn)換成時頻域的信號X(t,k),并對其進(jìn)行SSP檢測和歸一化處理,得到觀測信號X*(k);分別利用K-means,HC和本文方法對X*(k)進(jìn)行聚類分析從而估計出欠定的混合矩陣;最后計算出估計的角度偏差和均方誤差的指標(biāo)值,并把15次實驗的性能指標(biāo)(角度偏差和均方誤差)取平均值,得到結(jié)果如表2所示。

圖7 經(jīng)過歸一化處理后觀測信號的散點圖Fig.7 Scatter plot of observed signals after normalization

表13種方法的角度偏差和均方誤差值Tab.1 Angular deviation and MSE values of three methods

表23種方法的平均角度偏差和平均的均方誤差值Tab.2 Average angular deviation and average MSE values of three methods

從表2可看出:K-means算法對于混合矩陣的第3，4，5列向量估計的平均角度偏差較大(大于1),其中第4列的平均角度偏差高達(dá)13.7576;而HC算法對第1，3，4，5列向量估計的平均角度偏差都很大,其中第5列的平均角度偏差竟然為18.0434,而且HC的平均MSE是3種算法中最大的。這說明利用K-means和HC算法對欠定混合矩陣估計的誤差是很大的,因此會直接造成信源分離的精度下降。本文方法估計的混合矩陣列向量的平均角度偏差和平均MSE值都是這3種方法中最小的,除了第4列向量的角度偏差小于0.4之外,其他列向量的角度偏差都在0.07之下,這證明了本文方法對于欠定混合矩陣的估計精度很高。

在重復(fù)進(jìn)行的15次實驗中,對3種方法估計的混合矩陣,利用最短路徑法分離(恢復(fù))出信源,并對3種不同的方法，分別計算出每個源信號與對應(yīng)的恢復(fù)信號的相似度。圖8給出了3種方法每個信號在15次實驗中相關(guān)系數(shù)的變化曲線。

圖8 每次實驗中3種方法估計的信號相關(guān)系數(shù)Fig.8 Signal correlative coefficients estimated by three methods in each experiment

從圖8可以看出,對于全部的6個信源(sig1,sig2,sig3,sig4,sig5,sig6),由K-means算法進(jìn)行UBSS,所獲得信號相關(guān)系數(shù)的最小值為0.6517,由HC算法獲得相關(guān)系數(shù)的最小值為0.8407,而由本文方法獲得相關(guān)系數(shù)的最小值為0.9015。因此,用K-means和HC算法恢復(fù)的信源與原始信源的相似度較差,也就是說,本文方法具有更高的信源分離(恢復(fù))精度。同時還計算出了本文方法對每個信號在15次實驗中的平均相關(guān)系數(shù),它們的值分別是:0.9517,0.9508,0.9505,0.9616,0.9740和0.9673,把這6個信號的平均相關(guān)系數(shù)再取數(shù)學(xué)期望,則得到全部源信號的平均相關(guān)系數(shù)為0.9593,這是一個相當(dāng)好的UBSS性能指標(biāo)。

從以上的實驗結(jié)果可以看出，本文提出的算法具有較好的欠定混合矩陣的估計性能。而算法在對噪聲、信源相關(guān)性及欠定程度的適應(yīng)性方面,具有以下的特性。

(1)在實際應(yīng)用中,利用傳感器采集信源的過程,不可避免地會使觀測信號中包含有噪聲。在UBSS中,噪聲分為傳感器噪聲和信號源噪聲兩類,一般的BSS模型考慮的是傳感器噪聲[20]。對于本文算法,也進(jìn)行了含噪聲模型的實驗,即在采集信號中加入高斯噪聲,并通過降噪處理后再進(jìn)行UBSS。在具有噪聲的觀測信號的信噪比(Signal to noise ratio,SNR)分別為20,15和10 dB的情況下,本文算法能有效地估計出混合矩陣并實現(xiàn)信源的分離。這就是說,本文算法對噪聲的容忍度可達(dá)SNR=10 dB的惡劣環(huán)境。

(2)盲信源分離的主要方法是獨立分量分析(Independent component analysis,ICA)[20],即假設(shè)信源是相互獨立并且服從非高斯分布。然而,在對信源采集中,每個傳感器獲得的觀測信號與鄰近傳感器信號不可避免是互相關(guān)的,這種相關(guān)性也可能發(fā)生在觀測信號與非鄰近傳感器之間。對于具有互相關(guān)的信源,傳統(tǒng)BSS/ICA利用信源高階或二階統(tǒng)計性質(zhì)的方法將失效。為此,相關(guān)分量分析(Dependent component analysis,DCA)得到了廣泛的應(yīng)用。實際上,本文算法就是一種DCA方法,它是通過開發(fā)信源的稀疏性和非負(fù)性特征來實現(xiàn)對混合矩陣的估計;本文的單源點檢測就是對信源稀疏性的增強;而對數(shù)據(jù)的歸一化處理就是開發(fā)信源的非負(fù)特性。因此,本文算法對信源相關(guān)性具有較高的容忍度。

(3)對于欠定盲信源分離問題,傳感器數(shù)量決定了BSS的欠定程度。本文算法的仿真實驗中,在信源數(shù)量為6的情況下,分別取傳感器數(shù)量為4,3,2,1進(jìn)行了相應(yīng)的實驗。在傳感器為4,3,2這3種環(huán)境下算法對混合矩陣的估計均取得很好的性能。而當(dāng)傳感器數(shù)量為1時,即單通道UBSS,本文算法對混合矩陣的估計性能急劇下降,恢復(fù)的信源產(chǎn)生了較大的誤差。因此,本文算法對欠定程度的容忍度是要求最少有2個傳感器。

(4)盲源分離中的信源一般是非高斯分布。在實際應(yīng)用中,所遇到的非高斯信源,根據(jù)信號的峭度值可分為亞高斯信源和超高斯信源。本文討論的音頻信號屬于典型的超高斯信源,而一些圖像信號可能屬于亞高斯信源。本文的BSS算法既可以分離超高斯信源,也可以分離亞高斯信源。對于屬于亞高斯的圖像信號盲分離應(yīng)用,在文獻(xiàn)[21,22]中有介紹。

4 結(jié)束語

本文基于稀疏信號的UBSS問題,首先研究了增強信號稀疏性的方法,利用STFT把觀測信號從時域變換到時頻域,并采用SSP檢測剔除多源點數(shù)據(jù),通過鏡像映射把稀疏信號的直線聚類轉(zhuǎn)變成致密聚類;然后,采用聚類分析的方法在密集的數(shù)據(jù)堆中尋找代表其方向的關(guān)鍵數(shù)據(jù)。為了克服K-means算法需要預(yù)先設(shè)置聚類數(shù)目的缺點,本文采用DBSCAN算法實現(xiàn)對數(shù)據(jù)的自動分類從而得到源信號數(shù)目以及初始的聚類中心;在此基礎(chǔ)上,利用CFSFDP算法對DBSCAN獲得的聚類中心進(jìn)一步進(jìn)行修正,以提高混合矩陣的估計精度。由于把DBSCAN的聚類數(shù)量作為CFSFDP的輸入?yún)?shù),也克服了CFSFDP需要人為干預(yù)的缺陷。