射电天文台址常驻电磁干扰检测与识别方法*

蔡明辉 刘 奇 王 玥 苏晓明 党振伟 张国福

(1中国科学院新疆天文台 乌鲁木齐 830011)(2新疆微波技术重点实验室 乌鲁木齐 830011)

1 引言

大口径射电望远镜具有极高的灵敏度和很宽的工作频率范围，因此观测时极易受到台址内外各类电磁信号的干扰，这些干扰主要来自人类生产生活所产生的各类无线电信号.国际电信联盟(International Telecommunication Union，ITU)所属的无线电通信部门协调全球的无线电频谱资源，在建议框架中制订了频率分配和运行规则的综合清单，包括电台允许的广播频段和功率、射电天文受保护的频段等.国际电联只为射电天文业务(Radio Astronomy Service，RAS)分配了非常有限的频谱资源.实际工作中，天文学家对宽带接收和记录数据的需求不断提升，于是RAS系统的工作带宽不断增加，而且射电天文业务是一种被动(接收)业务，很容易受到调频广播(Frequency Modulation，FM)、电视信号(TV)、公众移动通信(Mobile Global System for Mobile Communications，如Mobile GSM、Mobile Code Division Multiple Access(Mobile CDMA)、Mobile Long Term Evolution(Mobile LTE))、卫星定位系统(如Global Positioning System，GPS)和飞机导航通讯(如Inmarsat)等主动(发射)业务的影响，导致观测数据丢失和数据质量的降低[1-2].表1列举了主要的射频干扰(Radio Frequency Interference，RFI)来源.

表1 低频射电天文观测中部分主要RFI来源Table 1 Part of the main sources of RFI contamination at low radio frequencies

国内外学者对射电望远镜受到的电磁干扰开展了深入研究，Waterman[3]讨论了影响射电望远镜工作的各种电磁干扰源，其中卫星传输、电视发射塔、移动通信基站是影响射电观测最主要的射频干扰源，在实际射电天文观测过程中这些干扰信号使得干涉仪接收机安全饱和.大多数情况下，这些无线电设备的带外和杂散辐射超标是引起干扰的主要原因[4].Bolli等人分别采用了固定和移动监测手段，对撒丁岛射电望远镜台址附近的电磁干扰进行了监测，提出适当的连续性射频干扰监测对于确保射电望远镜的正常运行至关重要:在分配给RAS的频段内发生的射频干扰，必须连同(实验检测到的)有害性向国家无线电管理部门报告;根据对特定时间(如夜间或周末)预期“无干扰”的统计性评估，动态安排射电天文观测[5].在电磁干扰信号识别方面，王利云等人通过将频谱中的信号与历史信号进行相似性分析，判断是否为同一信号[6].但当历史信号未发射，但其信号频段被其他信号占用，例如与原信号调制方式一致的新发信号或其他信号的谐波、互调等，仅通过对信号包络的相关性计算判断两个信号是否为同一信号会造成误判.

综上所述，良好的无线电环境是射电天文观测科学产出的重要保证，为确保对望远镜台址各类电磁干扰的有效监测，目前新疆天文台已在110 m射电望远镜台址(110 m Qi Tai radio Telescope，QTT)安装部署了自动化、高可靠性电磁环境测试系统，实现了对射电望远镜周边0°-360°方位，水平极化、垂直极化电磁干扰的监测[7-8].但在电磁干扰信号的分析判断方面，仍需要占用技术人员较多时间与精力.

本文针对天文台址电磁环境频谱中的地面数字电视、公众移动通信等常驻大信号，提出了一种信号检测与识别方法.首先对常驻信号的业务类型、传播信道、调制方式进行分析，通过信噪分离、信号能量估值从电磁环境历史监测数据中提取移动通信、地面数字电视等常驻大信号的中心频率、极化方式、来波方向等特征，经对无线信号传播信道的衰落方式进行分析，提出信号能量模型服从的分布函数假设并采用K-S检验验证该假设，合理设置能量阈值，建立信号模板库;然后将实时频谱数据与信号模板库中的信号进行双门限能量检测与相关性判断，统计信号在监测周期内出现的次数及变化情况.经实际测试该方法可提高常驻电磁干扰检测与识别的准确率与时效性，为望远镜建设及运行阶段系统电磁兼容性设计、干扰缓解策略制定、台址无线电管理提供重要依据.

2 电磁环境测试系统

