一種基于K-means聚類的跳頻信號(hào)快速檢測(cè)方法*

(火箭軍工程大學(xué) 導(dǎo)彈工程學(xué)院，西安 710025)

0 引 言

無人機(jī)產(chǎn)業(yè)的興起為社會(huì)提供了諸多便利，但是“黑飛”無人機(jī)亦帶來了一系列的安全隱患，進(jìn)一步引起了國內(nèi)外對(duì)反無人機(jī)技術(shù)的極大關(guān)注。跳頻通信因具有較強(qiáng)的抗干擾、抗截獲和易組網(wǎng)等能力[1]被廣泛應(yīng)用于無人機(jī)的信息傳輸，而跳頻信號(hào)的檢測(cè)是實(shí)現(xiàn)對(duì)無人機(jī)信號(hào)估計(jì)與干擾反制的重要前提。

目前，針對(duì)跳頻信號(hào)的提取檢測(cè)方法眾多，可分為非盲檢測(cè)算法和盲檢測(cè)算法兩類。非盲檢測(cè)算法是指預(yù)先知道跳頻信號(hào)部分參數(shù)，通過已知參數(shù)反映出跳頻信號(hào)的某些特征，從而實(shí)現(xiàn)對(duì)跳頻信號(hào)的檢測(cè)的一類算法。文獻(xiàn)[2]采用單跳自相關(guān)技術(shù)進(jìn)行預(yù)檢測(cè)處理，在跳速預(yù)知時(shí)，算法具有優(yōu)越的抗噪性能。文獻(xiàn)[3]提出基于多跳自相關(guān)的跳頻信號(hào)檢測(cè)方法，在已知跳頻信號(hào)周期的前提下計(jì)算出跳頻信號(hào)的自相關(guān)函數(shù)，通過判別自相關(guān)函數(shù)在不同時(shí)延情況下的值來實(shí)現(xiàn)跳頻信號(hào)的檢測(cè)。該算法計(jì)算量小，易于工程實(shí)現(xiàn)。文獻(xiàn)[4-6]提出基于輻射計(jì)的信道化檢測(cè)方法，跳頻信道數(shù)已知時(shí)利用檢測(cè)模型對(duì)各信道并行分析處理，根據(jù)各信道的檢測(cè)結(jié)果判斷跳頻信號(hào)的有無，檢測(cè)性能良好，但計(jì)算復(fù)雜。文獻(xiàn)[7-8]提出的對(duì)跳頻信號(hào)進(jìn)行分層處理的檢測(cè)模型，使得利用信道化檢測(cè)時(shí)接收信道中的無效信道數(shù)量減少，降低了模型計(jì)算復(fù)雜度。文獻(xiàn)[9]將接收信號(hào)進(jìn)行主成分分析，降低數(shù)據(jù)維數(shù)，再進(jìn)行數(shù)字信道化處理，通過短時(shí)能量對(duì)消實(shí)現(xiàn)跳頻信號(hào)檢測(cè)，降低了計(jì)算復(fù)雜度。以上方法均屬于非盲檢測(cè)，但實(shí)際戰(zhàn)場環(huán)境中難以獲取先驗(yàn)信息，同時(shí)也無法實(shí)現(xiàn)跳頻信號(hào)的非盲檢測(cè)。

