李紅光 郭 英 張東偉 楊銀松 齊子森 眭 萍

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

②(空軍工程大學空管領航學院 西安 710051)

③(空軍通信士官學校綜合訓練系 大連 116100)

1 引言

跳頻(Frequency Hopping, FH)同步組網(wǎng)由于組網(wǎng)數(shù)目多、網(wǎng)間干擾小和抗偵察能力強等優(yōu)點，在軍事通信中廣泛應用。該組網(wǎng)中的電臺FH周期、起跳時刻和頻率集均相同，傳統(tǒng)參數(shù)估計方法[1]很難進行有效分選，而且在實際偵察中，受天線體積限制，偵收的混合信號數(shù)目一般小于網(wǎng)臺個數(shù)，因此欠定盲源分離(Under-determined Blind Source Separation, UBSS)技術[2–5]是解決該類問題的有效方法。

目前基于稀疏分量分析[6–9]的“兩步法”是解決UBSS問題的常用方法。文獻[10]利用源信號時頻域的正交性，通過最大后驗模型估計二進制掩碼，完成源信號盲分離。該方法計算復雜度低，對混合信號稀疏性要求較高，只適用于2維空間。文獻[11]利用混合信號稀疏性，采用最大后驗概率的 l1范數(shù)分離算法實現(xiàn)多維源信號盲分離。該算法的分離效果和 l0范數(shù)算法相當，但對源信號稀疏性要求較高。文獻[12]改進了文獻[11]的目標函數(shù)，將稀疏盲源分離問題簡化為特定協(xié)方差矩陣的次特征向量的迭代計算問題，降低了計算復雜度。文獻[13]利用基追蹤貪婪算法實現(xiàn)UBSS源信號恢復，算法構造字典較復雜。上述文獻算法均是在定頻源信號條件下實現(xiàn)UBSS，對于FH信號盲分離效果不佳。文獻[14]利用Gabor變換將混合FH信號變換至時頻域并進行單源點檢測，通過近鄰傳播聚類方法對時頻比矩陣聚類完成混合矩陣估計。該方法提高了單源點檢測概率，但沒有考慮FH混合矩陣時變性。文獻[15]通過對S變換的時頻點進行單源點檢測建立時頻比矩陣，采用時頻比矩陣各列方差最小化準則估計混合矩陣。該方法門限閾值取值缺乏理論指導。文獻[16]利用歸一化的時頻點系數(shù)檢測單源點，通過層次聚類算法完成復數(shù)混合矩陣和波達方向(Direction Of Arrival, DOA)估計。該算法在低信噪比條件下，容易將部分噪聲點誤認為是單源點，導致混合矩陣估計誤差較大。針對上述問題，本文提出一種欠定盲源分離的同步FH網(wǎng)臺分選算法，該算法首先對FH信號時頻變換，結合L型陣列特點建立欠定盲源分離數(shù)學模型，并通過自適應閾值算法濾除時頻矩陣背景噪聲，采用絕對方位差檢測算法進行單源點檢測，保證單源點的充分稀疏性，然后利用基于戴維森堡丁指數(shù)(Davies-Bouldin Index, DBI)的模糊C值聚類算法完成混合矩陣和2維DOA估計，最后采用改進的稀疏自適應匹配追蹤算法實現(xiàn)源信號盲分離。

2 FH網(wǎng)臺盲分離數(shù)學模型

假設L型天線接收K 個遠場FH信號s(t)=[s1(t),s2(t),···,sK(t)]T，如 圖1所 示，子 陣 X 和 Y互 成90o，每個子陣有 M個陣元，M

將x ?y 平 面原點陣元作為參考，則第m 個陣元接收的混合信號可表示為其中，μmk表示第m 個陣元接收的第k 個FH信號幅度衰減，1 ≤m ≤M , 1 ≤k ≤K, τmk表 示第m 個陣元接收的第k 個FH信號傳輸時延，vm(t)表 示第m 個陣元接收的噪聲，fk(t)為 第t 時刻sk(t)的載頻。

假設接收陣列中各陣元是各向同性的，則混合信號在子陣 X和 Y的導向矢量分別為

其中，xX(t)=[xX,1(t),xX,2(t),···,xX,M(t)]T, xY(t)=[xY,1(t),xY,2(t),···,xY,M(t)]T, xX,m(t) 表示子陣X 的 第m 個 陣元接收的混合信號。eX(t)=[eX,1(t),eX,2(t),···,eX,M(t)]T, eY(t)=[eY,1(t),eY,2(t),···,eY,M(t)]T, eX,m(t) 表 示子陣 X 的第m 個陣元接收的均值為0，方差為σ2的高斯白噪聲。將式(4)和式(5)合并可得

圖1 FH信號L型陣列接收示意圖