在射电望远镜台址安装部署的自动化电磁环境测试系统由测量天线(HL033、HL050)、射频接收模块(RF box)、天线云台(G5500)及控制器(Rotator controller)、升降装置、信号分析仪(R&S FSW13)、控制计算机(PC)、直流电源(DC power supply)及射频线缆(RF cable)、系统控制与监控线(Control&monitor)、云台控制线(Pan tilt control)、串口线(RS232)、网线(LAN)等组成，如图1所示.测量天线采用对数周期天线，工作带宽分别覆盖100 MHz-2 GHz及2-12 GHz频段，天线云台满足360°方向水平转动及水平、垂直极化切换需求.将信号分析仪、控制计算机、云台控制器等安装于屏蔽机柜(Shield cabinet)内，避免自身电磁干扰影响测量结果[8].

图1 测试系统硬件框图[8]Fig.1 Block diagram of test system[8]

电磁环境监测系统的核心部分是射频接收模块，采用两个标准噪声源(Noise source 1，Noise source 2)对0.1-2 GHz和2-12 GHz微波链路系统增益和系统噪声温度进行测量，选用性能优异的低噪声放大器(LNA1、LNA2、LNA3、LNA4)确保监测系统拥有较好的灵敏度，该模块链路结构如图2所示(其中S1、S2、S3、S4、S5为射频开关，ATT1、ATT2、ATT3为衰减器，Hpf为高通滤波器)[8].对系统整体性能进行测试时，在重点测试频段100 MHz-6 GHz要求噪声系数小于2.5 dB、系统增益大于40 dB、动态范围达到-127-37 dBm即90 dB.

图2 射频接收模块链路[8]Fig.2 RF receiver link[8]

3 常驻信号分析

经对天文台址电磁环境测试频谱序列进行统计，电视发射塔、移动通信基站、卫星传输等是影响射电观测最主要的固定干扰源，其所发射的信号占总信号数量60%以上而且射频特性稳定.为了有效地对这些常驻信号进行检测与识别，需要对其业务类型、传播信道进行分析，掌握信号的发射频率、占用带宽、场强的变化规律等特征信息.

3.1 常驻信号业务类型

移动通信业务:3大运营商2G、3G、4G相关业务均已覆盖，工作频段在800-2700 MHz范围内，多种制式并存，调制方式主要为高斯最小频移键控、正交相移键控、8移相键控、正交幅度调制等[9].通常2G网络发射功率为60 W，3G网络发射功率为40 W，4G网络发射功率为20 W;对射电天文观测影响较大的为2G业务，对于山区为达到信号覆盖范围更广，存在增大功率发射情况.随着5G网络建设的加快，后续密切关注广电与移动共建共享的700 MHz、C波段的5G通信信号发射情况.

广播电视业务:地面电视信号频率范围跨度广，通常在40-80 MHz和160-230 MHz的甚高频以及470-960 MHz的超高频频带.经长期频谱监测结合广电总局于2020年发布的、关于按规划关停地面模拟电视有关工作安排的通知，要求各地于2020年12月31日前完成关停中央、省、市、县地面模拟电视信号工作，当前在台址所监测到的电视信号均为数字信号，频段分布在400-800 MHz频段，采用时域同步正交频分复用技术，支持多种调制方式:64种符号、32种符号、16种符号、4种符号的正交振幅调制[9].其辐射源特点为功率较大，覆盖距离远(几十公里至几百公里)，对射电天文低频观测影响较大.

卫星通信业务:频谱监测数据显示卫星通信业务主要包括卫星电话和导航业务，电磁干扰信号特征及影响分析如下:卫星电话主要有海事卫星(Inmarsat)和 欧 星(Thuraya)，工 作 频 段1525-1559 MHz，来自120°-180°方向;铱星(Inridium)，工作频段1610-1626.5 MHz.卫星导航电磁干扰为下行通信业务，其信号特点为覆盖范围极广，为右旋圆极化，信号来自120°--240°方向，工作频段覆盖1559-1604 MHz，包括美国全球定位系统(GPS)、俄罗斯格洛纳斯(Global Navigation Satellite System，GLONASS)、欧洲伽利略(GALILEO)、中国北斗导航卫星.

3.2 传播信道分析

无线信号的发射功率决定了信号的电场强度，直接关系着信号的覆盖范围、覆盖质量、可靠性等，通常常驻信号的发射功率按照标准配置基本保持恒定[9-10].

