张亚州，张海龙,3,4*，张 萌，王 杰，冶鑫晨，王万琼，李 嘉，杜 旭

(1. 中国科学院新疆天文台，新疆 乌鲁木齐 830011；2. 中国科学院大学，北京 100049；3. 中国科学院射电天文重点实验室，江苏 南京 210033；4. 国家天文科学数据中心，北京 100101)

随着射电天文技术的进步，射电望远镜向大口径及阵列方向发展。中国500 m口径球面射电望远镜(Five-hundred-meter Aperture Spherical radio Telescope, FAST)[1]、荷兰低频阵列(the Low-Frequency Array, LOFAR)[2]以及美国新墨西哥州的甚大天线阵(Very Large Array, VLA)[3]等观测设备具备高灵敏度和高分辨能力。望远镜在观测宇宙微弱天文信号的同时，也接收到干扰信号，强度比微弱天文信号强几个数量级的人为射电干扰信号，在射电天文数据处理过程中必须得到有效抑制。

射频干扰来源多种多样，文[4]将射频干扰分为外部射频干扰和内部射频干扰，外部射频干扰特指天文台址外产生的无线干扰，常见的有人造卫星[5](北斗和全球定位系统导航卫星等)、飞机、附近的基站信号以及广播信号等；内部射频干扰是指在天文台址内使用的电子设备，如各类数字终端设备、计算机、微波炉、咖啡机、荧光灯以及其他电子仪器等产生的射频干扰。

射电天文学发展过程中，射频干扰抑制一直是研究热点。文[6]提到天文信号通常为宽频带且在时间尺度上平滑变化，绝大部分射频干扰与天文信号相比有明显的差异。针对射频干扰抑制，文[7]提出了CUSUM(Cumulative Sum)方法，该方法简单高效，但依赖于单点采样，对噪声比较敏感。为了克服这一局限性，文[8]提出了基于计算相邻采样点的组合阈值表示算法SumThreshold，文[9]利用奇异值分解来解决巨型米波射电望远镜(Giant Metrewave Radio Telescope, GMRT)宽带射频干扰问题。这些算法本质上基于阈值实现，处理幅值远大于天文信号的射频干扰信号时效果较好，但对于特征与天文信号接近的射频干扰信号，阈值法可能误删天文信号，导致信息丢失，影响天文研究的准确性和有效性。

瞬时射频干扰是某些设备或系统正常运行时释放的射频信号。继电器、电动机和开关电源都会产生瞬时干扰。由于瞬时干扰信号是宽频、脉冲的，而且持续时间非常短，难以识别来源。文[10]提出了基于字典的瞬时射频干扰分类算法，使用隐马尔科夫模型将瞬态射频干扰识别为子序列，并能从观测数据中提取瞬时射频干扰序列。

1 射频干扰信号特征提取算法

1.1 信号周期计算

大部分射频干扰信号在时间尺度上呈现明显的周期波动。为了分析射频干扰信号的特征信息，我们需要对信号分窗处理，窗口大小与信号周期信息密切相关，计算信号周期是分析特征信息的重要前提。

信号周期计算算法有传统的倒谱法、自相关函数法、平均幅度差函数法和小波变换法等[11]。自相关函数法是一种非常有效的低信噪比信号周期估计算法[12]。基于自相关函数法，本文设计了一种自相关卷积滑动算法计算射频干扰信号的周期。算法过程为

(1)对一段长度为N的射频干扰信号振幅取绝对值；

(2)从射频干扰信号上截取一段长度为M的信号作为滑动窗口信号；

(3)滑动窗口与射频干扰信号从起始位置进行卷积，每次向后滑动d个采样点；

(4)重复步骤(3)，直至滑动窗口移动到射频干扰信号尾部；

(5)利用卷积后结果的峰值信息计算信号的周期。

图1 (a)k7ing-3545 kHz噪声信号时域图；(b)k7ing-3545 kHz噪声信号周期图

1.2 射频干扰信号特征提取

峰值是射频干扰信号重要的特征信息，能够刻画信号波动的变化细节。射频干扰信号特征提取过程的具体步骤为

(1)选取射频干扰信号窗口；

(2)对该窗口信号的振幅取绝对值；

(3)获取该信号窗口的峰值信息；

(4)对离散峰值信息进行平滑预处理；

(5)处理后的峰值信息作为射频干扰信号的特征信息。

算法第1步选取射频干扰信号窗口，对一段包含射频干扰的信号提取特征信息时，我们需要选取一段时间跨度尽可能短且至少包含一个完整射频干扰周期信息的信号段，根据1.1节提出的射频干扰信号周期计算方法确定信号周期。射频干扰信号相对于其他信号振幅较大，起始和终止位置明确，算法将信号振幅突增的位置作为窗口的起始点。射频干扰信号特征提取窗口大小与计算得到的周期相关，计算机显示器产生的干扰信号窗口选择如图2，窗口大小为D。第2步对信号的振幅求绝对值，以获取更多的峰值信息。第3步利用振幅特性，获取窗口的峰值信息。第4步对获取的峰值信号进行预处理，将相邻且振幅大小相差较小的值用均值代替，平滑峰值曲线，减少局部波动，缩小同类别信号差异，提高识别效率。图3(a)为射频干扰信号峰值信息折线图，波动趋势不平稳，不利于后续特征识别。平滑处理后射频干扰信号的峰值曲线如图3(b)，曲线相对平稳，能够缩小同类信号的差异。平滑处理后的峰值信息作为识别射频干扰信号的特征模板，可以进行相似度计算。

图2 射频干扰信号窗口选取

图3 (a)射频干扰信号的峰值曲线；(b)平滑后射频干扰信号的峰值曲线