為了保證源信號恢復精度，一般UBSS算法要求混合矩陣不變，考慮到FH載頻跳變特性，本文處理一幀數(shù)據(jù)為單跳周期數(shù)據(jù)，保證計算期間A(t)固定不變。通過短時傅里葉變換(Short-Time Fourier Transform, STFT)對混合FH信號時頻變換，采用文獻[17]跳時刻檢測算法完成同步FH信號的單跳分割，可得單跳FH信號在時頻域的UBSS數(shù)學模型，如式(7)所示

其中，X (t,f), S (t,f)和 E (t,f)分別表示接收的FH混合信號、源信號和噪聲的時頻變換。

3 FH混合矩陣估計

假設 sk(t)在第ζ 跳有 L個時頻單源點，其集合

由式(9)可知，在無噪聲條件下，觀測數(shù)據(jù)時頻比矩陣的列向量與混合矩陣AX(t)列向量只差一個復系數(shù)，則AX(t)列 向量估計值可表示為

由此可知，通過構造單源點時頻比矩陣即可完成子陣 X和Y 混合矩陣估計。

3.1 時頻矩陣去噪

為了避免噪聲對混合矩陣估計的影響，采用自適應閾值去噪算法，濾除時頻矩陣中噪聲和低能量時頻點。假設子陣 X接收的混合信號任一時頻點XX(th,fh), h ≥1 ， 通過式(11)來判斷XX(th,fh)是否被剔除

為了選取合適的閾值ε，提高去噪算法的魯棒性，ε 取值方法如下：

輸入：混合信號時頻矩陣XX,1。

輸出：最佳閾值εbest。

3.2 單源點選取

在同步非正交組網(wǎng)中會發(fā)生頻率碰撞情況，導致同一時頻點存在多個源FH信號，為了保證單源點的充分稀疏性，本文采用基于信號絕對方位差的檢測算法，對濾波后的時頻矩陣進行單源點檢測，從而得到單源點時頻矩陣X,1。 假設單源點P (tξ,fξ)存在一源信號Sχ(tξ,fξ)=0, 1 ≤χ ≤K，在不考慮噪聲條件下，子陣X 接收的觀測數(shù)據(jù)時頻矩陣XX(tξ,fξ)可表示為

對式(12)分別取實部和虛部，可得式(13)其中， εα=cos(?α) , ( tX,fX)∈X,1,X,1為去噪之后的時頻矩陣，0 <εα<1。

3.3 混合矩陣估計

模糊C均值 (Fuzzy C-Means, FCM) 聚類算法相比K均值算法更加精確合理，而且不易受初始聚類中心影響。傳統(tǒng)FCM目標函數(shù)如式(15)所示

由式(15)可知，F(xiàn)CM只是單純地將樣本點與聚類中心的歐氏距離作為類別分割標準，導致聚類對噪聲和離散點的影響較為敏感，聚類穩(wěn)定度和準確度降低。本文在式(15)基礎上引入懲罰項和模糊加權系數(shù)，降低噪聲和樣本集分布差異對聚類結果的影響。新目標函數(shù)如式(16)所示

其中，超參數(shù) ?j為懲罰項。由式(16)構建拉格朗日輔助目標函數(shù)

其中，λj表 示拉格朗日乘子，0 ≤λj≤1。對式(17)的ηi,j求 偏導，并令?(ηi,j,zi)/?ηi,j=0，可得=//fj?zi//2ηi,j

其中， ，由式(18)可求得

由式(19)和式(21)可得

式(22)即可更新FCM算法的 ηi,j，通過對?j簡化，將 N 個?j超 參數(shù)減少至對1個參數(shù)?′的設置。對式(17)zi求偏導，令? J?(ηi,j,zi)/?zi=0，可得

由式(23)可求得zi，如式(24)

將?′=?j/λj代 入式(21)，可得?j的更新式(25)

輸入：單源點矩陣X,m各時頻點所對應的頻率集 F，聚類個數(shù)設置C =M ，模糊系數(shù)β =2，最大迭代次數(shù)gmax， 目標函數(shù)收斂閾值為?，初始聚類中心 Z0為 M 個頻率點，初始化?j=0，迭代次數(shù)g =0。

輸出：最佳聚類個數(shù)Cbest， 聚類中心Zbest，分類集F′=[F1,F2,···,FCbest]。

步驟 1 根據(jù)式(22)計算更新ηg。

步驟 2 根據(jù)式(25)和式(24)計算更新Zg。

步驟 3 根據(jù)式(26)更新分類集合 F′。

步驟 5 根據(jù)聚類結果和式(27)計算DBI值I 。

步驟 6 令聚類個數(shù) C =M+1, g =0 , ?j=0，跳轉至步驟1，重新計算DBI值 I ，如果I 滿足IC>IC+1且IC+1