盲檢測(cè)類算法是指在不需跳頻信號(hào)任何參數(shù)為先驗(yàn)信息的情況下對(duì)其進(jìn)行檢測(cè)的一類算法。在盲檢測(cè)算法中，時(shí)頻分析可以清晰地反映跳頻信號(hào)的時(shí)頻特性，因而被廣泛應(yīng)用。文獻(xiàn)[1]利用短時(shí)傅里葉變換和 Wigner-Ville分布組合時(shí)頻變換獲取時(shí)頻矩陣，去除背景噪聲后進(jìn)行二值化處理，利用連通域標(biāo)記算法對(duì)同一信號(hào)像素點(diǎn)進(jìn)行標(biāo)記并提取信息，通過聚類算法實(shí)現(xiàn)信號(hào)的識(shí)別檢測(cè)。文獻(xiàn)[10]提出了一種基于局部自適應(yīng)閾值的跳頻信號(hào)檢測(cè)方法，可以有效地對(duì)時(shí)頻矩陣進(jìn)行噪聲處理，但抗噪性能較差，無法實(shí)現(xiàn)低信噪比條件下的跳頻信號(hào)檢測(cè)。文獻(xiàn)[11-12]將時(shí)頻圖看作二維圖像，根據(jù)圖像邊緣提取檢測(cè)算子，有效去除噪聲，但計(jì)算復(fù)雜度高，難以實(shí)時(shí)處理。文獻(xiàn)[13]對(duì)跳頻信號(hào)進(jìn)行稀疏重構(gòu)，利用形態(tài)學(xué)濾波去除噪聲，算法抗噪性好但復(fù)雜度高，難以工程實(shí)現(xiàn)。文獻(xiàn)[14]提出了一種基于滑動(dòng)相關(guān)與小波變換結(jié)合的跳頻信號(hào)檢測(cè)方法，可較好地克服信道衰落和傳播損耗的影響，適用于復(fù)雜電磁環(huán)境，但是檢測(cè)效果受限于小波基的選取。以上算法均面臨著計(jì)算復(fù)雜度高或者抗噪性差的問題，難以應(yīng)用于復(fù)雜的戰(zhàn)場環(huán)境。

針對(duì)低信噪比環(huán)境下跳頻信號(hào)快速檢測(cè)的需求，本文根據(jù)時(shí)頻矩陣中隨機(jī)噪聲幅值在一定信噪比條件下通常弱于跳頻信號(hào)且噪聲點(diǎn)區(qū)域數(shù)據(jù)總是大于信號(hào)區(qū)域數(shù)據(jù)點(diǎn)數(shù)[1,12,15]的特點(diǎn)，提出了一種基于時(shí)頻矩陣每一頻率分量最大值進(jìn)行K-means聚類的閾值方法，極大地減少了閾值選取以及數(shù)據(jù)截?cái)嗟挠?jì)算復(fù)雜度，并且實(shí)現(xiàn)了低信噪比條件下的跳頻信號(hào)檢測(cè)。

1 信號(hào)預(yù)處理

1.1 信號(hào)模型

跳頻信號(hào)的數(shù)學(xué)表達(dá)式為

(1)

式中：k≥1，k為整數(shù)；Th為跳頻周期；A為信號(hào)幅度;θ為初始相位;fk為第k個(gè)時(shí)隙的跳頻頻率;rectTh表示矩形窗，且滿足

(2)

由于跳頻頻帶較寬，故接收機(jī)接收的信號(hào)中往往摻雜著高斯白噪聲，因此實(shí)際環(huán)境中跳頻信號(hào)的接收模型表達(dá)式為

x(t)=s(t)+n(t)=

(3)

式中：n(t)是加性高斯白噪聲。

1.2 時(shí)頻變換

目前時(shí)頻變換方法主要有短時(shí)傅里葉變換[16](Short Time Fourier Transform，STFT)、Wigner-Ville分布[17](Wigner-Ville Distribution,WVD)、偽Wigner-Ville分布[18](Peudo Wigner-Ville Distribution,PWVD)和平滑偽Wigner-Ville分布[19](Smooth PWVD,SPWVD) 。STFT 計(jì)算復(fù)雜度低，且無交叉項(xiàng)影響；WVD 交叉項(xiàng)嚴(yán)重；PWVD 和 SPWVD 通過增加窗函數(shù)，對(duì)交叉項(xiàng)有一定的抑制作用，但計(jì)算復(fù)雜，難以工程實(shí)現(xiàn)?；诠こ虒?shí)際的需求，本文選取STFT作為時(shí)頻分析工具。

