一种基于K-means聚类的跳频信号快速检测方法*

(火箭军工程大学 导弹工程学院，西安 710025)

0 引 言

无人机产业的兴起为社会提供了诸多便利，但是“黑飞”无人机亦带来了一系列的安全隐患，进一步引起了国内外对反无人机技术的极大关注。跳频通信因具有较强的抗干扰、抗截获和易组网等能力[1]被广泛应用于无人机的信息传输，而跳频信号的检测是实现对无人机信号估计与干扰反制的重要前提。

目前，针对跳频信号的提取检测方法众多，可分为非盲检测算法和盲检测算法两类。非盲检测算法是指预先知道跳频信号部分参数，通过已知参数反映出跳频信号的某些特征，从而实现对跳频信号的检测的一类算法。文献[2]采用单跳自相关技术进行预检测处理，在跳速预知时，算法具有优越的抗噪性能。文献[3]提出基于多跳自相关的跳频信号检测方法，在已知跳频信号周期的前提下计算出跳频信号的自相关函数，通过判别自相关函数在不同时延情况下的值来实现跳频信号的检测。该算法计算量小，易于工程实现。文献[4-6]提出基于辐射计的信道化检测方法，跳频信道数已知时利用检测模型对各信道并行分析处理，根据各信道的检测结果判断跳频信号的有无，检测性能良好，但计算复杂。文献[7-8]提出的对跳频信号进行分层处理的检测模型，使得利用信道化检测时接收信道中的无效信道数量减少，降低了模型计算复杂度。文献[9]将接收信号进行主成分分析，降低数据维数，再进行数字信道化处理，通过短时能量对消实现跳频信号检测，降低了计算复杂度。以上方法均属于非盲检测，但实际战场环境中难以获取先验信息，同时也无法实现跳频信号的非盲检测。

盲检测类算法是指在不需跳频信号任何参数为先验信息的情况下对其进行检测的一类算法。在盲检测算法中，时频分析可以清晰地反映跳频信号的时频特性，因而被广泛应用。文献[1]利用短时傅里叶变换和 Wigner-Ville分布组合时频变换获取时频矩阵，去除背景噪声后进行二值化处理，利用连通域标记算法对同一信号像素点进行标记并提取信息，通过聚类算法实现信号的识别检测。文献[10]提出了一种基于局部自适应阈值的跳频信号检测方法，可以有效地对时频矩阵进行噪声处理，但抗噪性能较差，无法实现低信噪比条件下的跳频信号检测。文献[11-12]将时频图看作二维图像，根据图像边缘提取检测算子，有效去除噪声，但计算复杂度高，难以实时处理。文献[13]对跳频信号进行稀疏重构，利用形态学滤波去除噪声，算法抗噪性好但复杂度高，难以工程实现。文献[14]提出了一种基于滑动相关与小波变换结合的跳频信号检测方法，可较好地克服信道衰落和传播损耗的影响，适用于复杂电磁环境，但是检测效果受限于小波基的选取。以上算法均面临着计算复杂度高或者抗噪性差的问题，难以应用于复杂的战场环境。

针对低信噪比环境下跳频信号快速检测的需求，本文根据时频矩阵中随机噪声幅值在一定信噪比条件下通常弱于跳频信号且噪声点区域数据总是大于信号区域数据点数[1,12,15]的特点，提出了一种基于时频矩阵每一频率分量最大值进行K-means聚类的阈值方法，极大地减少了阈值选取以及数据截断的计算复杂度，并且实现了低信噪比条件下的跳频信号检测。

1 信号预处理

1.1 信号模型

跳频信号的数学表达式为

(1)

式中：k≥1，k为整数；Th为跳频周期；A为信号幅度;θ为初始相位;fk为第k个时隙的跳频频率;rectTh表示矩形窗，且满足

(2)

由于跳频频带较宽，故接收机接收的信号中往往掺杂着高斯白噪声，因此实际环境中跳频信号的接收模型表达式为

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

(3)

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

1.2 时频变换

