基于平行因子分析的欠定混合矩陣估計(jì)算法

王少波，郭 英，眭 萍，李紅光，楊 鑫

(空軍工程大學(xué)信息與導(dǎo)航學(xué)院, 陜西 西安 710077)

0 引言

盲源分離是指從觀測(cè)到的多源混合信號(hào)中分離并恢復(fù)出相對(duì)獨(dú)立的源信號(hào)的過程[1]。盲源分離是現(xiàn)代信號(hào)處理的重要前沿研究領(lǐng)域之一，已經(jīng)在通信、語音處理、機(jī)械故障分析、圖像處理、生物醫(yī)學(xué)、雷達(dá)及經(jīng)濟(jì)數(shù)據(jù)分析等領(lǐng)域得到廣泛應(yīng)用[2-4]。

盲源信號(hào)分離技術(shù)根據(jù)源信號(hào)數(shù)目和混合信號(hào)數(shù)目之間的關(guān)系可以分為：超定、正定和欠定。超定即觀測(cè)信號(hào)數(shù)目大于源信號(hào)數(shù)目，正定即觀測(cè)信號(hào)數(shù)目等于源信號(hào)數(shù)目，欠定即觀測(cè)信號(hào)數(shù)目小于源信號(hào)數(shù)目，其分別對(duì)應(yīng)著不同的解決方法，目前解決超定、正定盲源分離的方法[5]已經(jīng)相當(dāng)成熟，而欠定盲分離問題的解決方法還有待進(jìn)一步的研究?！皟刹椒ā笔墙陙斫鉀Q欠定盲源分離問題的一個(gè)有效方法，即首先估計(jì)混合矩陣,然后在混合矩陣已知的基礎(chǔ)上恢復(fù)源信號(hào)。估計(jì)混合矩陣的精度將直接影響分離信號(hào)的質(zhì)量，因此混合矩陣估計(jì)算法非常關(guān)鍵。

欠定條件下混合矩陣的估計(jì)，主要采用稀疏聚類方法[6-8]和張量分解法[9-13]。文獻(xiàn)[8]選擇時(shí)頻域做為稀疏變換域，通過計(jì)算接收混合信號(hào)時(shí)頻比的方差來檢測(cè)單源區(qū)域，再利用k-均值聚類來完成混合矩陣的估計(jì)?；谙∈杈垲惖幕旌暇仃嚬烙?jì)方法，計(jì)算相對(duì)簡(jiǎn)單，估計(jì)精度較高，但是需要滿足一定的稀疏條件，這限制了該方法的應(yīng)用范圍。文獻(xiàn)[9]將壓縮感知理論與平行因子分析相結(jié)合，實(shí)現(xiàn)了欠定盲源分離及DOA估計(jì)，但其只適用于具有特殊表示形式的源信號(hào)。文獻(xiàn)[10]和文獻(xiàn)[12]將平行因子分析應(yīng)用到盲源分離問題中，將混合矩陣的估計(jì)問題轉(zhuǎn)化為張量分解問題，該方法提高了估計(jì)的準(zhǔn)確性和魯棒性，但是進(jìn)行CP分解的ALS算法對(duì)初值敏感且容易陷入局部極值，運(yùn)算時(shí)間較長(zhǎng)。文獻(xiàn)[13]引入塔克分解先把張量壓縮為較低維的核張量，然后運(yùn)用交替最小二乘對(duì)核張量進(jìn)行標(biāo)準(zhǔn)分解，得到混合矩陣的估計(jì)，減少了運(yùn)算時(shí)間，但是初值敏感問題依然存在?；趶埩糠纸獾墓烙?jì)算法，克服了稀疏聚類算法稀疏性要求，通常需要利用信號(hào)的高階統(tǒng)計(jì)量構(gòu)造張量，將混合矩陣的估計(jì)問題轉(zhuǎn)化為張量分解問題，適用于信源相互獨(dú)立并且非高斯分布的情況，精度較高，但是算法的計(jì)算量較大，運(yùn)算時(shí)間長(zhǎng)，且對(duì)初值敏感，容易陷入局部收斂。目前研究一種有效的欠定條件下混合矩陣估計(jì)方法仍是一個(gè)技術(shù)難點(diǎn)。