(4)

式中：ω(t,τ)表示窗函數(shù)，τ表示時(shí)延量。STFT(t,f)表示信號(hào)x(t)時(shí)間和頻率的二維分布。

信號(hào)經(jīng)過短時(shí)傅里葉變換，假設(shè)產(chǎn)生一個(gè)M×N維矩陣以及對(duì)應(yīng)的時(shí)頻圖，如圖1所示。

(a)時(shí)頻矩陣部分示意圖

在時(shí)頻矩陣和二維時(shí)頻圖中，橫向表示時(shí)間，縱向表示歸一化頻率。時(shí)頻矩陣的每一行元素表示信號(hào)在時(shí)間窗內(nèi)經(jīng)過傅里葉變換產(chǎn)生的對(duì)應(yīng)頻率分量的幅值。圖1(a)中紅色、黃色分別表示信號(hào)以及信號(hào)的能量泄露，綠色部分表示噪聲。由圖1(a)可以看出，在一定信噪比條件下，時(shí)頻矩陣中噪聲的數(shù)量遠(yuǎn)遠(yuǎn)大于信號(hào)，分布廣泛且幅值較小?；诖?，本文通過K-means算法實(shí)現(xiàn)對(duì)噪聲與信號(hào)的聚類區(qū)分，并設(shè)置閾值濾除噪聲及信號(hào)泄露分量，進(jìn)而基于駐留時(shí)間實(shí)現(xiàn)對(duì)跳頻信號(hào)的檢測(cè)。

2 基于K-means聚類的跳頻信號(hào)檢測(cè)

2.1 K-means聚類原理

聚類分割是把給定的樣本集合按照某種準(zhǔn)則分割成若干個(gè)不相交的子集，每個(gè)子集中的樣本之間相似性較大，不同子集樣本之間的相似性較小。K-means聚類由于其實(shí)現(xiàn)簡單、聚類效果較好而被廣泛應(yīng)用。

K-means聚類算法采用距離作為相似性的評(píng)價(jià)指標(biāo)，即認(rèn)為兩個(gè)對(duì)象的距離越近相似度越大，距離越遠(yuǎn)相似度越小。對(duì)于給定的樣本集，按照樣本之間的距離大小，將樣本集劃分為k個(gè)集合。讓集合內(nèi)的點(diǎn)盡量緊密地連在一起，而讓集合間的距離盡量地大，集合內(nèi)的相似性越大，集合間的差別越大，則聚類效果越好。

假設(shè)將樣本集劃分為k個(gè)集合(C1,C2,C3,…Ck)，則K-means的誤差平方和準(zhǔn)則函數(shù)表達(dá)式為

(5)

式中：ui是集合Ci的均值向量，有時(shí)也稱為質(zhì)心，其表達(dá)式為

(6)

式中：Ni是第i集合Ci的樣本數(shù)目。在K-means聚類過程中，E度量了所有樣本子集的總誤差平方和。對(duì)于樣本集的不同分類，導(dǎo)致不同Ci中的不同的樣本子集以及對(duì)應(yīng)的均值ui，從而會(huì)得到不同的E值，而當(dāng)E最小時(shí)聚類效果最佳。

2.2 基于K-means聚類的閾值選取與數(shù)據(jù)截?cái)?/h3>

假設(shè)時(shí)頻矩陣閾值為T，將矩陣中大于閾值的元素保留，小于閾值的置0，則有

(7)

式中：i和j分別表示時(shí)頻矩陣的第i行和第j列，且1≤i≤M,1≤j≤N。

由時(shí)頻矩陣可以看出，跳頻信號(hào)具有時(shí)頻聚集性，并且幅值遠(yuǎn)遠(yuǎn)高于噪聲。因此，可以認(rèn)為當(dāng)時(shí)頻矩陣中縱軸的歸一化頻率有信號(hào)時(shí)，其最大值為信號(hào)幅值；而沒有信號(hào)時(shí)，因?yàn)榫鶠樵肼?，則所有元素均遠(yuǎn)小于信號(hào)幅值?；诖耍x取時(shí)頻矩陣中每一行的最大值，組成一個(gè)新的向量數(shù)組作為K-means聚類的向量集，極大地減少了聚類所產(chǎn)生的迭代次數(shù)。

