基于源數(shù)估計(jì)的無(wú)約束欠定盲源分離算法

付永慶，郭慧，，蘇東林，劉焱

(1.哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江哈爾濱150001;2.北京航空航天大學(xué)電磁兼容技術(shù)研究所，北京100191)

盲源分離(blind signal separation，BSS)起源于雞尾酒問題［1］，可在源信號(hào)和傳輸信道參數(shù)均未知的情況下僅根據(jù)源信號(hào)的統(tǒng)計(jì)獨(dú)立性和觀測(cè)信號(hào)來(lái)恢復(fù)源信號(hào)，被廣泛應(yīng)用于生物醫(yī)學(xué)［2］、地質(zhì)信號(hào)處理和分析［3］、圖像處理［4］和移動(dòng)通信［5］等領(lǐng)域。根據(jù)觀測(cè)信號(hào)和源信號(hào)的數(shù)目關(guān)系，將盲源分離問題分為正定、超定和欠定3種情況。上述情況都以信源數(shù)目的正確估計(jì)為前提實(shí)現(xiàn)信號(hào)的分離。目前，在信源數(shù)目估計(jì)領(lǐng)域中，以前兩種情況下的盲源分離為主，基于高階累積量擴(kuò)展了很多成熟的算法，而欠定盲源分離(undetermined blind signal separation，UBSS)方面的研究尚少。目前欠定盲源分離中信源數(shù)目的估計(jì)大多基于稀疏混合信號(hào)的線性聚類特點(diǎn)，如:2001年Bofill利用勢(shì)函數(shù)的局部最大值的個(gè)數(shù)對(duì)信源數(shù)目進(jìn)行估計(jì)［6］;2008年譚北海等直接利用概率統(tǒng)計(jì)知識(shí)，通過劃分區(qū)域并統(tǒng)計(jì)觀測(cè)點(diǎn)在不同區(qū)域的概率分布實(shí)現(xiàn)信源數(shù)目估計(jì)［7］;2009年張燁等提出采用拉普拉斯勢(shì)函數(shù)判斷其局部最大值達(dá)到源數(shù)估計(jì)的目的［8］。但這些方法都存在一些問題，如算法的抗噪聲性和信號(hào)的稀疏性易受異常值的影響從而造成誤判。針對(duì)上述問題，提出一種Hough加窗法，利用信號(hào)的稀疏特性，先確定經(jīng)平滑后的聚類區(qū)域，減小異常值的影響，再尋找區(qū)域最大值，避免陷入局部最大，在此基礎(chǔ)上對(duì)混合矩陣進(jìn)行估計(jì)。針對(duì)信號(hào)分離問題，提出一種無(wú)約束分離算法，采用“內(nèi)點(diǎn)法”選擇合適的混合矩陣列向量作為初始迭代值，修正了現(xiàn)有算法對(duì)混合矩陣列向量的選取標(biāo)準(zhǔn)。

1 稀疏欠定盲源分離模型

通常欠定盲源分離的線性瞬時(shí)混合模型寫作:

式中:X(t)為m個(gè)觀測(cè)信號(hào)矢量，A為m×n維的混合矩陣，S(t)為n個(gè)源信號(hào)矢量，t=1，2，…，T為觀測(cè)點(diǎn)時(shí)刻，當(dāng)m＜n時(shí)即為欠定盲源分離。

為了實(shí)現(xiàn)欠定情況下的盲源分離，必須對(duì)源信號(hào)提出進(jìn)一步的約束條件，假設(shè)源信號(hào)是稀疏信號(hào)，即源信號(hào)在絕大多數(shù)采樣點(diǎn)的取值為零或接近于零，只有少數(shù)采樣點(diǎn)的取值遠(yuǎn)離零［9］。假設(shè)在采樣時(shí)刻t，只有s1(t)起主導(dǎo)作用，則有