信号在传播过程中遇到山体、树木等会产生反射、散射等现象，都会在链路的发射端和接收端之间形成多条传输路径，而这些传输路径对无线通信会造成严重的影响，简称为多径效应.经过不同路径到达接收端的信号具备不同的传输时延，相位旋转和幅度增益，因而在接收端接收到的信号是多条路径的复合信号，当上述复合信号的各个多径信号分量相抵消时，就会造成严重的信号衰落.

在台址区域中无线电信号的衰落方式主要为大尺度衰落，可以分为:(1)随着传播距离增加导致信号衰减的路径损耗;(2)由于障碍物遮挡、地形特性导致的电磁波阴影区域中，理想路径损耗模型产生较大变化的阴影衰落.大尺度衰落现象可以用对数正态模型来描述，这一模型的信号幅值遵循正态分布[11-12].

3.3 信号识别

调制信号的调制样式是区分不同通信信号的重要特征.近年来国内外学者对调制样式的自动识别技术进行了大量的研究，处理调制类型识别问题一般采用两种方法:判决理论法和统计模式识别法.判决理论方法根据信号的统计特性，得到检验统计量，然后与一个合适的门限进行比较.但在实际应用中存在识别过程中要求的参数太多，表达式计算复杂且难于处理等问题.统计模式识别算法是从已知信号样本中提取特征参数并给定误差对信号进行识别，该方法在低信噪比的情况下难以提取信号特征参数，识别能力低[13].

综上所述，常驻信号的发射频率和调制方式固定，发射功率按照标准配置基本保持恒定，传播信道存在大尺度衰落现象，但信号幅值遵循正态分布.因此本文针对实时频谱序列提出了一种常驻干扰信号检测及识别方法，处理流程如下:首先，根据历史频谱数据提取信号频率、能量、频谱等特征信息，存储于信号模板库;然后，将实时频谱数据与信号模板库中的信号进行双门限能量检测与相关性判断，统计信号在监测周期内出现的次数，实现台址常驻干扰信号的检测与识别.

4 信号检测与识别方法

实现频谱序列中常驻信号的检测与识别需解决以下信号处理问题:(1)准确提取监测频谱序列中的常驻信号频谱，以便进一步获得信号特征;(2)根据信号的发射和传播特性，获得信号能量的分布函数，合理设定能量阈值并结合其他特征信息建立信号模板库;(3)选取信号相关性识别方法，分析频谱序列中干扰信号与信号模板库中信号的相似性，实现实时监测频谱序列中常驻干扰信号的快速检测和识别.

4.1 信号提取

依据无线电频谱监测中应用广泛的频谱噪声统计法[14]及无线电业务频率划分[15]确定常驻信号的起始、截止频率，进而从频谱序列中提取出2维频谱信号P{f，E}，其中f表示信号频率信息，E表示幅度信息.由于电磁环境监测依据测试频段内业务的特点进行分段测试分段存储，需要将信号P与所在频段中的存储序列进行匹配，以便提取出信号对应的频率点分布和幅度信息:由于f采取等间隔采样与存储顺序α成线性关系，采用简单线性回归方程来描述:

其中，β0是截距、β1是斜率(即回归系数)、ϵ是误差项.定义n为频率点列表的长度，将信号P所在频段的频率点列表及次序写成矩阵形式:

4.2 信号能量阈值估计

定义信号频谱能量:

其中i为信号样本在频谱中存储的序列，E i为幅度采样值，G为信号能量估值.依据一轮(在水平、垂直两个极化下，0°-360°方位间隔60°轮巡)监测过程中信号能量估值的最大值为信号能量值，对应的极化及方位为该信号的极化方式、来波方向.

根据§3.2分析，常驻信号的发射功率按照标准配置基本保持恒定，信号传播过程中衰落方式为大尺度衰落，信号能量遵循正态分布.严格意义上应当对以上假设进行检验，由于信号能量总体的分布形态未知，需要使用与总体分布统计无关的假设方法即非参数方法来检验信号能量值是否满足正态分布.K-S检验是一种极其重要的拟合优度非参数检验方法，它利用样本数据推断总体是否服从某一理论分布，具体是将样本的实际累计分布和假设的理论分布进行对比，计算两个分布函数之间距离的最大值D.若D值在要求对应的置信区间内，就可以判定这个样本服从这个假设的分布，否则认定样本与假设之间具有显著性差异.

采集不同时间段内信号能量值的多轮监测结果，利用K-S检验能量样本观测值是否近似服从正态分布，检验统计为:

