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

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

0 引言

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

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

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

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

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

1.1 欠定盲源分離模型

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

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

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

(1)

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

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

1.2 張量理論基礎(chǔ)

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

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

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

C=a·b=abT

(2)

定義2(秩1張量)[14]：如果三階張量χ等于三個向量的外積，則它的秩為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積是一個維的矩陣，并有如下表達式：

(4)

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

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

(5)

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

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

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

(6)

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

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

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

2.1 構(gòu)造張量

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

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

(7)

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

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

(8)

將其改寫為：

(9)

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

張量沿三個方向的切片展開為矩陣，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)，張量的三個切片展開矩陣可以寫為

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

(11)

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

(12)

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

(13)

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

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

至此，將欠定混合矩陣的估計問題轉(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的估計值可以由最小二乘結(jié)果得到：

(15)

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

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

(16)

(17)

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

2.3 直接三線性分解

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

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

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個切片矩陣。

第三步 使用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算法運算時間較長，本文使用標準線搜索加速收斂。

2.4 算法流程

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

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

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

本文采用的線加速方案為只對矩陣D進行線搜索：

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

(25)

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

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

具體算法流程：

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

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

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

步驟4 分別利用A(new)、D(new)和A(it-1)，D(it-1)及式(14)計算誤差函數(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)；計算得A(it)和D(it)：

(26)

(27)

3 仿真實驗與分析

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

3.1 混合矩陣估計誤差

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

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

評價混合矩陣估計性能采用歸一化均方誤差：

(28)

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

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

3.2 迭代次數(shù)

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

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

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

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

4 結(jié)論

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