觀測(cè)信號(hào)可以看作:在以x1(t)…xm(t)為坐標(biāo)軸的m維空間中，所有以源信號(hào)s1(t)為主導(dǎo)的采樣時(shí)刻確定的一條直線，其斜率取決于混合矩陣A的第一個(gè)列向量 [a11，…，am1]T。如果混合矩陣A的任意m×m的子矩陣都滿足可逆性，則每個(gè)源信號(hào)可確定一條直線。因此，源信號(hào)的數(shù)目估計(jì)就轉(zhuǎn)化為對(duì)觀測(cè)空間中直線的條數(shù)估計(jì)。

實(shí)際應(yīng)用中，很多信號(hào)具有稀疏特性，或者可以通過適當(dāng)變換(如短時(shí)傅里葉變換，STFT)［10］、Gabor變換［11］、小波變換［12］等)使得信號(hào)在變換域中滿足稀疏特性。

2 Hough加窗法

2.1 Hough 變換

Hough變換原用于圖像的匹配，將原圖中給定形狀的所有點(diǎn)都集中到變換空間中的某些位置，形成峰值點(diǎn)?；谠撍枷?，把觀測(cè)信號(hào)空間中直線的檢測(cè)問題轉(zhuǎn)化為尋找變換空間中峰值點(diǎn)的問題。

觀測(cè)數(shù)據(jù)可看作m維空間中經(jīng)過零點(diǎn)的一條直線［13］，據(jù)此構(gòu)造轉(zhuǎn)換關(guān)系，將觀測(cè)信號(hào)空間中的直線轉(zhuǎn)換成變換空間中的角度變換量［14］:

式中:k=1，…，π/h，表示高維直方圖的分區(qū)數(shù)目，h是Hough變換的量化步長(zhǎng)，k越大，即h越小，變換量的分類精度越高;在變換量的每一個(gè)分類中，i=1，…，m－1，包含了變換量的所有維數(shù)，[·]表示比·小的最大整數(shù)。則落入高維直方圖某一區(qū)域的變換量數(shù)目為:

2.2 加窗法

得到變換量的直方圖后，若直接檢測(cè)峰值點(diǎn)容易陷入局部最大值。為此對(duì)直方圖做加窗處理。由于m維的觀測(cè)信號(hào)經(jīng)Hough變換后變?yōu)閙－1維變換量，故對(duì)其加m－1維的超方體窗。加窗的目的是平滑直方圖中的異常值，得到聚類區(qū)域。加窗后重新計(jì)算變換量的分類統(tǒng)計(jì)數(shù)為

式中:j=0，…，[π/hd]－[L/d]+1表示重新計(jì)算后的變換量分類數(shù)，L表示超方體窗口的邊長(zhǎng)，d表示超方體窗口沿觀測(cè)空間每個(gè)坐標(biāo)軸方向的移動(dòng)步長(zhǎng)表示權(quán)重，由落入直方圖某一區(qū)域中的變換量數(shù)占變換量總數(shù)的比率獲得。

對(duì)新的直方圖搜索峰值即可得到源數(shù)的估計(jì)。峰值對(duì)應(yīng)的角度變換估計(jì)量與混合矩陣元素有式(7)的關(guān)系，由此得到混合矩陣的估計(jì)值:

式中:ami表示混合矩陣第m行i列的元素。

3 無(wú)約束分離算法

3.1 無(wú)約束問題

文獻(xiàn)［2］中分離混合信號(hào)有如下的約束條件:

由于式(1)可寫為

式中:Am稱作混合矩陣Am的滿秩因子，Amn－m稱作混合矩陣Am的余式因子。又有

式中:sign(·)表示取·的符號(hào)。因此，可以將式(8)帶約束的優(yōu)化問題轉(zhuǎn)化為上述無(wú)約束優(yōu)化問題，要滿足，令J(Sn－m)對(duì)Sn－m求導(dǎo)，即

3.2 初始迭代值選取

