何選森,徐 麗,夏 娟

(1.廣州商學(xué)院信息技術(shù)與工程學(xué)院，廣東 廣州 511363；2.湖南大學(xué)信息科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙 410082)

1 引言

在無(wú)噪聲情況下，盲信源分離BSS(Blind Source Separation)[1]的線性時(shí)域模型為x(t)=As(t)，其中s(t)=[s1(t),s2(t),…,sm(t)]T是信源向量，x(t)=[x1(t),x2(t),…,xn(t)]T是觀測(cè)信號(hào)向量，A∈Rn×m是體現(xiàn)傳感器采集信道特性的混合矩陣。一般地，信源與傳感器數(shù)量不相等，即m≠n。在BSS中，信源是不可觀測(cè)的，即向量s(t)所包含的信源數(shù)未知，而傳感器數(shù)量取決于信源的采集條件。在實(shí)際中，為節(jié)約成本，希望傳感器數(shù)盡可能地少，這就出現(xiàn)了n

很多自然信號(hào)(如音頻、圖像和視頻等)都是稀疏的。所謂稀疏信號(hào)是指非零元素很少，大部分元素都是零值的信號(hào)。實(shí)際應(yīng)用中，充分利用信號(hào)本身稀疏性的方法是解決UBSS問(wèn)題的首要方法[4]。對(duì)于充分稀疏的信號(hào)來(lái)說(shuō)，它的數(shù)據(jù)分布具有線性聚類的特性[5]。因此，聚類方法就成為了分析稀疏信號(hào)的重要工具。常用的聚類方法有K-均值(K-means)[6]、霍夫(Hough)變換[7]和層次聚類(Hierarchical Clustering)[8]等，其中K-means在UBSS中應(yīng)用廣泛。然而，K-means具有2個(gè)明顯的缺陷：(1)需要預(yù)先給定數(shù)據(jù)的聚類數(shù)量；(2)聚類結(jié)果對(duì)算法的隨機(jī)初始劃分非常敏感[9]。除了使用單獨(dú)的聚類方法外，把2種不同的聚類算法進(jìn)行組合，在估計(jì)欠定混合矩陣中也取得了很好的效果。例如，利用霍夫變換對(duì)K-means的聚類結(jié)果進(jìn)行修正，從而提高欠定混合信道的估計(jì)精度[10]；采取具有噪聲的基于密度空間聚類DBSCAN(Density-Based Spatial Clustering of Applications with Noise)算法[11]確定數(shù)據(jù)的聚類數(shù)量和聚類中心，再利用Hough變換對(duì)DBSCAN的結(jié)果進(jìn)行修正[12]，實(shí)現(xiàn)了對(duì)欠定系統(tǒng)混合矩陣的有效估計(jì)；在DBSCAN聚類形成密度基數(shù)據(jù)的基礎(chǔ)上，采用快速搜索與尋找密度峰值聚類CFSFDP(Clustering by Fast Search and Find of Density Peaks)[13]進(jìn)一步計(jì)算出對(duì)應(yīng)的密度峰值，以解決欠定混合矩陣的估計(jì)問(wèn)題[14]。

正是基于組合聚類的思想，為充分體現(xiàn)信號(hào)的稀疏性，本文提出在頻域中估計(jì)UBSS混合矩陣的分析方法。首先，利用傅里葉變換將時(shí)域中的觀測(cè)信號(hào)轉(zhuǎn)變成頻域中的稀疏信號(hào)，并通過(guò)數(shù)據(jù)的歸一化將稀疏信號(hào)的線性聚類轉(zhuǎn)變成致密聚類；然后，利用相似性(姻親)傳播AP(Affinity Propagation)聚類[15]搜索觀測(cè)數(shù)據(jù)的聚類數(shù)量和對(duì)應(yīng)的聚類中心，并把AP聚類的結(jié)果作為K-means算法的初始參數(shù)，進(jìn)一步修正數(shù)據(jù)的聚類中心，以提高欠定混合矩陣的估計(jì)精度；最后，利用線性規(guī)劃(最短路徑)法恢復(fù)信源。

本文其它部分組織如下：第2節(jié)介紹觀測(cè)信號(hào)的預(yù)處理，包括傅里葉變換以及數(shù)據(jù)的歸一化；第3節(jié)致力于稀疏信號(hào)的聚類分析，利用AP聚類實(shí)現(xiàn)對(duì)信源數(shù)目的自動(dòng)搜索，并采用K-means進(jìn)一步估計(jì)欠定混合矩陣的各個(gè)列向量；第4節(jié)介紹利用最短路徑法分離信源；第5節(jié)采用音頻信號(hào)進(jìn)行仿真，以驗(yàn)證算法有效性；第6節(jié)為結(jié)束語(yǔ)。

2 觀測(cè)數(shù)據(jù)的預(yù)處理

在時(shí)域線性UBSS模型x(t)=As(t)中，混合矩陣A的元素aij(1≤i≤n,1≤j≤m)為常數(shù)，將A分解為列向量:

A=[a1,a2,…,am]

(1)

其中ai=[a1i,a2i,…,ani]T,i=1,2,…,m。則UBSS模型又可以表示為:

x(t)=a1s1(t)+a2s2(t)+…+amsm(t)

(2)

一般地，信號(hào)在時(shí)域中的稀疏性體現(xiàn)得并不充分[5]，因此通常在頻域中對(duì)信號(hào)進(jìn)行聚類分析。對(duì)模型x(t)=As(t)，在2邊分別進(jìn)行短時(shí)傅里葉變換STFT(Short Time Fourier Transform):

(3)

