欠定條件下同步組網(wǎng)跳頻信號(hào)盲源分離方法

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

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

0 概述

跳頻擴(kuò)頻(Frequency-Hopping Spread Spectrum,FHSS)是指使用偽隨機(jī)碼進(jìn)行頻移鍵控,使得信號(hào)發(fā)射的載波頻率不斷跳變同時(shí)頻譜得到擴(kuò)展的方法[1]。FHSS技術(shù)具有保密性和組網(wǎng)能力強(qiáng)的優(yōu)點(diǎn),將其用于軍用通信,能夠有效避開干擾,保證通信效能。

欠定盲源分離是指在信源數(shù)多于觀測(cè)數(shù)的情況下,從接收到的混合信號(hào)中分離、恢復(fù)出源信號(hào)的過程[2]。跳頻信號(hào)的盲源分離是進(jìn)行信號(hào)偵察處理及電子對(duì)抗的一個(gè)關(guān)鍵步驟,是開展情報(bào)獲取、敵我識(shí)別、干擾引導(dǎo)和威脅等級(jí)分析等軍事行動(dòng)的前提。因此,研究跳頻信號(hào)盲源分離方法具有重要的理論意義和應(yīng)用價(jià)值。

通過判斷跳頻網(wǎng)臺(tái)之間跳變的時(shí)間是否同步,可以將跳頻網(wǎng)臺(tái)組網(wǎng)方式分為同步組網(wǎng)和異步組網(wǎng)兩種。異步組網(wǎng)是指在各跳頻網(wǎng)臺(tái)之間沒有時(shí)間同步的情況下,對(duì)跳頻信號(hào)用時(shí)頻分析工具得到每跳的到達(dá)時(shí)刻,然后根據(jù)不同電臺(tái)到達(dá)時(shí)刻和持續(xù)時(shí)長(zhǎng)進(jìn)行分離[3]。該方法簡(jiǎn)單有效,但其適用范圍僅局限于異步跳頻網(wǎng)臺(tái)。同步組網(wǎng)則是基于通信系統(tǒng)內(nèi)各跳頻電臺(tái)之間具有相同的時(shí)間同步,即電臺(tái)遵循相同跳變時(shí)間規(guī)律的情況?，F(xiàn)有欠定條件下同步組網(wǎng)跳頻信號(hào)的盲源分離方法主要采用“兩步法”,即第1步估計(jì)混合矩陣,第2步恢復(fù)源信號(hào)[4]。

混合矩陣的估計(jì)結(jié)果將直接影響后續(xù)源信號(hào)的恢復(fù)效果。文獻(xiàn)[5-7]方法采用基于稀疏聚類的混合矩陣估計(jì)方法,計(jì)算相對(duì)簡(jiǎn)單,估計(jì)精度較高,但是需要滿足一定的稀疏條件,應(yīng)用范圍有所限制。文獻(xiàn)[8-10]基于張量分解的估計(jì)算法,利用信號(hào)的高階統(tǒng)計(jì)量構(gòu)造張量,將混合矩陣的估計(jì)問題轉(zhuǎn)化為張量分解問題,適用于信源相互獨(dú)立并且非高斯分布的情況,精度較高,但是計(jì)算量較大。其中,文獻(xiàn)[8]方法將壓縮感知理論與張量分析相結(jié)合,在實(shí)現(xiàn)盲源分離的同時(shí)完成了波達(dá)方向(Direction of Arrival,DOA)估計(jì),但只適用于具有特殊表示形式的源信號(hào)。文獻(xiàn)[9-10]方法把張量分析應(yīng)用到盲源分離問題求解中,將混合矩陣的估計(jì)問題轉(zhuǎn)化為張量分解問題,提高了估計(jì)的準(zhǔn)確性和魯棒性,但其中用于CP(Canonical-Polyadic)分解的最小二乘(Alternating Least Square,ALS)算法存在對(duì)初值敏感、易陷入局部極值和運(yùn)算時(shí)間較長(zhǎng)的不足。