Kmeans={STFT(1,max),STFT(2,max),…,STFT(M,max)}。

(8)

式中：Kmeans表示用于K-means聚類的含有M個(gè)元素的向量集，STFT(i,max)表示時(shí)頻矩陣中第i行的最大值。

對(duì)向量集樣本進(jìn)行聚類，步驟如下：

Step1 把Kmeans中的M個(gè)樣本數(shù)據(jù)分成k個(gè)聚類的初始劃分，計(jì)算每個(gè)聚類的均值u1,u2,…uk和E。在本文中需要實(shí)現(xiàn)信號(hào)的檢測(cè)，故取k=2，即將樣本STFT(i,max)聚為噪聲和信號(hào)兩類。

Step2 從數(shù)據(jù)集中隨機(jī)選擇k個(gè)數(shù)據(jù)點(diǎn)作為質(zhì)心。

Step3 對(duì)數(shù)據(jù)集中每一個(gè)點(diǎn)，計(jì)算其與每一個(gè)質(zhì)心的距離，并劃分到距離最近的質(zhì)心所屬集合。

Step4 把所有數(shù)據(jù)歸為k個(gè)集合，重新計(jì)算每個(gè)集合的質(zhì)心以及ui的值，并修改E。

Step5 若連續(xù)迭代多次E不變則停止，否則轉(zhuǎn)到Step 3。

樣本STFT(i,max)聚類前后如圖2所示。

(a)樣本數(shù)據(jù)聚類前

由圖2(a)可以看出，樣本數(shù)據(jù)在聚類前就有明顯的分類跡象，下方密集數(shù)據(jù)表示的是噪聲，上方零星分布的是跳頻信號(hào)，通過聚類進(jìn)行準(zhǔn)確判別，可以認(rèn)為圖2(b)中的藍(lán)色部分為噪聲部分，紅色為跳頻信號(hào)。將聚類后的兩類數(shù)據(jù)分別存放于矩陣Noise和矩陣Sig中，因?yàn)閷oise矩陣中的所有元素默認(rèn)為噪聲，并且所有元素均為每一行的最大值，故取Noise矩陣中最大值Noise(max)作為閾值，時(shí)頻矩陣中所有小于等于Noise(max)的值均判別為噪聲，并置0，則式(7)可重新表示為

(9)

由于STFT(i,max)表示每一行中的最大值，當(dāng)最大值小于閾值時(shí)，其所在行的每一個(gè)元素均小于閾值，表示該行沒有跳頻信號(hào)，全部置0；當(dāng)最大值大于閾值時(shí)，對(duì)該行的每一個(gè)元素同閾值進(jìn)行比較，小于閾值的依然判定為噪聲并置0，大于閾值的元素保留，則式(9)可以進(jìn)一步改寫為

(10)

通過時(shí)頻矩陣以及二維時(shí)頻圖可以發(fā)現(xiàn)，跳頻信號(hào)對(duì)應(yīng)的頻率分量在縱軸頻率域占比很小。因此，基于式(10)可以對(duì)時(shí)頻矩陣中所有元素快速地判別截?cái)?，再次大幅度減少了計(jì)算量，提高了檢測(cè)效率。