其中,X(k)=[X1(k),X2(k),…,Xn(k)]T和S(k)=[S1(k),S2(k),…,Sm(k)]T分別為x(t)∈Rn和s(t)∈Rm的STFT系數(shù)向量，變量k是頻率抽樣點(diǎn)，而ai是矩陣A的第i個(gè)列向量。假設(shè)信源S(k)中所有信號(hào)都是完全稀疏的，即對(duì)于任一個(gè)頻域的樣本點(diǎn)k，僅有一個(gè)信源Si(k) (i=1,2,…,m)是非零值，其余信源Sj(k) (j≠i)均為零值，即:

X(k)=aiSi(k),i=1,2,…,m

(4)

或者

(5)

式(5)給出的是直線方程，且每個(gè)信源對(duì)應(yīng)于一條直線，其方向向量ai=[a1i,a2i,…,ani]T就是混合矩陣A的第i個(gè)列向量。對(duì)完全稀疏的信號(hào)，其數(shù)據(jù)具有線性聚類(Linear Clustering)的特性。在UBSS模型中，可以把對(duì)混合矩陣的估計(jì)轉(zhuǎn)變成對(duì)聚類形成的數(shù)據(jù)直線方向向量的估計(jì)。

實(shí)際上，信號(hào)的稀疏性在頻域中并未得到充分體現(xiàn)。即在每個(gè)離散頻率點(diǎn)k上，可能會(huì)有2個(gè)或2個(gè)以上信號(hào)具有非零值，使得信號(hào)的數(shù)據(jù)點(diǎn)不能形成理想的直線聚類特性。故采用單源點(diǎn)SSP(Single-Source-Point)檢測(cè)[14,16,17]增強(qiáng)稀疏信號(hào)的線性聚類特性。

一般地，一條直線的方向可以由2個(gè)方向向量來(lái)描述。例如，在三維空間中，一條經(jīng)過(guò)原點(diǎn)的直線可以由方向向量(1,1,1)和另一個(gè)方向向量(-1,-1,-1)來(lái)描述。為了把一條直線用唯一的方向向量來(lái)描述，可以采用一種稱為鏡像映射(Mirroring Mapping)[4]的方式，將負(fù)方向的向量(-1,-1,-1)映射到對(duì)應(yīng)正方向的向量(1,1,1)上，即把平面或空間上的一條直線映射成對(duì)應(yīng)的正方向上的一條射線。在信號(hào)處理領(lǐng)域，上述鏡像映射的過(guò)程實(shí)際上就是對(duì)信號(hào)歸一化。頻域中歸一化的觀測(cè)信號(hào)為:

(6)

其中,X(k)表示頻譜值，‖·‖表示歐幾里得(Euclidean)范數(shù)。顯然，觀測(cè)數(shù)據(jù)的歸一化是把負(fù)頻率方向的數(shù)據(jù)點(diǎn)映射到對(duì)應(yīng)正頻率方向的數(shù)據(jù)點(diǎn)上，形成了在上(正)半個(gè)單位圓(或單位球)上密集的數(shù)據(jù)堆。每一條聚類形成的數(shù)據(jù)直線都被映射成在上半個(gè)單位圓(球)上唯一的一堆數(shù)據(jù)。即鏡像映射是把稀疏信號(hào)的線性聚類轉(zhuǎn)變成觀測(cè)數(shù)據(jù)在正方向單位圓(球)上的致密聚類(Compact Clustering)[4]。

通過(guò)觀測(cè)數(shù)據(jù)歸一化，對(duì)欠定混合矩陣的估計(jì)就轉(zhuǎn)變成頻域中對(duì)上半個(gè)單位圓(球)觀測(cè)數(shù)據(jù)的聚類分析。由致密聚類特性[9]，在每一個(gè)密集的數(shù)據(jù)堆中，都存在一個(gè)關(guān)鍵的數(shù)據(jù)(聚類中心)，該數(shù)據(jù)的方向就對(duì)應(yīng)于線性聚類的直線方向。因此，通過(guò)搜索每一堆數(shù)據(jù)的聚類中心，就可以估計(jì)欠定混合矩陣的列向量。

3 聚類分析法搜索關(guān)鍵數(shù)據(jù)

聚類分析能夠?qū)⒕哂幸欢ㄌ刭|(zhì)的數(shù)據(jù)劃分為一族，而屬于不同族的數(shù)據(jù)則易于區(qū)分。平方誤差和準(zhǔn)則(Sum-of-Squared-Error Criterion)是聚類分析中使用最廣泛的準(zhǔn)則函數(shù)，利用該準(zhǔn)則能夠?qū)χ旅軘?shù)據(jù)集進(jìn)行很好的分類[9]。

在眾多聚類算法中，K-means采用平方誤差和準(zhǔn)則，以迭代優(yōu)化步驟搜索數(shù)據(jù)最優(yōu)劃分[9]。由于它易于實(shí)現(xiàn)，被認(rèn)為是經(jīng)典的聚類分析法，特別是對(duì)于密集數(shù)據(jù)堆和球(超球)的數(shù)據(jù)分布形狀，K-means能獲得理想的聚類效果[18]。盡管如此，K-means的弱點(diǎn)同樣也很明顯：使用算法的前提是要給定數(shù)據(jù)的聚類數(shù)量；聚類效果對(duì)算法隨機(jī)產(chǎn)生的初始劃分非常敏感，即要求有一個(gè)合理的初始化。顯然，UBSS無(wú)法滿足這2個(gè)要求。因此，本文采用AP聚類作為K-means的輸入?yún)?shù)，以彌補(bǔ)K-means需要人為干預(yù)的缺陷，從而提高欠定混合矩陣的估計(jì)精度。