文獻(xiàn)[11-13]針對(duì)源信號(hào)恢復(fù)問題進(jìn)行了研究。文獻(xiàn)[11]采用稀疏重構(gòu)的方法恢復(fù)源信號(hào)。該方法可以實(shí)現(xiàn)稀疏信號(hào)的盲分離,但無法有效分離非稀疏的源信號(hào)。文獻(xiàn)[12]采用二元掩蔽法進(jìn)行源信號(hào)的恢復(fù)。該方法計(jì)算簡(jiǎn)單,適用于線性混合方式,是目前較為常用的欠定盲分離算法,源信號(hào)越稀疏其恢復(fù)效果越好,但在實(shí)際應(yīng)用中解決跳頻信號(hào)分離問題的性能并不出色。針對(duì)非充分稀疏的混合信號(hào),文獻(xiàn)[13]采用子空間投影法進(jìn)行源信號(hào)的重構(gòu),只要任意時(shí)頻點(diǎn)同時(shí)存在的源信號(hào)個(gè)數(shù)小于陣元數(shù),即可利用子空間投影原理確定任意時(shí)頻點(diǎn)上源信號(hào)對(duì)應(yīng)的混合矩陣。

本文針對(duì)欠定條件下跳頻信號(hào)的盲源分離問題,按照上文所述的“兩步法”,提出一種基于平行因子分析模型與子空間投影法的盲源分離方法。將混合矩陣估計(jì)問題轉(zhuǎn)化為張量CP分解問題,利用直接三線性分解改進(jìn)交替ALS算法對(duì)其求解,同時(shí)在迭代過程中采用標(biāo)準(zhǔn)線搜索加速收斂得到混合矩陣。在此基礎(chǔ)上,通過子空間投影法完成盲源分離,并剔除離散噪點(diǎn)進(jìn)一步優(yōu)化分離效果。

1 跳頻信號(hào)模型

跳頻信號(hào)sn(t)的數(shù)學(xué)表達(dá)式[14]描述為:

(1)

其中,Tn為跳頻信號(hào)sn(t)的跳周期,K為觀測(cè)時(shí)間Δt內(nèi)總跳數(shù),fnk為第k跳(k=1,2,…,K)的載頻,φnk為初始相位,起始的非完整跳在觀測(cè)時(shí)間內(nèi)的持續(xù)時(shí)長(zhǎng)為Δt0n,vn(t)是跳頻信號(hào)sn(t)的基帶復(fù)包絡(luò),t′=t-(k-1)Tn-Δt0n表示瞬時(shí)刻,rect表示單位矩形窗函數(shù)。

盲源分離是指從接收到的混合信號(hào)中分離、恢復(fù)出混合前源信號(hào)的過程。盲源分離系統(tǒng)模型如圖1所示。

圖1 盲源分離系統(tǒng)模型Fig.1 Model of blind source separation system

設(shè)M個(gè)觀測(cè)信號(hào)是來自于N個(gè)源信號(hào)的線性混合,則盲源分離信號(hào)模型描述為:

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

(2)

其中,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),N(t)表示與源信號(hào)獨(dú)立的加性噪聲,A為M×N維的混合矩陣。

在實(shí)際應(yīng)用中,根據(jù)信源數(shù)目和觀測(cè)信號(hào)數(shù)目的關(guān)系,可以將盲源分離模型分為超定、正定和欠定3類,超定即觀測(cè)數(shù)目大于源信號(hào)數(shù)目的情況,正定即觀測(cè)數(shù)目等于源信號(hào)數(shù)目的情況,欠定即觀測(cè)數(shù)目小于源信號(hào)數(shù)目的情況[15]。

2 混合矩陣估計(jì)

本文基于平行因子模型PARAFAC進(jìn)行混合矩陣估計(jì)。通過計(jì)算跳頻信號(hào)時(shí)延相關(guān)矩陣構(gòu)造三階張量,將混合矩陣估計(jì)問題轉(zhuǎn)化為張量CP分解問題。同時(shí)改進(jìn)ALS算法,利用直接三線性分解方法粗估加載矩陣,將其作為ALS算法的初始迭代矩陣,并在迭代過程中采用標(biāo)準(zhǔn)線搜索加速收斂得到混合矩陣。

2.1 張量分析相關(guān)概念

C=a°b=abT

(3)

定義2(秩-1張量)[16]如果三階張量χ等于3個(gè)向量的外積,則其秩為1。

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

(4)

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

(5)

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

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

(6)

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

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

(7)

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

2.2 張量構(gòu)造

根據(jù)文獻(xiàn)[9],對(duì)于互不相關(guān)的零均值非平穩(wěn)實(shí)信號(hào),源信號(hào)在t時(shí)刻的協(xié)方差矩陣可以表示為:

C1=E{xtxt+τ1}=A·D1·A

C2=E{xtxt+τ1}=A·D1·A

?

Ck=E{xtxt+τk}=A·Dk·A