截?cái)嗵幚頃?huì)使跳頻信號(hào)的部分邊緣丟失，但是由圖1(b)不難發(fā)現(xiàn)，在一定的信噪比條件下，跳頻信號(hào)依然保留絕大部分的有用信息，即存在一定的駐留時(shí)間，在時(shí)頻矩陣中就表現(xiàn)為時(shí)間維度上一段幅值不為0的連續(xù)序列，因此可以根據(jù)時(shí)域上的連續(xù)性來區(qū)別信號(hào)和噪聲。當(dāng)非0序列有一定長度并且在頻域上依次變化時(shí)，則認(rèn)為接收的信號(hào)中可能包含有跳頻信號(hào)，從而實(shí)現(xiàn)對(duì)跳頻信號(hào)的檢測(cè)。

2.3 算法流程

通過上述的理論分析，初步認(rèn)為本文所提算法在低信噪比條件下實(shí)現(xiàn)跳頻信號(hào)的快速檢測(cè)是可行的。算法流程如下：

Step1 采集一段足夠長時(shí)間的數(shù)據(jù)，計(jì)算短時(shí)傅里葉變換，得到時(shí)頻矩陣。

Step2 取時(shí)頻矩陣每一行的最大值組成樣本集Kmeans。

Step3 通過K-means聚類算法將樣本集中的元素劃分為噪聲和信號(hào)兩類，并將各自樣本數(shù)據(jù)分別存放在矩陣Noise和Sig中。

Step4 選取矩陣Noise中元素最大值作為閾值。

Step5 根據(jù)式(10)，對(duì)時(shí)頻矩陣進(jìn)行截?cái)嗵幚恚瑸V除噪聲，得到新的時(shí)頻矩陣。

Step6 根據(jù)跳頻信號(hào)在駐留時(shí)間內(nèi)的連續(xù)性，對(duì)新的時(shí)頻矩陣進(jìn)行檢測(cè)。若非0序列均有一定的長度且在頻率上依次出現(xiàn)，則認(rèn)為存在跳頻信號(hào)。

跳頻信號(hào)檢測(cè)算法流程如圖3所示。

圖3 跳頻信號(hào)檢測(cè)算法流程圖

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

3.1 本文算法在低信噪比條件下的檢測(cè)概率

根據(jù)式(2)產(chǎn)生一段跳頻信號(hào)，包含8個(gè)跳頻周期，跳頻頻率集為{2.5,10,12.5,7.5,17.5,15,20,5}kHz，跳頻周期為5 ms。短時(shí)傅里葉變換采用長度為512的Hamming窗，采樣頻率設(shè)為10 MHz，對(duì)比分析本文算法和文獻(xiàn)[1]、文獻(xiàn)[10]算法在不同信噪比的檢測(cè)性能，信噪比范圍設(shè)置為-20～5 dB，間隔1 dB，并且進(jìn)行100次蒙特卡洛仿真，繪制檢測(cè)概率曲線，如圖4所示。

圖4 不同信噪比條件下三種算法的檢測(cè)概率

由圖4可知，三種算法的檢測(cè)概率均隨著信噪比的增加而增大，文獻(xiàn)[1]、文獻(xiàn)[10]以及本文算法實(shí)現(xiàn)90%以上檢測(cè)概率的信噪比分別為0 dB、0 dB、-8 dB，本文算法可以實(shí)現(xiàn)低信噪比環(huán)境下的高檢測(cè)概率。

在信噪比為-12 dB時(shí)，本文算法的檢測(cè)概率已達(dá)到60%；針對(duì)信噪比小于-12 dB時(shí)本文算法檢測(cè)概率偏低的問題進(jìn)行多次驗(yàn)證分析，發(fā)現(xiàn)當(dāng)信噪比低于-12 dB時(shí)，時(shí)頻矩陣中噪聲與信號(hào)的幅值差值較小，在利用K-means聚類時(shí)將部分信號(hào)判定為噪聲，進(jìn)而被截?cái)嗵幚?，使得時(shí)頻圖上跳頻信號(hào)部分缺失，影響了檢測(cè)概率；當(dāng)信噪比低于-16 dB時(shí)，信號(hào)已經(jīng)完全被噪聲淹沒，本文算法將不再適用。

3.2 本文算法的計(jì)算復(fù)雜度