基于相似性(Similarity)測(cè)度的數(shù)據(jù)聚類在數(shù)據(jù)分析和工程系統(tǒng)中具有重要的作用[15]。一種常見(jiàn)的方法是使用數(shù)據(jù)來(lái)學(xué)習(xí)一組(聚類)中心，使得數(shù)據(jù)點(diǎn)與其最近的中心之間的平方誤差和最小。當(dāng)聚類中心是從實(shí)際的數(shù)據(jù)點(diǎn)獲得，則稱其為范例(Exemplars)[15]。AP聚類最明顯的優(yōu)勢(shì)是將所有的數(shù)據(jù)點(diǎn)都等同地考慮為潛在的范例[15]，這樣就有效地避免了隨機(jī)選擇初始聚類中心的局限性。

對(duì)于數(shù)據(jù)點(diǎn)xi和xk，定義它們之間的相似性測(cè)度為:

sim(i,k)=-‖xi-xk‖2

(7)

為了尋找合適的范例，AP聚類在搜索過(guò)程中不斷積累2個(gè)證據(jù)：責(zé)任性(Responsibility)和可用性(Availability)。責(zé)任性r(i,k)表示從點(diǎn)xi的角度來(lái)看，點(diǎn)xk作為點(diǎn)xi的范例的適合程度；而可用性av(i,k)則表示從候選樣本點(diǎn)xk的角度來(lái)看，點(diǎn)xi選擇點(diǎn)xk作為它的范例的適合程度。顯然，r(:,k)+av(:,k)的數(shù)值越大，數(shù)據(jù)點(diǎn)xk作為最終聚類中心的概率就越大。根據(jù)不斷累計(jì)的證據(jù)，AP通過(guò)一個(gè)迭代過(guò)程搜尋數(shù)據(jù)的聚類，直到高質(zhì)量的范例(聚類中心)出現(xiàn)。證據(jù)信息r(i,k)和av(i,k)在AP聚類過(guò)程中可以被看作是對(duì)數(shù)概率比(Log-probability Ratios)。

當(dāng)AP迭代開(kāi)始時(shí)，將av(i,k)值初始化為零：av(i,k)=0；而責(zé)任性r(i,k)的計(jì)算規(guī)則如式(8)所示:

(8)

在第1次迭代中，由于可用性av(i,k)=0，責(zé)任性r(i,k)為sim(i,k)和xi與其它候選范例之間相似性最大值的差。而在以后的迭代過(guò)程中，當(dāng)一些數(shù)據(jù)點(diǎn)被有效地分配給其它的候選范例時(shí)，根據(jù)式(9)的更新規(guī)則，它們的可用性值av(:,:)將降到零以下(即為負(fù)值)。負(fù)的可用性數(shù)值將降低式(8)中一些相似性的有效值，從而在競(jìng)爭(zhēng)中消除相應(yīng)的候選范例。對(duì)k=i的情況，自責(zé)任性r(k,k)設(shè)置為點(diǎn)xk作為范例的偏好值減去點(diǎn)xi與其它候選范例之間相似性中的最大值，它反映了點(diǎn)xk作為一個(gè)范例的證據(jù)積累，用于糾正不合適分配造成的錯(cuò)誤。

責(zé)任性r(i,k)更新規(guī)則(式(8))是讓所有候選范例都爭(zhēng)奪數(shù)據(jù)點(diǎn)的所屬權(quán)，而式(9)所示的可用性av(i,k)更新規(guī)則是從數(shù)據(jù)點(diǎn)收集證據(jù)，用于說(shuō)明每個(gè)候選范例是否會(huì)成為一個(gè)好的范例:

av(i,k)=min{0,r(k,k)+

(9)

即av(i,k)設(shè)置為自責(zé)任性r(k,k)加上候選范例xk從其它點(diǎn)獲得的正值的責(zé)任性之和。只有所得責(zé)任性的正值部分才被加到一起，這是因?yàn)閷?duì)一個(gè)好的范例，必須(用正的責(zé)任性)解釋它如何適合一些數(shù)據(jù)點(diǎn)，而不必(用負(fù)的責(zé)任性)解釋它如何不適合一些數(shù)據(jù)點(diǎn)。當(dāng)以xk為范例的其它點(diǎn)具有正值的責(zé)任性，而自責(zé)任性r(k,k)是負(fù)值時(shí)，點(diǎn)xk作為一個(gè)范例的可用性得到增強(qiáng)。為限制強(qiáng)的正值責(zé)任性的影響，式(9)中的總和被設(shè)置了閾值，即av(i,k)不能大于零。自可用性av(k,k)的更新規(guī)則為:

(10)

基于其它點(diǎn)發(fā)送給候選范例xk正值的責(zé)任性，式(10)反映了點(diǎn)xk作為一個(gè)范例的累計(jì)證據(jù)。

在AP聚類中還有2個(gè)重要的參數(shù)：相似性矩陣Sim={sim(i,k)}對(duì)角元素構(gòu)成的偏好p(k)=sim(k,k)以及抑制算法陷入振蕩的阻尼因子λ。由式(8)可得自責(zé)任性的表示式如式(11)所示:

根據(jù)地域特點(diǎn)與現(xiàn)有資源，打造品牌概念，為游客提供具有本地特色的獨(dú)創(chuàng)性商品服務(wù)，提高鄉(xiāng)村的吸引力和影響力，讓身處繁雜都市的游客從自發(fā)尋找心靈休憩之所，變成慕名前來(lái)游玩度假。此外，還應(yīng)提升鄉(xiāng)村旅游的文化價(jià)值，保護(hù)旅游重點(diǎn)區(qū)的“鄉(xiāng)村性”，使其免于在城市化進(jìn)程中丟失自身的民俗文化。

(11)

由式(9)和式(11)可知，p(k)越大，則r(k,k)和av(i,k)就越大，點(diǎn)xk作為最終聚類中心的概率就越大。p(k)值將影響到哪個(gè)以及多少個(gè)范例將在競(jìng)爭(zhēng)中勝出，最終成為聚類中心。這意味著，調(diào)整p(k)的值，可以增加或減少聚類數(shù)量。在實(shí)際應(yīng)用中，一個(gè)較好的選擇方案是將所有的p(k)值都設(shè)置為所有對(duì)應(yīng)的數(shù)據(jù)點(diǎn)之間的相似性數(shù)值的中值(Median)。