2 射频干扰信号特征相似度算法

对未知信号进行特征提取后，计算与已知特征模板的相似程度，实现对未知射频干扰信号的识别和分类。算法采用分段打分策略计算两个序列的相似程度，每个序列根据相邻离散点可以分为多个区段，比较两个序列对应区段的趋势是否相同，进行记分，趋势相同加分，否则减分。

在相似度计算前，要确保这两个序列处于对齐状态，即确定两个序列的区段是否一一对应。待识别信号具体细节未知，提取峰值信息后，离散的峰值信息几乎不可能与特征模板对齐，即相位信息没有对齐，因此，峰值位置信息等都需要重新调整，为计算相似度做准备。

文[13]提出了一种测量时间序列相似性的方法，即动态时间规整，该方法能够比较不同长度的时间序列[14]，已广泛应用于语音信号处理。该算法计算得到的欧几里得距离越小，表明两种声音模式的相似度越高[15]。利用动态时间规整算法计算两个长度不同序列的欧几里得距离最小时，两个序列的相位点对应情况如图4。

图4 使用动态时间规整算法后两个序列的相位点对应情况

基于动态时间规整算法，本文设计了射频干扰信号特征相似度计算算法，具体过程为

(1)通过动态时间规整算法计算两组序列在最短欧几里得距离条件下离散点的对应情况；

(2)计算未知信号离散点的权重值；

(3)分段比较，如果对应段的变化趋势相同，相似度分数增加，否则相似度分数减小；

(4)累计各段的分数，总分数代表两个序列的相似程度。

射频干扰信号特征相似度计算算法第1步利用动态时间规整算法求得待比较序列在最小欧几里得距离下离散点的对应情况，为后续计算做准备。第2步计算被比较序列的离散点权值，计算方法为

(1)

其中，D为序列的离散点点数；S(t)为序列第t个离散点的幅值；W(t)为序列第t个离散点的权重，离散值越大，权重越大，在信号特征中的代表性越强。为了使相似度介于-1到1之间，第一个和最后一个的权重设置为0。

第3步进行序列相似度计算，采用分段打分累计的方式。通过动态时间规整算法计算每一区段的对应情况，采用分段比较方式，如图4。两个序列的区段有趋势相同、趋势不同和一对多3种情况，两个区段的趋势相同时加分，计算公式为

score=score+M[i] ，

(2)

其中，score是两个序列相似度的累积分数；M[i]是第i个区段的权重，大小为区段的两个离散点权重的平均值，即0.5×(W[i]+W[i+1])。两个区段的趋势不同时减分，计算公式为

score=score-M[i].

(3)

一对多情况不进行加减分。计算并累计各段的分数后得到两个序列相似度值。

3 实验结果与分析

本文选择ARRL官网提供的干扰源作为实验数据，选择7种射频干扰源进行交叉测试，基本信息如表1。为方便对射频干扰信号进行特征提取及后续识别测试，利用1.1节提出的射频干扰信号周期计算法得到各实验数据周期采样点数，以便下一步特征提取及识别测试。

表1 射频干扰源的基本信息

随机选取每个信号的10个周期作为测试数据，提取相应候选特征模板与其他信号进行交叉测试，选择与相同信号相似度大且与其他信号相似度较低的候选模板为最优模板。如图5(a)选取k7_2作为k7ing-3545 kHz的特征模板，(b)选取ks_9作为ks2am-streetlight的特征模板，(c)选取mo_4作为monitor的特征模板，(d)选取n6_10作为n6rce-sps-carrier的特征模板。图6(a)选取pl_8作为plc-4的特征模板，(b)选取au_7作为ausoth的特征模板，(c)选取ot_4作为18120oth25的特征模板。

图5 (a)选择k7ing-3545 kHz不同模板的相似度结果；(b)选择ks2am-streetlight不同模板的相似度结果；(c)选择monitor不同模板的相似度结果；(d)选择n6rce-sps-carrier不同模板的相似度结果

图6 (a)选择plc-4不同模板的相似度结果；(b)选择ausoth不同模板的相似度结果；(c)选择18120oth25不同模板的相似度结果

提取每个信号最优特征模板所在组的数据，结果如表2。相同信号的相似度在所在行最大，表明提取的射频干扰特征模板能够较正确地识别射频干扰来源。算法的相似度计算基于特征模板得出，权值计算依赖于特征模板，即计算两个序列的相似度时，权值的大小只与特征模板有关，因此与对角线互为对称的数值不同。同一信号在不同周期内存在差别，造成同一信号的相似度小于1。

表2 射频干扰信号相似度验证平均值结果

表2中个别数据偏高，如n6rce-sps-carrier行plc-4列的数值为0.726。通过分析得知，待比较的一组序列振幅波动趋势比较一致，导致相似度偏高，后续工作将继续优化算法，以提高精度。

4 总 结

本文设计了滑动卷积周期计算算法，利用卷积后的峰值区间完成了信号平均周期计算；基于设计的峰值提取算法实现了信号特征模板提取；基于动态时间规整算法和打分策略设计了信号特征识别算法，实现了未知信号的识别和分类。本文分别对从ARRL官网下载的射频干扰数据进行互相关与自相关计算，实验结果表明，信号生成的模板与原始信号的相似度明显高于其他信号，说明本文提出的算法可以有效提取特征信息并生成信号特征模板，且利用特征识别算法能对射频干扰进行正确分类。本文提出的方法可以对射频干扰进行细粒度识别和标记，有望为射频干扰特征识别和标记提供一种新的解决方案。

本文中算法实现代码以及实验数据已在码云仓库(https://gitee.com/zyazhou/rfi-feature-recognition.git)开源。