針對(duì)現(xiàn)有算法在解決非稀疏信號(hào)的欠定混合矩陣估計(jì)中，存在的計(jì)算時(shí)間長(zhǎng)、初值敏感且容易陷入局部收斂的問題，提出了基于平行因子分析的欠定混合矩陣估計(jì)算法。

1 欠定盲源分離模型及張量理論基礎(chǔ)

1.1 欠定盲源分離模型

欠定盲源分離是指在觀測(cè)數(shù)目小于源信號(hào)數(shù)目的條件下，從觀測(cè)到的多源混合信號(hào)中分離并恢復(fù)出相對(duì)獨(dú)立的源信號(hào)的過程。圖1為欠定盲源分離系統(tǒng)原理圖。

設(shè)M個(gè)相互獨(dú)立且非高斯分布觀測(cè)信號(hào)是來自于N個(gè)源信號(hào)的線性混合，則欠定盲源分離信號(hào)模型描述為：

X(t)=AS(t)+N(t)

(1)

式(1)中，X(t)=[x1(t),x2(t),…,xM(t)]T表示接收的M個(gè)觀測(cè)信號(hào)，S(t)=[s1(t),s2(t),…,sN(t)]T表示輸入的N個(gè)未知源信號(hào)，M

圖1 欠定盲源分離系統(tǒng)模型Fig.1 Blind source separation system model

1.2 張量理論基礎(chǔ)

一維數(shù)組，我們稱之為向量；二維數(shù)組，我們稱之為矩陣；三維數(shù)組以及多維數(shù)組，我們稱之為張量[12]。張量由于其強(qiáng)大的數(shù)據(jù)分析能力，近年來在數(shù)據(jù)挖掘、機(jī)器學(xué)習(xí)等領(lǐng)域得到廣泛應(yīng)用。

下面介紹與張量有關(guān)的定義與定理。

定義1(向量的外積)[14]：向量a∈RI1×1與b∈RI2×1的外積是一個(gè)I1×I2維的矩陣C，定義為：

C=a·b=abT

(2)

定義2(秩1張量)[14]：如果三階張量χ等于三個(gè)向量的外積，則它的秩為1。

定義3(CP分解)[14]：任何張量χ都允許分解為秩-1張量的總和。在三階張量的情況下，CP分解定義為：

(3)

式(3)中，a∈RI1×1，b∈RI2×1，c∈RI3×1，F(xiàn)為張量χ的秩。

定義4(Kronecker積)[14]：I1×I2維矩陣A=[a1,…,aI2]與I3×I4維矩陣B的Kronecker積是一個(gè)維的矩陣，并有如下表達(dá)式：

(4)

定義5(Khatri-Ra積)[12]：I1×I0維矩陣A=[a1,a2,…,aI0] 和I2×I0維矩陣B=[b1,b2,…,bI0]的Khatri-Rao積A⊙B是一個(gè)I1I2×I0維的矩陣，并有如下表達(dá)式：

A⊙B=[a1?b1,…,aI0?bI0]

(5)

平行因子分析(PARAFAC)最重要的一個(gè)特征就是CP分解是具有唯一性的。它是采用平行因子模型進(jìn)行數(shù)據(jù)分析的首要條件。下面給出PARAFAC模型分解唯一性的一個(gè)重要定理。

定理1(PARAFAC模型的分解唯一性)[13]：若χ∈RI1×I2×I3是一個(gè)具有如式的CP分解張量，其中A∈RI1×F、B∈RI2×F和C∈RI3×F，如果滿足：

k(A)+k(B)+k(C)≥2F+2

(6)

則矩陣A，B和C在存在尺度模糊和列模糊的條件下是唯一的(也稱為本質(zhì)唯一)。

2 基于平行因子分析的欠定混合矩陣估計(jì)