當(dāng)更新r(i,k)和av(i,k)的信息時(shí)，為避免出現(xiàn)數(shù)值振蕩(Numerical Oscillations)，對(duì)它們的迭代過(guò)程還需要利用阻尼因子進(jìn)行抑制。在每一個(gè)迭代步驟i中，r(i,k)和av(i,k)的值要與最近一次迭代(步驟i-1)中的值r(i-1,k)和av(i-1,k)一起進(jìn)行更新?？紤]阻尼因子λ的迭代規(guī)則為:

av(i,k)=(1-λ)av(i,k)+λav(i-1,k)

r(i,k)=(1-λ)r(i,k)+λr(i-1,k)

(12)

其中λ∈[0,1]，其默認(rèn)值為0.5。λ還可以改進(jìn)算法的收斂性能。當(dāng)由于振蕩(或識(shí)別出的范例處于周期變化中)使得AP聚類無(wú)法收斂時(shí)，通過(guò)增加λ的值可消除振蕩。

AP更新規(guī)則僅需要執(zhí)行簡(jiǎn)單的局部計(jì)算，即只需在一對(duì)已知相似性測(cè)度的點(diǎn)之間交換信息。在AP過(guò)程的任何點(diǎn)(迭代步驟i)，用av(i,k)和r(i,k)的組合即可識(shí)別出范例。對(duì)于點(diǎn)xi，找到使av(i,k)+r(i,k)最大化的點(diǎn)xk的值，當(dāng)k=i時(shí)辨識(shí)出點(diǎn)xi作為范例，k≠i時(shí)，辨識(shí)出點(diǎn)xk是點(diǎn)xi的范例。即AP聚類總能找出數(shù)據(jù)的聚類數(shù)量和對(duì)應(yīng)的聚類中心。顯然，把AP聚類的結(jié)果作為K-means初始條件，不僅克服了K-means的缺陷，同時(shí)還保持了經(jīng)典K-means的優(yōu)勢(shì)。

4 線性規(guī)劃方法分離信源

在獲得混合矩陣估計(jì)之后，“兩步法”的第2步是分離信源。考慮到環(huán)境噪聲，UBSS的時(shí)域模型為x(t)=As(t)+v(t)，其中v(t)∈Rn是加性高斯噪聲向量。對(duì)于充分稀疏的信源，其觀測(cè)信號(hào)在頻域中每個(gè)離散的頻點(diǎn)上，僅有一個(gè)分量具有非零值，其余分量均為零值，即給定的非零分量為稀疏分布。所謂稀疏分布是指在自變量為零時(shí)有尖峰，且隨著自變量增大尖峰拖了一個(gè)很長(zhǎng)的尾部，例如雙邊指數(shù)(拉普拉斯)分布。稀疏性等價(jià)于非高斯性。

對(duì)于稀疏信號(hào)，頻域中的盲源分離則是通過(guò)如式(13)所示的優(yōu)化問(wèn)題[3]實(shí)現(xiàn)的:

(13)

其中,σ2是噪聲的方差。式(13)的第1項(xiàng)為均方重建誤差的和(高斯噪聲的對(duì)數(shù)似然性)，而第2項(xiàng)為對(duì)非稀疏性(假設(shè)獨(dú)立的拉普拉斯信源或超高斯信源)的處罰。

式(13)是個(gè)多變量?jī)?yōu)化問(wèn)題，它需要同時(shí)估計(jì)混合矩陣A和信源S。由于利用組合聚類方法(見(jiàn)第3節(jié))已經(jīng)估計(jì)出混合矩陣A，因此，優(yōu)化問(wèn)題(13)簡(jiǎn)化成式(14)所示的形式:

(14)

在BSS問(wèn)題中，含噪聲BSS(Noisy-BSS)和UBSS模型都是BSS的難點(diǎn)[2]。一般地，BSS模型僅考慮無(wú)噪情況(in the absence of noise)[3,4,12]，因此，式(14)的優(yōu)化問(wèn)題可進(jìn)一步簡(jiǎn)化為式(15)[3]:

k=1,2,…,F

(15)

其中F為離散頻率范圍。式(15)是只有1個(gè)變量的優(yōu)化問(wèn)題,即線性規(guī)劃問(wèn)題，求解較容易。

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

為驗(yàn)證本文算法的有效性，考慮在惡劣環(huán)境中僅用2個(gè)傳感器采集多信源進(jìn)行仿真。仿真中，對(duì)本文算法(Proposed)、K-means算法[6,19]、層次聚類HC(Hierarchical Clustering)算法[8,20]、模糊C均值FCM(Fuzzy C-means)算法[21]、混合矩陣方向估計(jì)DEMIX(Direction Estimation of MIXing matrix)算法[22]的性能進(jìn)行比較。為得到公平的結(jié)果，各算法仿真環(huán)境都完全一樣。

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

仿真的PC機(jī)配置為：Intel(R) Celeron(R) 1007U-1.5 GHz的CPU，4 GB內(nèi)存，操作系統(tǒng)為Windows 10，所有仿真都在Matlab 9 (R2016a)上運(yùn)行。測(cè)試的音頻信源為SixFlutes (長(zhǎng)笛演奏音樂(lè))數(shù)據(jù)集[3]，信源向量為s(t)=[s1(t),s2(t),s3(t),s4(t),s5(t),s6(t)]T，信號(hào)采樣率為44.1 kHz，每個(gè)信號(hào)長(zhǎng)度為216=65536個(gè)樣本。每次仿真從s(t)中隨機(jī)抽取不同數(shù)量和不同位置的信號(hào)。