經過FCM-DBI聚類之后，根據(jù) F′即可得到所對應的時頻單源點集合′X,1，聚類中心Zbest即為各信源的跳頻率估計值集合X,1，而且每個跳頻率,1均與各單源點集合一一對應。最終根據(jù)式(11)可估計出混合矩陣列向量，經過對各陣元的重新組合排列，即可估計出子陣 X的混合矩陣， 同理可得子陣 Y 的混合矩陣。由式(2)和式(3)可得

對式(28)進行變換可得

由式(29)即可估計各FH信源方位角和俯仰角。

4 FH源信號恢復

輸入：第q 跳觀測信號的X (t,f)和。

輸出：源信號時頻矩陣(t,f)和時域數(shù)據(jù)。

步驟 1 初始化各參量， r0=X(t,f)，索引集Λ0=?， 列標號集Jlable=?， 支撐集Φ =?，迭代次數(shù)g =1， 步長Lsize=L0， 閾值Sstop和 Sth, Sstop?Sth。

在分離每跳源信號之后，可根據(jù)DOA信息對各跳源信號進行拼接，假設第1跳的第k 個源信號的DOA估計為,1和,1，第2跳的各源信號DOA估計為i,2和,2, i =1,2,···,Cbest，可根據(jù)式(32)完成信號拼接

由式(32)可知，第2跳的第 ic個FH源信號與第1跳的第k 個FH源信號來自同一FH電臺。

5 實驗仿真與分析

仿真實驗中，接收天線陣元數(shù) M =3，F(xiàn)H源信號數(shù) K =4 ，陣元間距d 為最大頻率信號波長的1/2，采樣率為20 MHz，觀測時長為6000采樣點，蒙特卡羅實驗次數(shù)均為1000，4個FH信源碼速率均為100 kbps，跳速均為1000 hop/s，調制方式均為BPSK，跳圖案和方位設置如表1所示。

5.1 單源點檢測

單源點數(shù)量和有效性決定混合矩陣估計的速度和精度，本文采用STFT方法，窗函數(shù)為1023點漢明窗，步進512點。圖2是信噪比0 dB時，STFT時頻變換去噪前后的時頻圖，由圖2(b)可知，本文自適應閾值去噪算法能夠去除時頻矩陣背景噪聲。

為了驗證算法有效性，對比分析了文獻[7,14,16]和本文4種算法的單源點檢測性能。圖3是4種檢測算法在信噪比0 dB時，對應混合矩陣列向量a3,n的時頻比實部和虛部散點圖，本次試驗εth=0.935。

由圖3(a)和圖3(b)可知，文獻[7]檢測的單源點數(shù)量較大，且沒有聚類特性，文獻[14]檢測的單源點具有一定的聚類特性，但具有很多噪聲點，主要是由于文獻[7]將去噪之后的時頻點均看作單源點，導致引入大量的多源點，文獻[14]在單源點檢測之前并沒有進行去噪處理，導致當信噪比較低時，會將一些能量較大的噪聲點誤判為單源點。由圖3(c)和圖3(d)可知，本文檢測的單源點相比文獻[16]的聚類特性更加明顯，而且單源點數(shù)量更少。主要由于文獻[16]是以時頻點能量大小作為單源點的判斷準則，信噪比變化和信道傳輸損耗對于該準則閾值的選擇具有較大影響，從而導致能量較小的多源點會誤判為單源點。本文是以時頻點的DOA信息作為判斷準則，只要信源方位具有一定的差異性，均能較準確地檢測出單源點，而且增加了自適應閾值的去噪預處理，使檢測算法具有較強的噪聲魯棒性。

5.2 混合矩陣估計

表1 各FH信源跳圖案和方位參數(shù)

圖2 STFT時頻變換去噪前后的時頻圖

為了對比分析混合矩陣的估計性能，本文采用均方誤差作為評價標準，如式(33)所示其中， A為實際混合矩陣，為估計混合矩陣，K 為源信號個數(shù)，// ·//F表示Frobenius范數(shù)。

圖4是單源點檢測閾值 εα在 信噪比?10 ～20 dB條件下對混合矩陣的估計性能影響，信噪比間隔3 dB, εα取值分別為0.90, 0.93, 0.96和0.99。實驗中FCM-DBI聚類個數(shù)初始值 C =3，最大迭代次數(shù)gmax=100 ，收斂閾值? =0.00001。由圖4可知，當檢測閾值 εα=0.99時，由于檢測的單源點數(shù)量過少，低信噪比條件下混合矩陣估計誤差較大，當閾值εα≤0.93時，由于引入大量多源點，同樣使混合矩陣估計不準。通過分析，檢測閾值 εα=0.945較為合理。