本文利用通信信號(hào)的協(xié)方差矩陣構(gòu)造三階張量，將混合矩陣估計(jì)問題轉(zhuǎn)化為了張量分解問題，重點(diǎn)針對(duì)現(xiàn)有張量分解算法存在的局部收斂和初始矩陣選擇的問題提出改進(jìn)，采用非迭代的方法確定交替最小二乘算法的初始迭代矩陣，為了減少運(yùn)算時(shí)間，在迭代過程中采用標(biāo)準(zhǔn)線搜索加速收斂。

2.1 構(gòu)造張量

本文采用文獻(xiàn)[10]中的方法來構(gòu)造三階張量，具體方法如下:

對(duì)于零均值且互不相關(guān)的非平穩(wěn)實(shí)信號(hào)，源信號(hào)在t時(shí)刻的協(xié)方差矩陣可表示為：

(7)

式(7)中，Dk=E{stst+tk}是對(duì)角矩陣，k=1,…，K。為了簡(jiǎn)化，先忽略噪聲的影響。解決的問題是在M

可以將矩陣C1,…,Ck按如下方式壓縮為張量C∈M×M×K：(C)ijk=(Ck)ij，i=1,…,M，j=1,…,M，k=1,…,K。通過(D)kf=(Dk)ff定義一個(gè)矩陣Dkf，得到：

(8)

將其改寫為：

(9)

式(9)中，af和df分別代表矩陣A和D的列向量。

張量沿三個(gè)方向的切片展開為矩陣，C(1)∈M2×K，C(2)∈KM×M，C(3)∈MK×M：

(C(1))(i-1)M+j,k=(C(2))(k-1)K+i,j=(C(3))(j-1)M+k,i=Cijk

(10)

根據(jù)式(10)，張量的三個(gè)切片展開矩陣可以寫為

C(1)=(A⊙A)·DT

(11)

C(2)=(D⊙A)·AT

(12)

C(3)=(A⊙D)·AT

(13)

由定理1知，滿足一定條件，張量的CP分解唯一存在。文獻(xiàn)[10]給出了傳感器數(shù)目和源信號(hào)數(shù)目使CP分解唯一的關(guān)系，如表 1所示。

表1 觀測(cè)數(shù)M與所允許的源信號(hào)最大數(shù)目Nmax的關(guān)系Tab.1 The relationship between the number of observations M and the maximum number Nmaxof source signals

至此，將欠定混合矩陣的估計(jì)問題轉(zhuǎn)化為三階張量的CP分解問題。

2.2 ALS算法

作為解決PARAFAC模型的經(jīng)典算法，ALS算法算法步驟如圖 2所示。

圖2 ALS算法流程Fig.2 The ALS algorithm flowchart

算法通過交替最小二乘誤差γ，來獲得混合矩陣。

(14)

式(14)中，‖·‖F(xiàn)代表F-范數(shù)，當(dāng)矩陣A固定，D的估計(jì)值可以由最小二乘結(jié)果得到：

(15)

式(15)中，(·)+表示廣義逆矩陣。

同理，分別利用式(12)與式(13)，通過同樣的方式可得

(16)

(17)

由于初始矩陣的選擇對(duì)ALS算法的準(zhǔn)確性及收斂路徑有比較大的影響，本文采用非迭代的直接三線性分解法確定ALS算法初始矩陣。

2.3 直接三線性分解

直接三線性分解[16]算法是一種非迭代的求解平行因子模型的三維分解方法，但是其精度和可靠性比ALS算法要差，可以作為ALS算法的初始矩陣。

第一步 根據(jù)式(11)、式(12)和式(13)，分別對(duì)C(1)、C(2)和C(3)進(jìn)行奇異值分解，其中U、V取前F個(gè)奇異值矢量，W取前兩個(gè)左奇異值矢量。

C(1)=UISIVIU=UI(I×R)

(18)

C(2)=UJSJVJV=UJ(J×R)

(19)

C(3)=UKSKVKW=UK(K×2)