為了度量各算法對(duì)混合矩陣的估計(jì)精度，采用如式(16)所示的角度偏差(Angular Deviation)[5]指標(biāo)進(jìn)行度量:

(16)

其中,a是原始混合矩陣A的某個(gè)列向量，b是估計(jì)的混合矩陣B對(duì)應(yīng)的列向量，而〈a,b〉表示向量a與b的內(nèi)積。角度偏差d(a,b)的值越小，說(shuō)明混合矩陣的估計(jì)精度越高。

在獲得混合矩陣的估計(jì)之后，利用最短路徑法[3]即可分離出信源。為度量原始信源與估計(jì)的信源之間的相似性，采用如式(17)所示的相關(guān)系數(shù)(Correlative Coefficients)[14]指標(biāo)進(jìn)行度量:

(17)

其中，si(t)是時(shí)域中第i個(gè)信源，rj(t)是時(shí)域中第j個(gè)被分離的信源，T是時(shí)域中信號(hào)的樣本數(shù)。若分離的信源與原始信源越相似，其相關(guān)系數(shù)ρii的絕對(duì)值就越接近于1，而ρij(i≠j)越接近于0。

5.2 仿真過(guò)程

首先，在時(shí)域中進(jìn)行信號(hào)混合。利用Matlab隨機(jī)地產(chǎn)生混合矩陣A，把音頻信源s(t)轉(zhuǎn)變成觀測(cè)信號(hào)x(t)=As(t)。

對(duì)于時(shí)域中的觀測(cè)信號(hào)x(t)，利用Hanning窗截?cái)鄟?lái)實(shí)現(xiàn)音頻信號(hào)分幀。一般地，需要將一個(gè)完整信號(hào)分成整數(shù)幀，幀的數(shù)量取2的整數(shù)次冪，本文取23=8。由于每個(gè)信號(hào)的長(zhǎng)度為65 536，可得每個(gè)音頻幀長(zhǎng)為L(zhǎng)=65536/8=8192。為避免音頻信號(hào)所含信息的損失，要求連續(xù)(相鄰)2幀之間的樣本要重疊。實(shí)際中，連續(xù)2個(gè)音頻幀之間重疊樣本數(shù)為d=round(0.15×L×4)，其中round(x)為對(duì)x取整。于是可得到d=4915，即連續(xù)2幀的重疊率為60%。

然后，在頻域中估計(jì)混合矩陣。利用STFT對(duì)時(shí)域信號(hào)x(t)的每一幀(采用同樣的Hanning窗)進(jìn)行變換，并且只取正頻率范圍內(nèi)頻譜(實(shí)部與虛部)的L個(gè)STFT系數(shù)，將連續(xù)(相鄰)的頻域幀銜接起來(lái)就形成一個(gè)信號(hào)，從而獲得頻域信號(hào)向量X(k)。為便于聚類分析，對(duì)X(k)進(jìn)行鏡像映射獲得歸一化的數(shù)據(jù)X*(k)。

對(duì)于歸一化信號(hào)X*(k)，本文算法首先利用AP聚類自動(dòng)搜索聚類數(shù)量和對(duì)應(yīng)的聚類中心，然后把該結(jié)果作為K-means的輸入?yún)?shù)，實(shí)現(xiàn)對(duì)聚類中心(矩陣列向量)的精確估計(jì)。

最后，在得到混合矩陣估計(jì)之后，利用線性規(guī)劃方法分離頻域的信源。對(duì)每個(gè)分離出的頻域信源，STFT的系數(shù)向量被分解成幀，其實(shí)部和虛部分量又被重新組合成復(fù)系數(shù)，頻譜也拓展到負(fù)頻率范圍；對(duì)該頻譜利用逆STFT(使用同樣的Hanning窗)再變換到時(shí)域中。將時(shí)域信號(hào)每一幀乘以對(duì)應(yīng)的逆Hanning窗，在連續(xù)2幀重疊部分的2邊各移除50%的樣本，僅保留幀的中間部分，把這些時(shí)域幀按先后順序銜接起來(lái)，得到分離的信源，并計(jì)算相關(guān)系數(shù)。

5.3 仿真結(jié)果

5.3.1 2個(gè)傳感器采集3個(gè)信源

由Matlab命令randi([1,6],1,3)隨機(jī)抽取3個(gè)信源s(t)=[s1(t),s2(t),s6(t)]T，由命令A(yù)=2*rand(2,3)-1隨機(jī)地生成代表白噪聲信道的矩陣為A2×3=[-0.0798,0.3703,0.1458; 0.8008,-0.1234,0.9501]，由x(t)=As(t)=[x1(t),x2(t)]T得到時(shí)域觀測(cè)信號(hào)，其時(shí)域散點(diǎn)圖如圖1所示。

從圖1知，時(shí)域中的數(shù)據(jù)點(diǎn)在坐標(biāo)原點(diǎn)附近形成密集的數(shù)據(jù)族，沒(méi)有體現(xiàn)信號(hào)的稀疏性，無(wú)法進(jìn)行聚類分析。為此，將時(shí)域觀測(cè)信號(hào)x(t)利用STFT變換成頻域信號(hào)X(k)=[X1(k),X2(k)]T。頻域中信號(hào)X(k)的散點(diǎn)圖如圖2所示。

Figure 1 Time-domain scatter plot of the observed signals圖1 觀測(cè)信號(hào)的時(shí)域散點(diǎn)圖

從圖2可看出，雖然頻域中觀測(cè)信號(hào)X(k)的稀疏性得到了一定的體現(xiàn)，隱約有數(shù)據(jù)點(diǎn)形成直線，但線性聚類并不明顯。還需要對(duì)頻域信號(hào)采用單源點(diǎn)SSP檢測(cè)[14,16]技術(shù)進(jìn)行預(yù)處理，以提高信號(hào)的線性聚類特性。

