欠定條件下同步組網跳頻信號盲源分離方法

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

(空軍工程大學 信息與導航學院,西安 710077)

0 概述

跳頻擴頻(Frequency-Hopping Spread Spectrum,FHSS)是指使用偽隨機碼進行頻移鍵控,使得信號發射的載波頻率不斷跳變同時頻譜得到擴展的方法[1]。FHSS技術具有保密性和組網能力強的優點,將其用于軍用通信,能夠有效避開干擾,保證通信效能。

欠定盲源分離是指在信源數多于觀測數的情況下,從接收到的混合信號中分離、恢復出源信號的過程[2]。跳頻信號的盲源分離是進行信號偵察處理及電子對抗的一個關鍵步驟,是開展情報獲取、敵我識別、干擾引導和威脅等級分析等軍事行動的前提。因此,研究跳頻信號盲源分離方法具有重要的理論意義和應用價值。

通過判斷跳頻網臺之間跳變的時間是否同步,可以將跳頻網臺組網方式分為同步組網和異步組網兩種。異步組網是指在各跳頻網臺之間沒有時間同步的情況下,對跳頻信號用時頻分析工具得到每跳的到達時刻,然后根據不同電臺到達時刻和持續時長進行分離[3]。該方法簡單有效,但其適用范圍僅局限于異步跳頻網臺。同步組網則是基于通信系統內各跳頻電臺之間具有相同的時間同步,即電臺遵循相同跳變時間規律的情況。現有欠定條件下同步組網跳頻信號的盲源分離方法主要采用“兩步法”,即第1步估計混合矩陣,第2步恢復源信號[4]。

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

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

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

1 跳頻信號模型

跳頻信號sn(t)的數學表達式[14]描述為:

(1)

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

盲源分離是指從接收到的混合信號中分離、恢復出混合前源信號的過程。盲源分離系統模型如圖1所示。

圖1 盲源分離系統模型Fig.1 Model of blind source separation system

設M個觀測信號是來自于N個源信號的線性混合,則盲源分離信號模型描述為:

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

(2)

其中,X(t)=[x1(t),x2(t),…,xM(t)]T表示接收的M個觀測信號,S(t)=[s1(t),s2(t),…,sN(t)]T表示輸入的N個未知源信號,N(t)表示與源信號獨立的加性噪聲,A為M×N維的混合矩陣。

在實際應用中,根據信源數目和觀測信號數目的關系,可以將盲源分離模型分為超定、正定和欠定3類,超定即觀測數目大于源信號數目的情況,正定即觀測數目等于源信號數目的情況,欠定即觀測數目小于源信號數目的情況[15]。

2 混合矩陣估計

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

2.1 張量分析相關概念

C=a°b=abT

(3)

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

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

(4)

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

(5)

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

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

(6)

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

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

(7)

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

2.2 張量構造

根據文獻[9],對于互不相關的零均值非平穩實信號,源信號在t時刻的協方差矩陣可以表示為:

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}是對角矩陣,k=1,2,…,K。為簡化運算,首先忽略噪聲的影響,而需要解決的問題是在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)

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

C1=(A⊙A)·DT

(12)

C2=(D⊙A)·AT

(13)

C3=(A⊙D)·AT

(14)

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

2.3 基于CP分解的混合矩陣估計算法

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

(15)

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

(16)

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

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

(17)

(18)

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

2.3.1 直接三線性分解

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

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

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

(19)

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

(20)

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

(21)

步驟2利用式(22)、式(23)構造樣本矩陣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隨機確定初始矩陣Ait-1、Dit-1,利用直接三線性分解算法粗估加載矩陣,其中,it表示迭代次數。

ALS算法對初值敏感且容易陷入局部極值,同時運算時間較長,因此,本文使用標準線搜索加速收斂。

2.3.3 算法流程

在t次迭代,利用算法預測一定的迭代步長ρ,對矩陣A進行更新,可以表示為:A=A+ρΔA。

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

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

(26)

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

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

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

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

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

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

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

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

(27)

(28)

3 源信號恢復

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

(29)

其中,‖·‖2代表二范數,Lk表示源信號的時頻單源點集合所含有的向量個數。

假設時頻點(t′,f′)上對應跳頻源信號的混合矩陣列向量為AR=[an1,an2,…,anR],R為每個時頻點上同時存在的信源個數,則在該時頻點上的信號模型為:

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

(30)

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

當時頻點上的源信號個數小于觀測陣元個數時,混合矩陣列向量AR的正交投影矩陣H為:

(31)

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

(32)

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

當R=M′時,利用式(33)估計源信號S(t′,f′),實現盲源分離:

(33)

源信號恢復算法的具體步驟如下:

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

4 仿真與結果分析

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

表1 跳頻信號參數Table 1 Parameters of frequency-hopping signals

跳頻源信號的時域波形圖和STFT時頻圖如圖2和圖3所示,混合跳頻信號的時域波形圖如圖4所示,混合加噪跳頻信號的STFT時頻圖如圖5所示。

圖2 跳頻源信號的時域波形圖Fig.2 Time domain waveforms of frequency-hopping source signals

圖3 跳頻源信號的STFT時頻圖Fig.3 STFT time-frequency diagrams of frequency-hopping source signals

圖4 混合跳頻信號的時域波形圖Fig.4 Time domain waveform of mixed frequency-hopping signal

圖5 混合加噪跳頻信號的STFT時頻圖Fig.5 STFT time-frequency diagram of mixed and noisy frequency-hopping signal

4.1 混合矩陣估計

混合矩陣的估計結果將直接影響后續源信號的恢復效果。本文采用均方誤差評價混合矩陣的估計性能,如式(34)所示:

(34)

文獻[1]采用K均值聚類算法迭代估計出混合矩陣,文獻[2]采用ALS算法進行估計,文獻[3]采用的AP聚類方法將每個樣本都作為可能的聚類中心,最后通過迭代確定聚類中心。本文算法先用非迭代的方法粗估混合矩陣,再進行迭代,確定混合矩陣。4種算法均方誤差隨信噪比的變化曲線如圖6所示。可以看出:文獻[1]直接采用K均值聚類算法,當接收陣元信噪比達到5 dB時,E1極限值約為-15.5 dB,與其他算法相比估計準確性較差;在低信噪比條件下,本文算法性能明顯優于文獻[2]算法,與文獻[3]算法接近;在高信噪比條件下,本文算法性能明顯優于文獻[3]算法,略優于文獻[2]算法。

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

4.2 同步組網跳頻信號盲源分離

圖7和圖8分別給出了源信號數目為4、陣元數為3、接收陣元信噪比為3 dB時,文獻[20]基于聚類與子空間投影法的盲源分離結果與本文算法盲源分離結果的STFT時頻圖。可以看出,本文算法分選出的跳頻源信號時頻圖較文獻[20]算法噪點更少,跳頻信號的時頻圖中能量更集中,說明本文算法能更有效地實現盲源分離。

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

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

本文采用均方誤差衡量多跳頻信號網臺分選的性能,如式(35)所示:

(35)

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

5 結束語

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