其中Fζ(z)为信号频谱能量值的经验分布函数，G(z)为指定的分布函数.其原假设:信号频谱能量值服从正态分布.当D＞D(ζ，δ)则拒绝原假设，否则接受原假设.式中ζ、δ分别表示样本数量及检验的显著性水平，著名的统计学家R.A Fisher建议δ取0.05为宜[16]，D(ζ，δ)由K-S临界值表在指定的ζ、δ给出.

能量阈值λ1、λ2(λ1＜λ2)根据信号频谱能量所服从的正态函数合理设定，即λ1、λ2取正态分布的某一上下分位点.

4.3 信号模板库建立

根据业务类型对常驻信号进行分类，从历史频谱数据中提取完整的信号信息存入信号模板库.

目前常驻信号分为地面数字电视信号、移动通信信号、卫星业务信号(导航、卫星电话)3类.

信号特征信息包含特征信息表和数值表，特征信息表包括:信号的极化方式(Polarization)、来波方向(Azimuth)、中心频率(Center frequency)、能量阈值(Signal Threshold);数值表为信号频率信息及对应场强数值表(f&Elist).其中中心频率为起始频率和截止频率的均值，极化方式、来波方向、能量阈值依据§4.2所述确定;频率点及对应场强数值表按照§4.1中的信号提取方式获得，为进一步的信号相似性分析提供支撑.信号模板库结构如图3所示.

图3 常驻信号模板库结构Fig.3 The structure of resident signal template library

4.4 双门限能量检测

在电磁频谱检测方法中，能量检测算法具有复杂度低、检测速度快等特点，是一种可靠的盲源检测算法，在实际工程中应用广泛.双门限能量检测可解决信号能量值不确定的问题，提高判决结果的可靠性.双门限能量检测根据信号能量阈值λ1，λ2确定判决门限，然后将G与判决门限进行比较，若λ1≤G≤λ2，则判决假设H0成立:常驻信号存在;若G＜λ1，则判决假设H1成立:信道空闲，信号不存在;若G＞λ2，则判决假设H2成立:发出能量告警，信号能量值过大，可能会导致接收机饱和.

则能量检测的三元假设模型为:

4.5 相关性识别

由于常驻信号的调制样式保持不变，可以通过对信号包络的相关性识别判决实时采样信号与模板库中的信号是否为同一信号.皮尔逊系数广泛用于度量两个变量之间的相关程度，其优点在于:(1)能够区分正相关与负相关;(2)对相关程度的变化敏感.通过计算模板信号与实时频谱信号的皮尔逊系数来判断是否属同一信号.皮尔逊系数可表示为:

其中j为信号样本在频谱中存储的序列，ψ、ω分别为信号样本在频谱中存储的起始和截止次序，X j、Y j分别为信号采样值及对应的模板值，、分别为采样值均值、模板值均值.

皮尔逊系数的绝对值|γ|∈[0，1]，|γ|的值越趋近于1，表明信号X与信号Y相似程度越高;|γ|的值越趋近于0，表明信号X与信号Y相似度程度越低.通过确定相似阈值，若γ大于阈值，认为两个信号相关，反之则不相关.相关阈值的选取影响信号分析的精度，需要合理选取.

4.6 算法流程

根据历史监测数据确定常驻信号的起止频率并从频谱序列中提取出2维频谱信号.计算信号能量估值，其中的最大值为信号能量值，信号能量值的极化及方位为该常驻信号的极化方式、来波方向.根据能量值的分布函数确定信号能量阈值，结合上述信号特征建立常驻信号模板.读取实时频谱，将频谱特征与信号模板库中的特征信息表进行匹配.首先进行双门限能量检测，若频谱信号能量低于下门限时判断常驻信号消失(未发射)，高于上门限时发出能量告警(功率过高可能会导致监测系统饱和).当信号能量在能量阈值范围内时对信号包络进行相关性识别，若频谱信号的调制样式与模板信号一致则为常驻信号，算法流程如图4所示.

图4 常驻电磁干扰检测与识别算法流程Fig.4 Algorithm flow of resident EMI detection and identification

5 实例分析

本节拟从射电天文台址电磁环境历史实测数据中提取出常驻信号特征建立信号模板库，将实时频谱数据与模板库中的信号特征进行比对，实现信号检测与识别为例说明本方法的科学合理性.

台址电磁环境频谱数据按照极化、方向的不同分频段保存[17].经分析某地面电视信号频段为670-678 MHz、移动GSM900下行信号频段为935-954 MHz，从0.4-2 GHz频段监测数据中提取出以上信号，该频段共有53334个频率点即n=53334，结合频段频率点列表代入(2)式，得:

将电视信号、通信信号的起止频率代入(2)式及(6)式中，则电视信号、移动信号在频谱中存储的序列分别为9001、9002、···、9268和17834、17835、···、18468.根据(3)式计算的能量估值如表2所示(SN表示信号能量估值序列对应G1，···，G12;POL表示接收天线的水平和垂直极化分别对应H、V;AZ表示接收天线的指向方位对应0°，···，300°;TV代表电视信号，COM代表移动通信信号，NRG EST为能量估值.).

表2 不同极化和方向下的信号能量估值Table 2 Estimation of signal energy under different polarization and direction

表2中电视信号的能量值即能量估值最大值为105.33994，移动通信信号的能量值为191.43665，对应的极化方式、来波方向(DOA)如表3所示.

表3 电视与通信信号特征Table 3 Characteristics of energy value of TV and COM signals

选取连续4 d的历史频谱数据，即单一极化、单一方向下共计38组信号样本值.再由(3)式求得电视信号、移动信号的38组能量值如表4所示.

使用K-S检验方法，在显著性水平0.05下，表4中的地面电视信号能量近似服从均值109.9820、标准差为5.7927的正态分布;移动通信信号能量近似服从均值197.94、标准差为5.14的正态分布.取正态分布的上0.01分位点、下0.01分位点为能量阈值的上下限，则电视信号λ1、λ2分别为96.5、123.45;移动通信λ1、λ2分别为185.9826、209.8974.信号能量分布函数(x表示信号能量分布，F(x)表示信号能量分布函数，均为无量纲)如图5和图6所示.

表4 电视、通信信号能量值统计Table 4 Statistics of energy value of TV and COM signals

图5 电视信号能量分布函数Fig.5 Energy distribution function of TV signal

图6 通信信号能量分布函数Fig.6 Energy distribution function of COM signal

建立信号模板库，信号的特征信息表如表5所示(SIG表示不同的常驻信号对应TV和COM;fc表示信号的中心频率;Lower NRG THR和Upper NRG THR分别对应信号能量阈值的下限和上限.)，数值表取与统计值均值最邻近的频率/场强2维数组，即电视信号取序列号20(110.0617)，移动信号取序列号10(197.9639)所对应的P{f，E}.

表5 电视、通信信号特征信息表Table 5 Characteristic information table of TV and COM signals

使用该方法对射电天文台址电磁环境频谱进行连续一周测试统计，选取样本信号一天测试值即19组能量值进行举例，如表6所示.

表6 电视、通信信号能量值Table 6 Energy value of TV and COM signals

常驻信号能量值均在信号模板库能量阈值范围内，进行下一步相关性识别.经大量试验，信号样本与信号模板高度相关，相关性系数下限可取0.6，信号相关性如图7和图8所示.

图7 电视信号相关性Fig.7 Correlation coefficient of TV signal

图8 通信信号相关性Fig.8 Correlation coefficient of COM signal

对信号调制样式进行相关性识别，样本信号19组相关系数如表7所示.相关性系数均大于0.6，判断为常驻信号，与人工判断吻合.

表7 电视、通信信号相关系数Table 7 Correlation coefficient of TV and COM signals

对台址内所有常驻电磁干扰连续监测一个月的统计结果为:存在3个电视信号凌晨期间停发(OffTime)，其他时间持续发射(Run Time)，频道占用度(Channel Occupancy)接近80%，如表8所示;其余电视信号、移动通信信号及卫星业务信号均为连续发射，频道占用度为100%，统计结果与人工判断吻合.

表8 间断发射电视信号统计结果Table 8 Statistical results of discontinuous TV signals

6 结论

为了提高射电天文台址电磁环境监测频谱分析的时效性，本文提出了一种常驻电磁干扰检测与识别方法.经实际测试，该方法可提高信号判别的准确率，实现监测频谱中常驻干扰信号的自动识别与统计，为台址干扰缓解策略制定提供重要依据:根据无线电宁静区保护办法，在无线电主管部门的牵头下协调用频单位产生干扰的发射设备，根据实际情况调整其发射功率，方位、俯仰角等.也可根据对特定时间(如夜间或周末)预期“无干扰”的统计性评估，动态安排射电天文观测.后续工作安排:(1)随着5G网络建设的加快，密切关注台址内常驻信号的变化情况，实时更新常驻信号模板库;(2)统计信号能量随时间慢衰落的变化特征;(3)探索短时猝发干扰信号的识别方法.