Figure 2 Frequency-domain scatter plot of the signals圖2 觀測(cè)信號(hào)頻域散點(diǎn)圖

在頻域中，SSP是指僅有一個(gè)信源具有較大能量而其余信源能量很小以致可被忽略的頻率點(diǎn)。然而在很多頻點(diǎn)上，2個(gè)或2個(gè)以上信源都具有較大的能量值，這些頻點(diǎn)就稱為多源點(diǎn)MSP(Multi-Source Point)。正是MSP影響了信號(hào)的線性聚類。SSP檢測(cè)就是利用信號(hào)頻譜實(shí)部與虛部夾角的余弦[14,16]作為相似性度量，通過(guò)合適的閾值，檢測(cè)并刪除MSP，從而體現(xiàn)由SSP數(shù)據(jù)點(diǎn)形成的線性聚類特性。經(jīng)過(guò)SSP檢測(cè)后信號(hào)X(k)的散點(diǎn)圖如圖3所示。

Figure 3 Frequency-domain scatter plot after SSP圖3 SSP檢測(cè)后的頻域散點(diǎn)圖

從圖3可看出，經(jīng)過(guò)SSP檢測(cè)，觀測(cè)數(shù)據(jù)在頻域散點(diǎn)圖上形成了明確的3條直線，即直觀地給出了不可觀察的信源數(shù)量；欠定混合矩陣的列向量則體現(xiàn)為數(shù)據(jù)直線的方向向量。然而，基于平方誤差和準(zhǔn)則的聚類方法都是把觀測(cè)數(shù)據(jù)劃分為不同形狀和不同大小的數(shù)據(jù)族，這就需要用鏡像映射將線性聚類轉(zhuǎn)變成致密聚類，形成在正半個(gè)單位圓上數(shù)據(jù)的密集聚類。經(jīng)過(guò)歸一化后的數(shù)據(jù)散點(diǎn)圖如圖4所示。

Figure 4 Normalized scatter plot of the signals圖4 信號(hào)的歸一化散點(diǎn)圖

從圖4可以看出，在正的半個(gè)單位圓上，歸一化數(shù)據(jù)形成密集的數(shù)據(jù)堆，而堆的數(shù)量就是線性聚類形成的直線數(shù)量(信源數(shù)量)。另外，圖4中數(shù)據(jù)并沒(méi)有被完全聚到一起，還有部分?jǐn)?shù)據(jù)需要進(jìn)一步歸類到對(duì)應(yīng)數(shù)據(jù)堆中，這正是各種聚類方法所要執(zhí)行的聚類任務(wù)。

對(duì)于歸一化的數(shù)據(jù)，分別采用5種聚類算法搜索聚類的數(shù)量和相應(yīng)的聚類中心，從而估計(jì)欠定混合矩陣。各方法的參數(shù)設(shè)置如下：在本文算法中，當(dāng)AP聚類的阻尼系數(shù)取默認(rèn)值λ=0.5時(shí)，形成的聚類中心數(shù)值周而復(fù)始地重復(fù)振蕩，為此本文取λ=0.9消除其周期性。DEMIX算法采用觀測(cè)數(shù)據(jù)對(duì)混合矩陣列向量進(jìn)行計(jì)數(shù)和方向估計(jì)。K-means算法的初始聚類中心采用歐氏距離算術(shù)平均值法[19]，聚類重復(fù)次數(shù)為1，最大迭代次數(shù)為100，算法的停止規(guī)則是相鄰2次聚類中心的迭代值相同。HC算法采用文獻(xiàn)[20]中的“Algorithm 1”，該算法已給出了R代碼，仿真中將R改寫(xiě)成Matlab命令執(zhí)行。FCM采用文獻(xiàn)[21]中的“Algorithm 2”，其中范數(shù)為歐氏距離，指定聚類數(shù)量C，以樣本點(diǎn)到聚類中心隸屬度大小劃分?jǐn)?shù)據(jù)歸類。

基于生成的混合矩陣A2×3，5種算法估計(jì)得到的混合矩陣分別為:

A1=[-0.9457,0.9487,0.0598; 0.3252,-0.3162,0.9908]

A2=[-0.0991,0.9487,0.1514; 0.9951,-0.3163,0.9885]

A3=[-0.1003,0.9477,0.1513; 0.9949,-0.3192,0.9885]

A4=[-0.9288,0.9520,0.0532; 0.3361,-0.3151,0.9911]

A5=[-0.9459,0.9491,0.0601; 0.3308,-0.3149,0.9881]

其中A1,A2,…,A5分別表示K-means、DEMIX、Proposed、HC和FCM算法估計(jì)的矩陣。

由此可得5種算法在估計(jì)混合矩陣各個(gè)列向量時(shí)，所引起的角度偏差值如表1所示。

Table 1 Angular deviation of 3 sources separatedfrom 2 observed signals

由表1可以看出，DEMIX算法的性能是最優(yōu)的，這是由于DEMIX是為估計(jì)混合矩陣列向量方向而專門設(shè)計(jì)的算法[22]，它的角度偏差指標(biāo)本就應(yīng)該是最好；其它幾種都是通用聚類算法，只是將它們用于估計(jì)混合矩陣。K-means對(duì)第2個(gè)列向量估計(jì)的角度偏差值是這5種算法中最小的，說(shuō)明K-means確實(shí)是經(jīng)典的聚類算法，然而K-means、HC、FCM 3種算法對(duì)其它列向量估計(jì)的角度偏差卻很大，尤其是K-means和HC算法對(duì)第1個(gè)列向量估計(jì)的角度偏差都超過(guò)62。盡管本文算法對(duì)混合矩陣各列向量估計(jì)的角度偏差指標(biāo)不是5種算法中最優(yōu)的，但也僅次于DEMIX的結(jié)果，并且其值都小于0.2，明顯優(yōu)于其它3種聚類算法。