(20)

第二步 利用式(21)、式(22)構(gòu)造為樣本矩陣G1、G2。

(21)

(22)

式中，wki為矩陣W中的元素，Ck∈RM×M為張量C的第k個(gè)切片矩陣。

第三步 使用QZ分解陣G1、G2組成的特征方程，L、M分別為方程的特征向量，可證明矩陣A和D分別可由式(23)、式(24)得到[14]。

G1Lλ2r=G2Lλ1rA0=UL-1

(23)

G1Mλ2r=G2Mλ1rA0=VM-1

(24)

至此，通過直接三線性分解得到了粗估的加載矩陣，將其作為ALS算法的初始矩陣，開始迭代。

ALS算法運(yùn)算時(shí)間較長(zhǎng)，本文使用標(biāo)準(zhǔn)線搜索加速收斂。

2.4 算法流程

本文采用直接三線性分解粗估混合矩陣，將其作為ALS算法初始矩陣，迭代過程中采用標(biāo)準(zhǔn)的線搜索加速收斂。

線搜索方法首先找到使目標(biāo)函數(shù)減少的反向，然后計(jì)算沿著下降方向移動(dòng)的步長(zhǎng)。下降方向可以由各種不同的方法計(jì)算，如梯度下降法。

例如在t次迭代，根據(jù)算法預(yù)測(cè)一定的迭代步長(zhǎng)ρ，對(duì)矩陣A進(jìn)行更新，可以表示為：A=A+ρΔA。

本文采用的線加速方案為只對(duì)矩陣D進(jìn)行線搜索：

D(new)=D(it-2)+RLS(D(it-1)-D(it-2))

(25)

針對(duì)矩陣A，利用式(17)采用最小二乘進(jìn)一步估計(jì)。

算法步長(zhǎng)的計(jì)算應(yīng)該十分簡(jiǎn)單，如果它比相應(yīng)的迭代需要更多的時(shí)間則沒有意義。最簡(jiǎn)單的情況是，RLS被賦予固定值(在1.2～1.3之間)[17]，或者被設(shè)置為(it)1/3[18]。本文的RLS設(shè)置為(it)1/3。

具體算法流程：

步驟1 根據(jù)直接三線性分解算法粗估加載矩陣A0，D0；

步驟2 利用式(15)和A0估計(jì)出D1，利用式(16)和A0及D1估計(jì)出A1，令it=2，則A(it-2)=A0，D(it-2)=D0，A(it-1)=A1，D(it-1)=D1；

步驟3 根據(jù)線加速公式(25)計(jì)算D(new)，利用式(17)及A(it-1)D(it-1)計(jì)算A(new)；

步驟4 分別利用A(new)、D(new)和A(it-1)，D(it-1)及式(14)計(jì)算誤差函數(shù)γ(new)和γ(it-1)；

步驟5 若γ(new)<γ(it-1)，A(it-1)=A(new)，D(it-1)=D(new)，γ(it-1)=γ(new)；否則，直接執(zhí)行步驟6；

步驟6 利用式(15)和式(16)及A(it-1)，D(it-1)；計(jì)算得A(it)和D(it)：

(26)

(27)

3 仿真實(shí)驗(yàn)與分析

為說明算法的性能，算法的效果評(píng)估采用混合矩陣估計(jì)誤差和迭代次數(shù)作為指標(biāo)。

3.1 混合矩陣估計(jì)誤差

本節(jié)通過采用不同調(diào)制方式信號(hào)來完成欠定條件下混合矩陣估計(jì)實(shí)驗(yàn)。

實(shí)驗(yàn)采用4路不同調(diào)制方式的輻射源信號(hào)：s1為跳頻信號(hào)；s2為L(zhǎng)FM信號(hào)；s3為QPSK信號(hào)；s4為FM信號(hào)。觀測(cè)信號(hào)數(shù)目設(shè)為3，信號(hào)采樣點(diǎn)數(shù)為4 096，信噪比變化范圍為-5～30 dB，步進(jìn)為5 dB，每個(gè)信噪比處蒙特卡洛仿真實(shí)驗(yàn)100次。