初始迭代值的選取直接影響算法的迭代次數(shù)和收斂速度，因此，需要合理選擇初始值?，F(xiàn)有算法通過選擇觀測(cè)點(diǎn)與混合矩陣列向量夾角最小的前m列作為分離矩陣，該方法并不適用于任何情況。

圖1 二維矢量分解示意圖Fig.1 Schematic diagram of 2D vector decomposition

如圖1所示，以兩個(gè)混合信號(hào)為例，在單位化對(duì)稱化后的散點(diǎn)圖中，混合信號(hào)x與列向量a1、a2的夾角最小，但是其分解矢量的長(zhǎng)度和遠(yuǎn)大于混合信號(hào)x在列向量a2、a3上的分解矢量。這是因?yàn)楫?dāng)兩個(gè)列向量分布于混合信號(hào)的同一側(cè)時(shí)，有一個(gè)分解矢量必在其中一列向量的反方向上，此時(shí)，分解矢量與觀測(cè)信號(hào)的夾角并不是最小。

根據(jù)散點(diǎn)圖分布采用“內(nèi)點(diǎn)法”選取初始迭代值。首先，令混合矩陣單位化和對(duì)稱化，且使混合矩陣第一行和混合信號(hào)第一行數(shù)據(jù)符號(hào)為正。然后，去處混合矩陣和混合信號(hào)第一行數(shù)據(jù)后繪制散點(diǎn)圖，以混合信號(hào)為原點(diǎn)，混合矩陣列向量分布在其四周。從n個(gè)點(diǎn)中選取m個(gè)點(diǎn)，這m個(gè)點(diǎn)圍成的多邊形應(yīng)能包含混合信號(hào)點(diǎn)。從滿足要求的組合中選擇與混合信號(hào)夾角和最小的作為滿秩因子Am，否則按照最小夾角法選取Am，其余列向量作為余式因子An－m，Sm則由式(10)得到，令Sn－m=0。

源信號(hào)越稀疏，J'(Sn－m)值越小，在迭代過程中，分離信號(hào)逐漸逼近源信號(hào)S，稀疏性增強(qiáng)，J'(Sn－m)逐漸減小，當(dāng)J'(S(k+1)n－m≥J'()時(shí)迭代終止，第k次迭代計(jì)算得到的即為分離信號(hào)。

4 仿真實(shí)驗(yàn)

4.1 算法性能評(píng)價(jià)標(biāo)準(zhǔn)

4.1.1 泛化交擾誤差

采用泛化交擾誤差(generalized crosstalking error，GCE)評(píng)價(jià)混合矩陣的估計(jì)精度。定義如下:

式中:Π表示所有n×n維可逆矩陣組成的集合，表示與一個(gè)尺度矩陣和置換矩陣的乘積，以消除的幅度不確定性和排序不確定性。當(dāng)且僅當(dāng)和A完全等價(jià)時(shí)，Err(A，)=0。

4.1.2 信干比

信號(hào)分離后第i個(gè)源信號(hào)si和分離信號(hào)的信干比(signal to interference radio，SIR)定義為:

盲源分離的幅度具有不確定性，在計(jì)算信干比前需先確定源信號(hào)si和分離信號(hào)的幅度因子λi:

通過式(15)使得si和的幅度盡量一致。在該度量下，信干比SIR越大，說(shuō)明分離信號(hào)越接近源信號(hào)。

4.1.3 互相關(guān)系數(shù)

互相關(guān)系數(shù)(cross correlation coefficient，CCC)描述兩個(gè)信號(hào)的相似性，定義為

式中:sj為源信號(hào)，為還原信號(hào)。在該度量下，ξij=1和0分別表示與sj完全相似和不相似。

在上述3個(gè)標(biāo)準(zhǔn)中，泛化交擾誤差和互相關(guān)系數(shù)都是與特定值做比較(GCE值與0比較，CCC值與1比較)，可以體現(xiàn)單獨(dú)算法的有效性，但因其計(jì)算值很小，在多算法性能比較時(shí)不利于突出算法的差異，因此采用信干比來(lái)衡量不同算法的性能。