利用各種算法估計(jì)出的混合矩陣B，采取最短路徑法在頻域分離信源，通過(guò)逆STFT轉(zhuǎn)換成時(shí)域信號(hào)，從而計(jì)算出原始信源與估計(jì)信源之間的相關(guān)系數(shù)，其結(jié)果如表2所示。

Table 2 Correlative coefficients of 3 sources separatedfrom 2 observed signals

從表2可以看出，對(duì)于恢復(fù)信源與原始信源的相似性指標(biāo)，本文算法與DEMIX算法的估計(jì)性能幾乎是相同的。盡管K-means對(duì)信源6的分離性能是這5種算法中最優(yōu)的，但它與HC、FCM對(duì)其它信源的分離效果較差。

對(duì)照表1和表2的結(jié)果可知，混合矩陣列向量的角度偏差指標(biāo)與分離信源的相似性指標(biāo)并不完全一致。這是因?yàn)榛旌暇仃囀窃陬l域中估計(jì)的，而信源的相關(guān)系數(shù)是在時(shí)域中計(jì)算的。從仿真過(guò)程可知，頻域信源正頻率的STFT系數(shù)向量采用Hanning窗分幀，相鄰2幀的樣本要重疊，將STFT實(shí)部和虛部重新組合成復(fù)系數(shù)且把頻譜拓展到負(fù)頻率的范圍；然后再利用相同長(zhǎng)度的Hanning窗對(duì)該頻譜進(jìn)行逆STFT變換到時(shí)域中；采用時(shí)域Hanning窗進(jìn)行分幀，并在相鄰2幀重疊部分的2邊各移除50%的樣本，僅保留幀中間部分，從而構(gòu)成時(shí)域信源。顯然在STFT和逆STFT過(guò)程中的數(shù)據(jù)分幀、樣本重疊、樣本移除操作，會(huì)造成部分樣本發(fā)生變化，使得2種指標(biāo)產(chǎn)生不一致現(xiàn)象。

為了比較5種算法對(duì)混合矩陣估計(jì)和信源分離的總體性能，對(duì)以上仿真過(guò)程重復(fù)進(jìn)行了10次，每次分別記錄各算法對(duì)混合矩陣各列向量估計(jì)的角度偏差值，其結(jié)果如圖5所示。同時(shí)，記錄各算法對(duì)各個(gè)分離出信源與原始信源的相似性指標(biāo)，其結(jié)果如圖6所示。

Figure 5 Angular deviation of 3 sources separated from 2 observed signals圖5 由2路觀測(cè)信號(hào)分離3路信源的角度偏差

從圖5可看出，在10次仿真中，DEMIX算法與本文算法的角度偏差值基本保持一致，都是小于1的，即這2種算法對(duì)混合矩陣的估計(jì)性能是優(yōu)秀的。而K-means角度偏差最大值竟達(dá)到149.881 9,HC的最大角度偏差值為66,FCM的最大角度偏差值為22.421 7,說(shuō)明這3種算法對(duì)欠定混合矩陣的估計(jì)性能較差。

Figure 6 Correlative coefficients of 3 sources separated from 2 observed signals圖6 由2路觀測(cè)信號(hào)分離3路信源的相關(guān)系數(shù)

從圖6及記錄的數(shù)據(jù)可知，本文算法在10次仿真中最小的相關(guān)系數(shù)為0.956 3,DEMIX的最小相關(guān)系數(shù)為0.952 2,這2種算法的絕大多數(shù)相關(guān)系數(shù)都在0.99以上，說(shuō)明它們具有很好的UBSS性能。而K-means在10次仿真中最小相關(guān)系數(shù)為0.493 0，HC最小相關(guān)系數(shù)為0.499 3，F(xiàn)CM最小相關(guān)系數(shù)為0.559 9,這3種算法大約有一半的相關(guān)系數(shù)小于0.9，說(shuō)明這3種算法恢復(fù)的信源相對(duì)于原始信源產(chǎn)生了較大失真。

由以上的仿真結(jié)果可知，DEMIX與本文算法在估計(jì)混合矩陣和分離信源方面的性能大致相同。其它3種算法的性能屬于同一級(jí)別，并且K-means是這3種算法中性能最差的。因此，在下面的仿真中，只對(duì)DEMIX、K-means和本文算法的性能進(jìn)行比較。

5.3.2 2個(gè)傳感器采集5個(gè)信源

5個(gè)信源s(t)=[s1(t),s2(t),s3(t),s4(t),s5(t)]T，由Matlab命令A(yù)=2*rand(2,5)-1隨機(jī)生成混合矩陣為A2×5=[0.8394,-0.7207,0.1236,-0.9114,0.3763;-0.5956,0.7679,-0.3332,0.3283,-0.6618]，時(shí)域觀測(cè)信號(hào)為x(t)=As(t)=[x1(t),x2(t)]T。然后，x(t)通過(guò)STFT變換成頻域信號(hào)X(k)=[X1(k),X2(k)]T，頻譜X(k)的散點(diǎn)圖如圖7所示。歸一化數(shù)據(jù)X*(k)的散點(diǎn)圖如圖8所示。

Figure 7 Frequency-domain scatter plot of 5 sources and 2 observed signals圖7 2路觀測(cè)信號(hào)(5信源)的頻域散點(diǎn)圖

Figure 8 Normalized frequency-domain scatter plot圖8 2路觀測(cè)信號(hào)(5信源)的歸一化散點(diǎn)圖

對(duì)于圖8的歸一化觀測(cè)數(shù)據(jù)，利用DEMIX、K-means和本文算法分別對(duì)混合矩陣各列向量進(jìn)行估計(jì)。這3種算法的參數(shù)設(shè)置如前所述，而且仿真過(guò)程與上一個(gè)仿真類似。