(8)

其中,Dk=E{stst+tk}是對(duì)角矩陣,k=1,2,…,K。為簡(jiǎn)化運(yùn)算,首先忽略噪聲的影響,而需要解決的問題是在M

(9)

將式(9)改寫為:

(10)

其中,af和df分別代表矩陣A和D的列向量。

(C1)(i-1)M+j,k=(C2)(k-1)K+i,j=(C3)(j-1)M+k,i=Cijk

(11)

根據(jù)式(11),張量的3個(gè)切片展開矩陣可表示為:

C1=(A⊙A)·DT

(12)

C2=(D⊙A)·AT

(13)

C3=(A⊙D)·AT

(14)

由定理1可知,滿足一定條件時(shí)張量的CP分解唯一存在。文獻(xiàn)[9]給出了傳感器數(shù)目和源信號(hào)數(shù)目使CP分解唯一的關(guān)系,即觀測(cè)數(shù)M為2～8時(shí),所允許的源信號(hào)最大數(shù)目Nmax分別為2、4、6、10、15、20、26。至此,可將欠定混合矩陣的估計(jì)問題轉(zhuǎn)化為三階張量的CP分解問題。

2.3 基于CP分解的混合矩陣估計(jì)算法

作為解決PARAFAC模型的經(jīng)典算法,ALS算法通過交替最小二乘誤差γ來獲得加載矩陣,如式(15)所示:

(15)

其中,‖·‖F(xiàn)代表F-范數(shù)。矩陣A固定時(shí),D的估計(jì)值可由最小二乘結(jié)果得到,如式(16)所示:

(16)

其中,(·)+表示廣義逆矩陣。

對(duì)于式(13)與式(14),通過同樣的方式可得:

(17)

(18)

由于初始矩陣的選擇對(duì)ALS算法的準(zhǔn)確性及收斂路徑有較大影響,因此本文通過直接三線性分解粗估加載矩陣,將其作為ALS算法初始矩陣。

2.3.1 直接三線性分解

直接三線性分解[18]是一種非迭代的求解PARAFAC模型的三維分解方法,但其精度和可靠性比ALS算法要差,在基于PARAFAC模型的ALS算法中可以用于粗估初始矩陣。算法具體步驟如下:

步驟1根據(jù)式(13)～式(14),分別對(duì)C1、C2和C3進(jìn)行奇異值分解,其中,U、V取前F個(gè)奇異值矢量,W取前兩個(gè)左奇異值矢量,下標(biāo)I、J、K分別表示張量的行、列、通道,R為信源數(shù)目。

C1=UISIVI,U=UI(I×R)

(19)

C2=UJSJV,V=UJ(J×R)

(20)

C3=UKSKVK,W=UK(K×2)

(21)

步驟2利用式(22)、式(23)構(gòu)造樣本矩陣G1、G2:

(22)

(23)

步驟3使用QZ分解G1、G2組成的特征方程,令L、M分別為方程的特征向量,矩陣A和D分別可由式(24)、式(25)得到[13]:

G1Lλ2r=G2Lλ1r,A0=UL-1

(24)

G1Mλ2r=G2Mλ1r,A0=VM-1

(25)

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

2.3.2 ALS算法

ALS算法通過交替最小二乘誤差來獲得加載矩陣,算法具體步驟如下:

步驟1隨機(jī)確定初始矩陣Ait-1、Dit-1,利用直接三線性分解算法粗估加載矩陣,其中,it表示迭代次數(shù)。

ALS算法對(duì)初值敏感且容易陷入局部極值,同時(shí)運(yùn)算時(shí)間較長(zhǎng),因此,本文使用標(biāo)準(zhǔn)線搜索加速收斂。

2.3.3 算法流程

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

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

Dnew=Dit-2+RLS(Dit-1-Dit-2)

(26)

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

算法步長(zhǎng)的計(jì)算應(yīng)較為簡(jiǎn)單,如果它比相應(yīng)的迭代需要更多的時(shí)間則沒有意義。最簡(jiǎn)單的情況是RLS被賦予固定值(在1.2和1.3之間)或者被設(shè)置為it1/3[19]。本文中RLS設(shè)置為it1/3。改進(jìn)ALS算法的具體步驟如下:

步驟1利用直接三線性分解算法粗估加載矩陣A0、D0。

步驟2利用式(16)和A0估計(jì)D1,利用式(16)和A0、D1估計(jì)A1,令it=2,則Ait-2=A0,Dit-2=D0,Ait-1=A1,Dit-1=D1。