圖5是文獻[7,16,19]和本文4種混合矩陣估計算法在信噪比 ? 10 ～20 dB條 件下的AMSE隨信噪比變化情況，信噪比間隔3 dB。由圖5可知，文獻[7]和文獻[19]算法的混合矩陣估計精度相對較低，主要是由于文獻[7]單源點檢測時引入了較多的多源點，導致混合矩陣估計不準，文獻[19]利用濾波之后的觀測信號比值來估計混合矩陣，濾波器種類和參數(shù)選擇對估計結果影響較大，算法魯棒性較差。文獻[16]和本文算法的估計精度隨著信噪比增加均逐漸提高，當信噪比大于14 dB 時，兩種算法的估計精度接近，當信噪比小于5 dB時，本文算法的估計精度高于文獻[16]算法，主要是由于本文FCM - DBI算法通過引入DBI指數(shù)，不需要預先已知聚類個數(shù)，而且對FCM的目標函數(shù)增加了懲罰項和模糊加權系數(shù)，有效降低了由噪聲引起的聚類樣本分布不均衡對聚類結果的影響，增強了FCM - DBI聚類的魯棒性。而文獻[16]聚類算法不僅需要已知聚類個數(shù)，而且噪聲對聚類閾值的選擇有較大影響，特別是低信噪比條件下，算法對噪聲和異常樣本點較為敏感，導致聚類結果不穩(wěn)定，混合矩陣估計精度下降。

圖3 4種檢測算法的時頻比散點圖

εα AMSE圖4 不同信噪比下閾值 對 的影響

AMSE圖5 不同信噪比下4種算法的

5.3 源信號恢復

為了對比分析信源時域信號的恢復性能，本文采用信干比作為評價標準，如式(34)所示其中， Sk(t) 為第k 個信源時域信號，k(t)為恢復第k 個時域信號，SSIR越大說明源信號恢復效果越好。

圖6是文獻[16]子空間投影、SAMP和SASP 3種算法恢復的信號 SSIR隨信噪比變化情況，信噪比范圍為? 5 ～20 dB，間隔3 dB，3種恢復算法均采用本文所估計的混合矩陣，SASP的初始步長Lsize=4。由圖6可知，3種恢復算法的恢復性能均隨著信噪比增加而提高，當信噪比大于1 dB時，本文SASP算法的恢復性能優(yōu)于子空間投影和SAMP算法，主要是由于子空間投影算法需預先已知源信號個數(shù)，當假設的源信號數(shù)大于真實個數(shù)時，算法會引入額外噪聲，特別是低信噪比下，恢復誤差較大。而本文基于稀疏重構的SASP算法對噪聲具有一定的抑制作用，同時引入變步長的稀疏度估計方法，相比需已知稀疏度的SAMP算法具有更高的恢復精度。

根據(jù)文獻[7,16,19]和本文4種算法估計的混合矩陣，均采用SASP算法恢復信號。圖7是恢復信號SSIR隨信噪比變化情況，信噪比范圍? 5 ～20 dB，間隔3 dB。由圖7可知，在相同恢復算法條件下，本文估計的混合矩陣恢復性能優(yōu)于其他3種算法，而且混合矩陣估計精度越高，源信號恢復誤差越小，因此提高混合矩陣估計精度，是改善源信號恢復質量的有效方法之一。

5.4 DOA估計與FH信號拼接

為了對比分析DOA的估計性能，本文采用均方根誤差作為評價標準，如式(35)所示

SSIR圖6 不同信噪比下3種恢復算法的時域信號

圖8是文獻[16]和本文估計的DOA均方根誤差隨信噪比變化情況，信噪比范圍? 5 ～20 dB，間隔3 dB。由圖8可知，兩種算法估計的DOA均方根誤差均隨信噪比增加逐漸減少，當信噪比大于10 dB時，兩種算法的估計誤差均在1o左右，但在低信噪比條件下，本文算法估計的DOA誤差更小，主要是由于本文的混合矩陣在低信噪比下估計精度更高。

表2是不同信噪比下估計的相鄰兩跳DOA(方位角/俯仰角)和均方根誤差。通過表1和表2可知，當信噪比大于10 dB時，DOA估計誤差小于1°。同時相鄰兩跳的同一信源DOA估計非常接近，通過式(31)可以完成FH信號拼接。

6 結束語

本文通過自適應閾值去噪和方位差單源點檢測方法，有效提高了混合矩陣估計精度，并利用SASP算法完成各跳源信號重構恢復，根據(jù)同一信源的DOA估計一致原則實現(xiàn)各跳源信號拼接。當SNR大于5 dB時，所提算法恢復的信號干信比達到20 dB, DOA估計精度小于1°，如何進一步降低SASP算法的計算復雜度，是下一步研究方向。

SSIR圖7 不同信噪比下4種混合矩陣恢復的時域信號

圖8 不同信噪比下兩種算法估計的DOA均方根誤差

表2 不同信噪比下相鄰兩跳的DOA和RMSE(°)