為了全面地比較這3種算法對(duì)混合矩陣的估計(jì)性能，同樣的仿真過(guò)程重復(fù)進(jìn)行10次，計(jì)算并記錄10次仿真中混合矩陣列向量估計(jì)的平均角度偏差值，其結(jié)果如表3所示。

Table 3 Average angular deviation of 5 sources separated from 2 observed signals

從表3可知，DEMIX和本文算法的平均角度偏差值都小于1，但DEMIX要優(yōu)于本文算法。而K-means的平均角度偏差值都大于1，且有3個(gè)列向量的平均角度偏差值都在25以上。

為了更直觀地觀察每次仿真中各算法估計(jì)混合矩陣列向量的角度偏差值的變化情況，圖9給出了10次仿真的結(jié)果。

從圖9及記錄的數(shù)據(jù)可知，DEMIX的最大角度偏差值為1.163 3,本文算法的最大偏差值為2.547 1；K-means的最大偏差值為132.428 6。顯然，K-means估計(jì)混合矩陣的誤差較大。

同樣地，在10次仿真中，3種算法分離5個(gè)信源的平均相關(guān)系數(shù)值如表4所示。

從表4可知，DEMIX和本文算法恢復(fù)信源的平均相關(guān)系數(shù)幾乎是相同的，而且只有2個(gè)信號(hào)的相關(guān)系數(shù)小于0.9。而K-means對(duì)4個(gè)信號(hào)的平均相關(guān)系數(shù)都小于0.9，其中最小的相關(guān)系數(shù)僅有0.733 1。

Figure 9 Angular deviation of separating 5 sources from 2 observed signals圖9 2路觀測(cè)信號(hào)(5信源)的角度偏差

算法信源1信源2信源3信源4信源5K-means0.931 40.817 90.862 60.733 10.846 1DEMIX0.934 20.840 30.901 20.853 60.914 6Proposed0.938 10.843 60.903 10.874 30.914 9

圖10給出了3種算法在每次仿真中恢復(fù)信源與原始信源相關(guān)系數(shù)的變化波形。

Figure 10 Correlative coefficients of 5 sources separated from 2 observed signals圖10 由2路觀測(cè)信號(hào)分離5路信源的相關(guān)系數(shù)

從記錄的數(shù)據(jù)可知，在5個(gè)信號(hào)的10次仿真中，本文算法的相關(guān)系數(shù)最小值為0.754 1，最大值為0.964 5；DEMIX的最小值為0.650 7，最大值為0.965 4；K-means的最小值為0.445 3，最大值為0.967 1。從圖10的相關(guān)系數(shù)波形也可以看出，DEMIX與本文算法分離信源的相似性指標(biāo)基本一致，而K-means的相關(guān)系數(shù)值波動(dòng)較大，說(shuō)明它對(duì)信源分離的性能不穩(wěn)定，產(chǎn)生了較大的分離誤差。

綜合考慮以上2個(gè)仿真，可以得出如下結(jié)論：(1)在固定傳感器數(shù)量(2個(gè))的情況下，隨著信源數(shù)目的增加，各種聚類算法對(duì)混合矩陣的估計(jì)性能以及對(duì)信源的分離性能都會(huì)下降。這是由于在UBSS系統(tǒng)中，信源之間會(huì)產(chǎn)生干擾，信源數(shù)量越多，采集傳感器遭受的干擾就越強(qiáng)，因而造成算法對(duì)UBSS的估計(jì)性能下降。(2)由于DEMIX是專門用于估計(jì)混合矩陣列向量方向的算法，因此在角度偏差的指標(biāo)上，它相對(duì)其它聚類算法具有一定的優(yōu)勢(shì)。盡管如此，本文算法仍能和DEMIX保持基本一致的(性能指標(biāo))水平。(3)單獨(dú)的聚類算法，如K-means、HC和FCM，對(duì)混合矩陣估計(jì)和信源分離的效果都不理想，產(chǎn)生了較大的估計(jì)誤差。基于此，本文將2種聚類算法相結(jié)合，首先利用AP聚類辨識(shí)出信源(類)數(shù)量和聚類中心，再利用K-means進(jìn)行聚類中心(關(guān)鍵數(shù)據(jù))的定位修正，從而提高混合矩陣以及信源的估計(jì)精度。

6 結(jié)束語(yǔ)

針對(duì)稀疏信源的UBSS問(wèn)題，本文提出了在頻域中估計(jì)混合矩陣的聚類分析算法。首先利用STFT將時(shí)域觀測(cè)信號(hào)轉(zhuǎn)變成頻域中的稀疏信號(hào)，以增強(qiáng)信源的線性聚類特性。然后利用鏡像映射將線性聚類轉(zhuǎn)變成致密聚類，以實(shí)現(xiàn)對(duì)觀測(cè)信號(hào)的歸一化。在此基礎(chǔ)上，利用聚類方法搜索數(shù)據(jù)族的數(shù)量和對(duì)應(yīng)的關(guān)鍵數(shù)據(jù)。為了克服K-means算法本身的缺陷，本文利用基于相似性傳播的AP聚類自動(dòng)搜索出不可觀測(cè)的信源數(shù)量和相對(duì)應(yīng)的數(shù)據(jù)聚類中心，并將AP聚類得到的結(jié)果作為K-means的輸入?yún)?shù)，對(duì)UBSS系統(tǒng)混合矩陣的列向量進(jìn)行精確的估計(jì)。把AP和K-means相結(jié)合的優(yōu)勢(shì)，不僅可以提高對(duì)UBSS系統(tǒng)的估計(jì)精度，而且也彌補(bǔ)了K-means算法需要人為干預(yù)(即需要指定聚類數(shù)量以及對(duì)初始化敏感)的不足。