步驟3根據(jù)線加速公式(式(26))計(jì)算Dnew,利用式(18)和Ait-1、Dit-1計(jì)算Anew。

步驟4利用式(15)和Anew、Dnew、Ait-1、Dit-1計(jì)算誤差函數(shù)γnew和γit-1。

步驟5若γnew<γit-1,則Ait-1=Anew,Dit-1=Dnew,γit-1=γnew;否則,直接執(zhí)行步驟6。

步驟6利用式(16)、式(17)和Ait-1、Dit-1計(jì)算Ait和Dit:

(27)

(28)

3 源信號(hào)恢復(fù)

文獻(xiàn)[14]提出了基于子空間投影法的源信號(hào)恢復(fù)算法,在實(shí)際源信號(hào)恢復(fù)過程中設(shè)定一個(gè)與跳頻信號(hào)信源功率有關(guān)的門限,將其與混合矩陣的正交投影作比較來判斷任意時(shí)頻點(diǎn)同時(shí)存在的源信號(hào)個(gè)數(shù)是否小于陣元數(shù),并通過對(duì)離散噪點(diǎn)的進(jìn)一步剔除,達(dá)到更好的分離效果。跳頻源信號(hào)的相對(duì)功率定義為:

(29)

其中,‖·‖2代表二范數(shù),Lk表示源信號(hào)的時(shí)頻單源點(diǎn)集合所含有的向量個(gè)數(shù)。

假設(shè)時(shí)頻點(diǎn)(t′,f′)上對(duì)應(yīng)跳頻源信號(hào)的混合矩陣列向量為AR=[an1,an2,…,anR],R為每個(gè)時(shí)頻點(diǎn)上同時(shí)存在的信源個(gè)數(shù),則在該時(shí)頻點(diǎn)上的信號(hào)模型為:

X(t′,f′)=ARSR(t′,f′)+V(t′,f′)

(30)

其中,SR(t′,f′)=[Sn1(t′,f′),Sn2(t′,f′),…,SnR(t′,f′)]。

當(dāng)時(shí)頻點(diǎn)上的源信號(hào)個(gè)數(shù)小于觀測(cè)陣元個(gè)數(shù)時(shí),混合矩陣列向量AR的正交投影矩陣H為:

(31)

其中,I為單位矩陣。由此可以得到:

(32)

將式(32)代入式(30),在AR已知的情況下可以計(jì)算出源信號(hào)S(t′,f′),從而實(shí)現(xiàn)時(shí)頻點(diǎn)上源信號(hào)個(gè)數(shù)小于接收陣元個(gè)數(shù)時(shí)的盲源分離。但當(dāng)時(shí)頻點(diǎn)(t′,f′)上同時(shí)存在的源信號(hào)個(gè)數(shù)等于接收陣元的個(gè)數(shù)M′時(shí),則不能按此來計(jì)算,而是需要利用與源信號(hào)相對(duì)功率偏差相關(guān)的門限ε來判斷何種情況下源信號(hào)數(shù)目等于接收陣元的個(gè)數(shù)。本文將門限設(shè)定為相對(duì)功率的0.2倍。

當(dāng)R=M′時(shí),利用式(33)估計(jì)源信號(hào)S(t′,f′),實(shí)現(xiàn)盲源分離:

(33)

源信號(hào)恢復(fù)算法的具體步驟如下:

步驟1設(shè)R=M-1,基于式(30)計(jì)算式(31)。

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

設(shè)定接收陣元數(shù)目為M=3,陣元間距為2.5 m,接收機(jī)處理帶寬為[1 MHz,50 MHz],采樣頻率為100 MHz。跳頻信號(hào)的調(diào)制方式均為BPSK,頻率集與跳頻周期如表1所示。

表1 跳頻信號(hào)參數(shù)Table 1 Parameters of frequency-hopping signals

跳頻源信號(hào)的時(shí)域波形圖和STFT時(shí)頻圖如圖2和圖3所示,混合跳頻信號(hào)的時(shí)域波形圖如圖4所示,混合加噪跳頻信號(hào)的STFT時(shí)頻圖如圖5所示。

圖2 跳頻源信號(hào)的時(shí)域波形圖Fig.2 Time domain waveforms of frequency-hopping source signals

圖3 跳頻源信號(hào)的STFT時(shí)頻圖Fig.3 STFT time-frequency diagrams of frequency-hopping source signals