4.2 實(shí)驗(yàn)及結(jié)果分析

4.2.1 實(shí)驗(yàn)一:驗(yàn)證Hough加窗法的性能

為了驗(yàn)證信號(hào)稀疏度對(duì)算法的影響，需構(gòu)建稀疏度可度量的稀疏信號(hào)。滿足廣義Gaussian分布的信號(hào)，其概率密度函數(shù)為:

據(jù)此構(gòu)建信噪比為40 dB，α＜1.5的稀疏度不同的廣義Gaussian分布的隨機(jī)信號(hào)和α=1，信噪比不同的信號(hào)，采樣值T=5 000，隨機(jī)選取2×6的混合矩陣得到混合信號(hào)，分別進(jìn)行100次蒙特卡洛運(yùn)算，分析算法的正確率，如圖2所示。

圖2 源數(shù)估計(jì)的正確性Fig.2 Accuracy of the source number estimation

Hough加窗法在稀疏性度量值為1.2及以上時(shí)，算法的正確率較低，在稀疏性度量值小于1.1和信噪比大于2 dB時(shí)算法正確率分別可達(dá)81%和85%以上，在稀疏性度量值小于0.8和信噪比大于15 dB時(shí)算法正確率達(dá)到100%?？梢娝惴▽?duì)信號(hào)的稀疏敏感性較低，且具有較好的抗噪聲性能，即使在低信噪比下也能以較高的正確率估計(jì)出信源數(shù)目。

源數(shù)估計(jì)的高正確率保證了算法的穩(wěn)定性。選取源數(shù)估計(jì)正確的情況估計(jì)混合矩陣，計(jì)算泛化交擾誤差取平均值，并與勢(shì)函數(shù)法、K均值法進(jìn)行比較，得到不同算法下噪聲和稀疏度對(duì)混合矩陣估計(jì)性能的影響圖，如圖3和圖4所示。

從兩幅圖中看出，無(wú)論是在低信噪比下還是低稀疏度下，Hough加窗法得到的混合矩陣的泛化交擾誤差都很小，趨近于0，尤其是在稀疏度為1.2以下和信噪比為8 dB以上時(shí)，混合矩陣估計(jì)精度更高，性能明顯優(yōu)于其他算法。勢(shì)函數(shù)法直接對(duì)累計(jì)量進(jìn)行峰值搜索，無(wú)法平滑奇異值對(duì)峰值檢測(cè)的影響，因此性能不如Hough加窗法，K均值法由于受初始點(diǎn)選取的影響較大，故估計(jì)性能最差。

圖3 噪聲對(duì)混合矩陣估計(jì)精度的影響Fig.3 The influence of noiseon the accuracy of the estimation of mixing matrix

圖4 稀疏度對(duì)混合矩陣估計(jì)精度的影響Fig.4 The influence of sparse on the accuracy of the estimation of mixing matrix

4.2.2 實(shí)驗(yàn)二:驗(yàn)證無(wú)約束分離算法的性能

在實(shí)驗(yàn)一得到的混合信號(hào)基礎(chǔ)上，假設(shè)混合矩陣已知，信噪比為40 dB，α=1的情況下，通過互相關(guān)系數(shù)和信干比兩個(gè)指標(biāo)來(lái)評(píng)價(jià)算法的性能，如圖5所示，并從分離信號(hào)的信干比、時(shí)間頻度、復(fù)雜度和運(yùn)行時(shí)間等方面，將本文算法與最短路徑法、最小l1范數(shù)法進(jìn)行比較，結(jié)果如表1所示。