評(píng)價(jià)混合矩陣估計(jì)性能采用歸一化均方誤差：

(28)

實(shí)驗(yàn)仿真結(jié)果如圖3所示，橫坐標(biāo)為信噪比，縱坐標(biāo)為歸一化均方誤差。參考算法1采用K均值聚類算法[19]，直接進(jìn)行迭代估計(jì)出混合矩陣。參考算法2采用ALS算法[10]，受初始值的選擇及局部收斂的影響，在低信噪比條件下混合矩陣估計(jì)性能一般。參考文獻(xiàn)[3]采用的AP聚類方法[20]，將每個(gè)樣本都視為潛在的聚類中心，通過迭代確定聚類中心。本文算法先用非迭代的方法確定初始矩陣，再進(jìn)行迭代確定混合矩陣。從圖 3可以看出，直接采用K均值聚類算法，接收陣元信噪比達(dá)到5 dB時(shí)， 就已達(dá)到極限，約為-15.5 dB，信噪比大于5 dB時(shí)，信噪比的增加無法提升混合矩陣的估計(jì)準(zhǔn)確性；采用本文算法、參考算法2以及參考算法3，混合矩陣估計(jì)的準(zhǔn)確性隨著信噪比的增加而增加，在低信噪比條件下，本文算法估計(jì)準(zhǔn)確性較參考算法2提升約3 dB，與參考算法3性能接近，在信噪比較高的條件下，本文算法的性能要明顯優(yōu)于參考算法3，略優(yōu)于參考算法2。對(duì)100次蒙特卡洛實(shí)驗(yàn)得到的歸一化均方誤差求方差，得到仿真結(jié)果，如圖4所示，橫坐標(biāo)為信噪比，縱坐標(biāo)為100次實(shí)驗(yàn)得到的均方誤差的方差。從圖中可以看出，本文算法均方誤差的方差位于ALS算法減小了約0.5，說明本文算法在估計(jì)的穩(wěn)定性上有了較大的提升。

圖3 均方誤差隨信噪比變化曲線Fig.3 Mean square error with SNR curve

3.2 迭代次數(shù)

本節(jié)通過設(shè)定不同的觀測(cè)數(shù)目與信源數(shù)來完成算法迭代次數(shù)估計(jì)實(shí)驗(yàn)。實(shí)驗(yàn)首先設(shè)置觀測(cè)數(shù)目及信源數(shù)，然后隨機(jī)生成A和D，進(jìn)行蒙特卡洛仿真實(shí)驗(yàn)100次，記錄歸一化均方誤差ENMSE達(dá)到10 dB時(shí)，所需的平均迭代次數(shù)。

圖4 均方誤差的方差隨信噪比變化曲線Fig.4 The variance of the mean square error varies with SNR curve

表2列出了兩種算法運(yùn)行一次的平均迭代次數(shù)，從表中可以看出在(M,N)取不同值條件下，本文算法迭代次數(shù)上較ALS算法減少約41.4%～84.3%，且混合矩形規(guī)模越小，減少幅度越大。

表2 (M,N)取不同值時(shí)算法的迭代次數(shù)Tab.2 Number of iterations of the algorithm when (M,N) take different values

4 結(jié)論

本文提出了基于平行因子分析的欠定混合矩陣估計(jì)算法。該算法利用信號(hào)的協(xié)方差矩陣構(gòu)造三階張量，采用直接三線性分解確定交替最小二乘算法的初始迭代矩陣，然后在迭代過程中采用標(biāo)準(zhǔn)線搜索加速收斂，最終實(shí)現(xiàn)張量分解得到混合矩陣。仿真實(shí)驗(yàn)表明，該方法不要求信源的稀疏性，較ALS算法估計(jì)精度可以提高約3 dB，迭代次數(shù)可以減少約41.4%～84.3%，是一種有效的欠定混合矩陣估計(jì)算法。