圖4 混合跳頻信號(hào)的時(shí)域波形圖Fig.4 Time domain waveform of mixed frequency-hopping signal

圖5 混合加噪跳頻信號(hào)的STFT時(shí)頻圖Fig.5 STFT time-frequency diagram of mixed and noisy frequency-hopping signal

4.1 混合矩陣估計(jì)

混合矩陣的估計(jì)結(jié)果將直接影響后續(xù)源信號(hào)的恢復(fù)效果。本文采用均方誤差評(píng)價(jià)混合矩陣的估計(jì)性能,如式(34)所示:

(34)

文獻(xiàn)[1]采用K均值聚類算法迭代估計(jì)出混合矩陣,文獻(xiàn)[2]采用ALS算法進(jìn)行估計(jì),文獻(xiàn)[3]采用的AP聚類方法將每個(gè)樣本都作為可能的聚類中心,最后通過迭代確定聚類中心。本文算法先用非迭代的方法粗估混合矩陣,再進(jìn)行迭代,確定混合矩陣。4種算法均方誤差隨信噪比的變化曲線如圖6所示?？梢钥闯?文獻(xiàn)[1]直接采用K均值聚類算法,當(dāng)接收陣元信噪比達(dá)到5 dB時(shí),E1極限值約為-15.5 dB,與其他算法相比估計(jì)準(zhǔn)確性較差;在低信噪比條件下,本文算法性能明顯優(yōu)于文獻(xiàn)[2]算法,與文獻(xiàn)[3]算法接近;在高信噪比條件下,本文算法性能明顯優(yōu)于文獻(xiàn)[3]算法,略優(yōu)于文獻(xiàn)[2]算法。

圖6 均方誤差隨信噪比的變化曲線1Fig.6 Variation curve 1 of mean square error with SNR

4.2 同步組網(wǎng)跳頻信號(hào)盲源分離

圖7和圖8分別給出了源信號(hào)數(shù)目為4、陣元數(shù)為3、接收陣元信噪比為3 dB時(shí),文獻(xiàn)[20]基于聚類與子空間投影法的盲源分離結(jié)果與本文算法盲源分離結(jié)果的STFT時(shí)頻圖?？梢钥闯?本文算法分選出的跳頻源信號(hào)時(shí)頻圖較文獻(xiàn)[20]算法噪點(diǎn)更少,跳頻信號(hào)的時(shí)頻圖中能量更集中,說明本文算法能更有效地實(shí)現(xiàn)盲源分離。

圖7 文獻(xiàn)[20]算法盲源分離結(jié)果Fig.7 Results of blind source separation by the algorithm proposed in reference[20]

圖8 本文算法盲源分離結(jié)果Fig.8 Results of blind source separation by the algorithm proposed in this paper

本文采用均方誤差衡量多跳頻信號(hào)網(wǎng)臺(tái)分選的性能,如式(35)所示:

(35)

圖9 均方誤差隨信噪比的變化曲線2Fig.9 Variation curve 2 of mean square error with SNR

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

在現(xiàn)代戰(zhàn)爭(zhēng)中,爭(zhēng)奪電磁頻譜使用和控制權(quán)的軍事斗爭(zhēng)日益激烈,通信組網(wǎng)趨向密集和頻繁。作為抗干擾通信的重要方式,欠定條件下同步組網(wǎng)跳頻信號(hào)的盲源分離具有重要意義。本文將欠定盲源分離中混合矩陣的估計(jì)問題轉(zhuǎn)化為觀測(cè)信號(hào)統(tǒng)計(jì)量所組成張量的CP分解問題,并通過改進(jìn)的ALS算法進(jìn)行求解,以滿足稀疏性要求。在源信號(hào)恢復(fù)時(shí)采用子空間投影法,并通過對(duì)離散噪點(diǎn)的進(jìn)一步剔除,達(dá)到更好的分離效果。仿真結(jié)果表明本文方法是一種有效的欠定盲源分離算法,但其在建立跳頻信號(hào)混合模型中采用了最簡(jiǎn)單的線性瞬時(shí)混合模型,適用范圍有限。后續(xù)將考慮時(shí)延、多徑干擾等更多影響因素,建立適用范圍更廣的混合模型。此外,本文方法雖然降低了運(yùn)算的迭代次數(shù),但是算法復(fù)雜度仍較高,運(yùn)算時(shí)間較長(zhǎng),下一步對(duì)此進(jìn)行改進(jìn)。