目前时频变换方法主要有短时傅里叶变换[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 计算复杂度低，且无交叉项影响；WVD 交叉项严重；PWVD 和 SPWVD 通过增加窗函数，对交叉项有一定的抑制作用，但计算复杂，难以工程实现。基于工程实际的需求，本文选取STFT作为时频分析工具。

(4)

式中：ω(t,τ)表示窗函数，τ表示时延量。STFT(t,f)表示信号x(t)时间和频率的二维分布。

信号经过短时傅里叶变换，假设产生一个M×N维矩阵以及对应的时频图，如图1所示。

(a)时频矩阵部分示意图

在时频矩阵和二维时频图中，横向表示时间，纵向表示归一化频率。时频矩阵的每一行元素表示信号在时间窗内经过傅里叶变换产生的对应频率分量的幅值。图1(a)中红色、黄色分别表示信号以及信号的能量泄露，绿色部分表示噪声。由图1(a)可以看出，在一定信噪比条件下，时频矩阵中噪声的数量远远大于信号，分布广泛且幅值较小。基于此，本文通过K-means算法实现对噪声与信号的聚类区分，并设置阈值滤除噪声及信号泄露分量，进而基于驻留时间实现对跳频信号的检测。

2 基于K-means聚类的跳频信号检测

2.1 K-means聚类原理

聚类分割是把给定的样本集合按照某种准则分割成若干个不相交的子集，每个子集中的样本之间相似性较大，不同子集样本之间的相似性较小。K-means聚类由于其实现简单、聚类效果较好而被广泛应用。

K-means聚类算法采用距离作为相似性的评价指标，即认为两个对象的距离越近相似度越大，距离越远相似度越小。对于给定的样本集，按照样本之间的距离大小，将样本集划分为k个集合。让集合内的点尽量紧密地连在一起，而让集合间的距离尽量地大，集合内的相似性越大，集合间的差别越大，则聚类效果越好。

假设将样本集划分为k个集合(C1,C2,C3,…Ck)，则K-means的误差平方和准则函数表达式为

(5)

式中：ui是集合Ci的均值向量，有时也称为质心，其表达式为

(6)

式中：Ni是第i集合Ci的样本数目。在K-means聚类过程中，E度量了所有样本子集的总误差平方和。对于样本集的不同分类，导致不同Ci中的不同的样本子集以及对应的均值ui，从而会得到不同的E值，而当E最小时聚类效果最佳。

2.2 基于K-means聚类的阈值选取与数据截断

假设时频矩阵阈值为T，将矩阵中大于阈值的元素保留，小于阈值的置0，则有

(7)

式中：i和j分别表示时频矩阵的第i行和第j列，且1≤i≤M,1≤j≤N。

由时频矩阵可以看出，跳频信号具有时频聚集性，并且幅值远远高于噪声。因此，可以认为当时频矩阵中纵轴的归一化频率有信号时，其最大值为信号幅值；而没有信号时，因为均为噪声，则所有元素均远小于信号幅值。基于此，选取时频矩阵中每一行的最大值，组成一个新的向量数组作为K-means聚类的向量集，极大地减少了聚类所产生的迭代次数。

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

(8)

式中：Kmeans表示用于K-means聚类的含有M个元素的向量集，STFT(i,max)表示时频矩阵中第i行的最大值。

对向量集样本进行聚类，步骤如下：

Step1 把Kmeans中的M个样本数据分成k个聚类的初始划分，计算每个聚类的均值u1,u2,…uk和E。在本文中需要实现信号的检测，故取k=2，即将样本STFT(i,max)聚为噪声和信号两类。

Step2 从数据集中随机选择k个数据点作为质心。

Step3 对数据集中每一个点，计算其与每一个质心的距离，并划分到距离最近的质心所属集合。

Step4 把所有数据归为k个集合，重新计算每个集合的质心以及ui的值，并修改E。

Step5 若连续迭代多次E不变则停止，否则转到Step 3。

样本STFT(i,max)聚类前后如图2所示。

(a)样本数据聚类前

由图2(a)可以看出，样本数据在聚类前就有明显的分类迹象，下方密集数据表示的是噪声，上方零星分布的是跳频信号，通过聚类进行准确判别，可以认为图2(b)中的蓝色部分为噪声部分，红色为跳频信号。将聚类后的两类数据分别存放于矩阵Noise和矩阵Sig中，因为将Noise矩阵中的所有元素默认为噪声，并且所有元素均为每一行的最大值，故取Noise矩阵中最大值Noise(max)作为阈值，时频矩阵中所有小于等于Noise(max)的值均判别为噪声，并置0，则式(7)可重新表示为

(9)

由于STFT(i,max)表示每一行中的最大值，当最大值小于阈值时，其所在行的每一个元素均小于阈值，表示该行没有跳频信号，全部置0；当最大值大于阈值时，对该行的每一个元素同阈值进行比较，小于阈值的依然判定为噪声并置0，大于阈值的元素保留，则式(9)可以进一步改写为

(10)

通过时频矩阵以及二维时频图可以发现，跳频信号对应的频率分量在纵轴频率域占比很小。因此，基于式(10)可以对时频矩阵中所有元素快速地判别截断，再次大幅度减少了计算量，提高了检测效率。

截断处理会使跳频信号的部分边缘丢失，但是由图1(b)不难发现，在一定的信噪比条件下，跳频信号依然保留绝大部分的有用信息，即存在一定的驻留时间，在时频矩阵中就表现为时间维度上一段幅值不为0的连续序列，因此可以根据时域上的连续性来区别信号和噪声。当非0序列有一定长度并且在频域上依次变化时，则认为接收的信号中可能包含有跳频信号，从而实现对跳频信号的检测。

2.3 算法流程

通过上述的理论分析，初步认为本文所提算法在低信噪比条件下实现跳频信号的快速检测是可行的。算法流程如下：

Step1 采集一段足够长时间的数据，计算短时傅里叶变换，得到时频矩阵。

Step2 取时频矩阵每一行的最大值组成样本集Kmeans。

Step3 通过K-means聚类算法将样本集中的元素划分为噪声和信号两类，并将各自样本数据分别存放在矩阵Noise和Sig中。

Step4 选取矩阵Noise中元素最大值作为阈值。

Step5 根据式(10)，对时频矩阵进行截断处理，滤除噪声，得到新的时频矩阵。

Step6 根据跳频信号在驻留时间内的连续性，对新的时频矩阵进行检测。若非0序列均有一定的长度且在频率上依次出现，则认为存在跳频信号。

跳频信号检测算法流程如图3所示。

图3 跳频信号检测算法流程图

3 仿真实验与分析

3.1 本文算法在低信噪比条件下的检测概率

根据式(2)产生一段跳频信号，包含8个跳频周期，跳频频率集为{2.5,10,12.5,7.5,17.5,15,20,5}kHz，跳频周期为5 ms。短时傅里叶变换采用长度为512的Hamming窗，采样频率设为10 MHz，对比分析本文算法和文献[1]、文献[10]算法在不同信噪比的检测性能，信噪比范围设置为-20～5 dB，间隔1 dB，并且进行100次蒙特卡洛仿真，绘制检测概率曲线，如图4所示。

图4 不同信噪比条件下三种算法的检测概率

由图4可知，三种算法的检测概率均随着信噪比的增加而增大，文献[1]、文献[10]以及本文算法实现90%以上检测概率的信噪比分别为0 dB、0 dB、-8 dB，本文算法可以实现低信噪比环境下的高检测概率。

在信噪比为-12 dB时，本文算法的检测概率已达到60%；针对信噪比小于-12 dB时本文算法检测概率偏低的问题进行多次验证分析，发现当信噪比低于-12 dB时，时频矩阵中噪声与信号的幅值差值较小，在利用K-means聚类时将部分信号判定为噪声，进而被截断处理，使得时频图上跳频信号部分缺失，影响了检测概率；当信噪比低于-16 dB时，信号已经完全被噪声淹没，本文算法将不再适用。

3.2 本文算法的计算复杂度

取信噪比5 dB，其他仿真条件与3.1节相同，对本文算法和文献[1]、文献[10]算法在时频分析、阈值选取以及截断处理过程中需要处理的数据数目进行分析，并且进行100次蒙特卡洛仿真，对比结果如表1所示。

表1 三种算法的计算复杂度对比

通过表1对比发现，文献[10]算法基于对阈值选取方法的改进，提高了检测概率却没有降低计算量；文献[1]算法进行了STFT以及WVD两种时频变换，并且利用哈达玛积运算对时频矩阵进行组合，进一步提高了较低信噪比条件下的检测概率却额外增加了计算量；本文算法在阈值选取以及截断处理过程中都基于第二节中的最大值进行处理，使计算量降低，运行时间比文献[1]、文献[10]算法降低了37.51%和12.34%，极大地提高了算法的检测效率。

3.3 本文算法的有效性

取信噪比-8 dB，其他仿真条件与3.1节相同，对本文算法的有效性进行仿真验证。由于时频变换具有能量聚集性，使得信号的能量总是沿着瞬时频率集中分布，即信号的能量集中分布在瞬时频率周围[20-22]，因此原始信号经短时傅里叶变换后的时频图如图5所示。

(a)二维时频图

由图5可以看出，经过短时傅里叶变换后，在信噪比较低的条件下，信号基本上被噪声所淹没，很难进行准确的判别检测，图5(a)中红色标记是淹没在噪声中的部分信号。

将仿真条件代入本文第2节，通过K-means聚类后选取自适应阈值，对时频矩阵进行截断处理，处理结果如图6所示。由图6可以看出，本文算法具有很好的去噪作用，在较低信噪比条件下依然可以去除噪声，实现跳频信号的检测。

(a)二维时频图

4 结束语

针对低信噪比环境下跳频信号的快速检测问题，本文提出了一种新的阈值选取算法。该算法通过提取时频矩阵每一频率分量最大值进行K-means聚类确定阈值和截断处理，完成噪声的去除，再依据跳频信号在驻留时间上的连续性与频域上的跳变性实现对跳频信号的检测。仿真结果表明，本文算法计算量远低于同类算法，并且可以实现-8 dB环境下90%以上的检测概率，为瞬变环境下的无人机跳频信号的快速高效检测提供了一种有效手段，具有一定的工程应用价值。下一步的研究重点是搭建硬件平台采集实测数据，并结合所提算法进行分析验证。