圖5中信噪比為40 dB時(shí)，互相關(guān)系數(shù)平均值達(dá)到0.991 8，對(duì)應(yīng)的信干比平均值為169.65。整體來(lái)說(shuō)，在信噪比較低的情況下，互相關(guān)系數(shù)可達(dá)到0.74以上，在信噪比為15 dB以上時(shí)，互相關(guān)系數(shù)達(dá)到0.9以上。在2～40 dB的范圍內(nèi)，信干比的數(shù)值波動(dòng)在37 dB左右，約占最佳信干比值的21.9%。表明無(wú)約束分離算法自身具有分離準(zhǔn)確性較高，抗噪聲性能較好的特點(diǎn)。

表1給出了無(wú)約束分離算法與其他2個(gè)有約束算法的性能比較結(jié)果。其中時(shí)間頻度以m個(gè)混合信號(hào)，n個(gè)源信號(hào)，采樣值為T的情況進(jìn)行計(jì)算，k表示無(wú)約束分離算法的迭代次數(shù)。3種分離算法得到的源信號(hào)信干比值均一樣，說(shuō)明了本文算法的正確性。從時(shí)間頻度上看，最短路徑法值最小，最小l1范數(shù)法最大。本文算法與最小l1范數(shù)法相差一個(gè)系數(shù)，本文算法的k取值一般為1～3，最小l1范數(shù)法在本實(shí)驗(yàn)中取值為15，運(yùn)行時(shí)間也證明了上述分析。在運(yùn)算時(shí)間上，本文算法不及勢(shì)函數(shù)法，但適用范圍要大于勢(shì)函數(shù)法，因?yàn)閯?shì)函數(shù)法只適用于兩個(gè)混合信號(hào)的情況。因此，綜合運(yùn)行時(shí)間和適用性可以得出:無(wú)約束分離算法具有較好的性能。

表1 算法比較結(jié)果(α=1，SNR=40 dB)Table 1 The comparison results of different algorithms(α =1，SNR=40 dB)

圖5 混合分離的分離性能Fig.5 Performance of signal separation

為了驗(yàn)證“內(nèi)點(diǎn)法”修正了最小夾角法對(duì)混合矩陣列向量的選取標(biāo)準(zhǔn)，將無(wú)約束分離算法與基于最小夾角法的分離算法進(jìn)行比較，結(jié)果如圖6所示。

圖6 內(nèi)點(diǎn)法與最小夾角法的性能比較Fig.6 Performance comparison of signal separation between interiorpoint method and minimum angle method

在混合矩陣選取的過程中，不同的混合信號(hào)出現(xiàn)列向量散點(diǎn)包含混合信號(hào)散點(diǎn)的情況是不確定的，因此，針對(duì)不同的混合信號(hào)而言，內(nèi)點(diǎn)法較最小夾角法的性能優(yōu)勢(shì)表現(xiàn)的也不盡相同。在圖6中，內(nèi)點(diǎn)法的優(yōu)勢(shì)較最小夾角法比較突出，且在信噪比較高時(shí)，這種優(yōu)勢(shì)更明顯。

4 結(jié)論

1)針對(duì)現(xiàn)有混合矩陣估計(jì)算法易受噪聲干擾，對(duì)信號(hào)稀疏性要求較高等問題，提出一種Hough加窗算法，利用稀疏信號(hào)的特性，用加窗的方法平滑因噪聲或稀疏性差造成的異常值，對(duì)得到的聚類區(qū)域分別尋找最大值，達(dá)到全局最大值搜索的目的，峰值的個(gè)數(shù)即為信源的數(shù)目。

2)針對(duì)現(xiàn)有分離算法帶約束的優(yōu)化問題，提出一種無(wú)約束的分離方法，并采用“內(nèi)點(diǎn)法”來(lái)選擇合適的初始迭代值，修正了現(xiàn)有算法對(duì)混合矩陣列向量的選取標(biāo)準(zhǔn)。

仿真結(jié)果不但驗(yàn)證了算法的有效性，還表明了算法具有較好的抗噪聲性、較低的稀疏敏感性。若信號(hào)不滿足稀疏特性，可對(duì)其采用適當(dāng)?shù)淖儞Q使其在變換域中具有稀疏性，該算法在變換域中仍適用。