取信噪比5 dB，其他仿真條件與3.1節(jié)相同，對(duì)本文算法和文獻(xiàn)[1]、文獻(xiàn)[10]算法在時(shí)頻分析、閾值選取以及截?cái)嗵幚磉^程中需要處理的數(shù)據(jù)數(shù)目進(jìn)行分析，并且進(jìn)行100次蒙特卡洛仿真，對(duì)比結(jié)果如表1所示。

表1 三種算法的計(jì)算復(fù)雜度對(duì)比

通過表1對(duì)比發(fā)現(xiàn)，文獻(xiàn)[10]算法基于對(duì)閾值選取方法的改進(jìn)，提高了檢測(cè)概率卻沒有降低計(jì)算量；文獻(xiàn)[1]算法進(jìn)行了STFT以及WVD兩種時(shí)頻變換，并且利用哈達(dá)瑪積運(yùn)算對(duì)時(shí)頻矩陣進(jìn)行組合，進(jìn)一步提高了較低信噪比條件下的檢測(cè)概率卻額外增加了計(jì)算量；本文算法在閾值選取以及截?cái)嗵幚磉^程中都基于第二節(jié)中的最大值進(jìn)行處理，使計(jì)算量降低，運(yùn)行時(shí)間比文獻(xiàn)[1]、文獻(xiàn)[10]算法降低了37.51%和12.34%，極大地提高了算法的檢測(cè)效率。

3.3 本文算法的有效性

取信噪比-8 dB，其他仿真條件與3.1節(jié)相同，對(duì)本文算法的有效性進(jìn)行仿真驗(yàn)證。由于時(shí)頻變換具有能量聚集性，使得信號(hào)的能量總是沿著瞬時(shí)頻率集中分布，即信號(hào)的能量集中分布在瞬時(shí)頻率周圍[20-22]，因此原始信號(hào)經(jīng)短時(shí)傅里葉變換后的時(shí)頻圖如圖5所示。

(a)二維時(shí)頻圖

由圖5可以看出，經(jīng)過短時(shí)傅里葉變換后，在信噪比較低的條件下，信號(hào)基本上被噪聲所淹沒，很難進(jìn)行準(zhǔn)確的判別檢測(cè)，圖5(a)中紅色標(biāo)記是淹沒在噪聲中的部分信號(hào)。

將仿真條件代入本文第2節(jié)，通過K-means聚類后選取自適應(yīng)閾值，對(duì)時(shí)頻矩陣進(jìn)行截?cái)嗵幚?，處理結(jié)果如圖6所示。由圖6可以看出，本文算法具有很好的去噪作用，在較低信噪比條件下依然可以去除噪聲，實(shí)現(xiàn)跳頻信號(hào)的檢測(cè)。

(a)二維時(shí)頻圖

4 結(jié)束語

針對(duì)低信噪比環(huán)境下跳頻信號(hào)的快速檢測(cè)問題，本文提出了一種新的閾值選取算法。該算法通過提取時(shí)頻矩陣每一頻率分量最大值進(jìn)行K-means聚類確定閾值和截?cái)嗵幚?，完成噪聲的去除，再依?jù)跳頻信號(hào)在駐留時(shí)間上的連續(xù)性與頻域上的跳變性實(shí)現(xiàn)對(duì)跳頻信號(hào)的檢測(cè)。仿真結(jié)果表明，本文算法計(jì)算量遠(yuǎn)低于同類算法，并且可以實(shí)現(xiàn)-8 dB環(huán)境下90%以上的檢測(cè)概率，為瞬變環(huán)境下的無人機(jī)跳頻信號(hào)的快速高效檢測(cè)提供了一種有效手段，具有一定的工程應(yīng)用價(jià)值。下一步的研究重點(diǎn)是搭建硬件平臺(tái)采集實(shí)測(cè)數(shù)據(jù)，并結(jié)合所提算法進(jìn)行分析驗(yàn)證。