［1］COMON P.Independent component analysis:a new concept［J］.Signal Processing，Special Issue on Higher Order Statistics，1994，36(3):287-314.

［2］MAMMONE N，La FORESTA F，MORABITO F C.Automatic artifact rejection from multichannel scalp EEG by wavelet ICA［J］.IEEE Sensors Journal，2012，12(3):533-542.

［3］TAKAHATA A K，NADALIN E Z，F(xiàn)ERRARI R，et al.Unsupervised processing of geophysical signals:a review of some key aspects of blind deconvolution and blind source separation［J］.IEEE Signal Processing Magazine，2012，29(4):27-35.

［4］BADAWI W K M，CHIBELUSHI C C，PATWARY M N，et al.Specular-based illumination estimation using blind signal separation techniques［J］.IET Image Processing，2012，6(8):1181-1191.

［5］YANG Xiaomei.The application of blind signal separation in mobile communication［C］//2011 3rd International Conference on Computer Research and Development(ICCRD).Shanghai，China，2011:274-277.

［6］BOFILL P，ZIBULEVSKY M.Underdetermined blind source separation using sparse representations［J］.Signal Processing，2001，81(11):2353-2362.

［7］譚北海，謝勝利.基于源信號(hào)數(shù)目估計(jì)的欠定盲分離［J］.電子與信息學(xué)報(bào)，2008，30(4):863-867.TAN Beihai，XIE Shengli.Underdetermined blind separation based on source signals’number estimation［J］.Jour-nal of Electronics＆ Information Technology，2008，30(4):863-867.

［8］張燁，方勇.基于拉普拉斯勢(shì)函數(shù)的欠定盲分離中源數(shù)的估計(jì)［J］.信號(hào)處理，2009，25(11):1719-1725.ZHANG Ye，F(xiàn)ANG Yong.A new method to estimate the number of the sources for underdetermined blind separation based on lapulacial potential function［J］.Signal Processing，2009，25(11):1719-1725.

［9］邱天爽，畢曉輝.稀疏分量分析在欠定盲源分離問題中的研究進(jìn)展及應(yīng)用［J］.信號(hào)處理，2008，24(6):966-970.QIU Tianshuang，BI Xiaohui.Spares component analysis and application for underdetemined blind source separation［J］.Signal Processing，2008，24(6):966-970.

［10］BOFILL P，ZIBULEVSKY M.Blind separation of more sources than mixtures using sparsity of their short-time fourier transform［C］//Second International Workshop on Independent Component Analysis and Blind Signal Separation.Espoo，F(xiàn)inland，2000:87-92.

［11］CORDERO E，NICOLA F，RODINO L.Sparsity of gabor representation of schrodinger propagators［J］.Applied and Computational Harmonic Analysis，2009，26(3):357-370.

［12］MA Jianwei，PLONKA G，CHAURIS H.A new sparse representation of seismic data using adaptive easy-path wavelet transform［J］.IEEE Geoscience and Remote Sensing Letters，2010，7(3):540-544.

［13］HIROYUKI S，YUZO H.Blind source separation by a geometrical method［C］//Proceedings of the International Joint Conference on Neural Networks.Honolulu，USA，2002:1108-1114.

［14］付寧，彭喜元.K-Hough欠定盲信道估計(jì)算法［J］.電子測(cè)量與儀器學(xué)報(bào)，2008，22(5):63-67.FU Ning，PENG Xiyuan.K-Hough underdetermined blind mixing model recovery algorithm［J］.Journal of Electronic Measurement and Instrument，2008，22(5):63-67.

［15］何昭水，謝勝利，傅予力.信號(hào)的稀疏性分析［J］.自然科學(xué)進(jìn)展，2006，16(9):1167-1173.HE Shaoshui，XIE Shengli，F(xiàn)U Yuli.Sparse analysis of signal［J］.Progress in Natural Science，2006，16(9):1167-1173.