基于Simple Thresholding和CUSUM联合算法的L波段太阳射电流量可观测频段分析

杨 航，董 亮，何乐生*

(1. 云南大学，云南 昆明 650500；2. 中国科学院云南天文台，云南 昆明 650216)

太阳活动是空间天气的源头，日冕物质抛射(Coronal Mass Ejection, CME)[1]事件是影响日地空间环境最主要的因素，射电观测反映不同日冕高度的等离子体信息[2]，包括密度、磁场等。等离子体抛射时，其密度随着远离太阳而降低，随着日冕物质抛射远离太阳，其发射的无线电频率也随之降低。为此，太阳物理界一致认为地基太阳射电望远镜能观测到的无线电爆发频率越低，日冕物质抛射离地球越近。在不同频率上的太阳射电观测可以反演日冕高度的动态信息[3-5]。

同时，长期太阳射电流量监测发现，不同频段的流量变化与太阳活动周高度相关，是重要的空间天气观测指标。但另一方面，射电的观测频段与民用无线电使用频段重合，极易干扰观测数据。在射电望远镜建设之前的无线电环境测试与评估尤为重要，与脉冲星、活动星系核(Active Galactic Nucleus, AGN)光变等需要宽带观测不同[6-7]。由于高时间分辨率要求，太阳射电精密流量观测时对一些偶发性的干扰并不敏感，这个短时 “坏点” 可以通过后期与太阳光学等数据进行比较，排除爆发可能。

L波段太阳射电爆发是导航通信潜在的干扰因素，与其他空间天气事件不同，太阳爆发的强信号以光速传播，是作用最快的空间天气事件，从理论上来说，整个迎日面的导航接收机均有可能受到干扰[8-10]。文[11-13]提出了采用L波段多通道精密太阳射电望远镜同步监测的办法，实时监测L波段太阳射电流量，在爆发早期发出预警信号。

但是L波段存在大量无线电干扰，在建设望远镜前，我们需要评估MHz级的观测带内无线电干扰信号的分布规律，判断可观测太阳时间，对于长时间被干扰占用的信道进行排查，以防发生误报等情况。针对即将在云南天文台建设的L波段精密流量太阳射电望远镜(针对中高日冕进行研究和太阳射电爆发干扰导航通信预警，并作为廊坊太阳射电望远镜的备份)，本文介绍了对台址的无线电测试结果和观测频段遴选方法以及遴选介于北斗B1, B2和B3频点，全球定位系统L1和L2频点之间7个5 MHz无线电干扰较少的无线电通带。

1 系统组成

1.1 系统组成

测试系统由HL050天线、低噪声放大器、稳压电源、FSU26频谱仪和计算机组成，测试系统框图如图1。

系统框图中，由HL050天线将采集的信号通过低噪声放大器，传输到FSU26频谱仪转化成频谱，并把数据传输给计算机进行显示和存储。其中稳压电源为放大器提供稳定电压，保证数据的准确性，计算机通过LabView图形化编程，实现控制频段范围和极化方式，主要分为两部分：HL050天线频段控制和频谱仪控制程序[14]。

图1 测试系统Fig.1 Test system

1.2 测试设备及参数

HL050天线应用广泛，主要参数见表1，其优点是结构简单、频带宽、功率容量大、调整和使用方便。

频谱仪采用德国罗德施瓦茨公司的FSU26，工作频率20 Hz～26.5 GHz，平均噪声电平为-158 dBm，参数见表1。

表1 测试设备及参数Table 1 Test equipment and parameters

计算机利用LabView软件图形化编程实现，在LabView软件图形化编程中，可以调节测试的频段起始频率、极化方式、分辨率带宽、积分时间和执行次数，并将采集的信号保存在预设地址中，便于后续数据的处理分析。

1.3 监测过程

测试地点为云南省昆明市中国科学院云南天文台(25°1′48.538″N，102°48′14.180″E)，我们在东、西、南、北4个方向进行L波段的测试，系统频率范围为1～2 GHz，频谱仪扫频宽度设置为1 000 MHz。频谱仪显示带宽和分辨率带宽设置为300 kHz，扫描时间30 s。4个方向都测试了4次且时间超过8 h，共测试16次，数据均超过800组。具体测试日期及其组数如表2。

2 数据处理

本节主要介绍采集数据的格式、数据阈值遴选和本次最终信道的遴选过程。

表2 测试日期及组数(2022年)Table 2 Date of test and number of groups (in 2022)

2.1 数据介绍

本文采集的数据记录格式为TXT，单位为dBi，缩小数据的绝对数值，方便记录。由于整个测试过程设置频谱仪输出两条线，所以每一组检测数据分为两列，如图2。第1列为所采集数据的最大保持值，第2列为扫描时间内的平均值，为每次检测时判断测试系统的稳定效果。数据储存为非线性的数据格式，这是一种对数关系，底数为10，数学关系为log10(P2/P1)，其中P2/P1表示功率比。后续数据处理需要对数据进行线性转化[15-16]。

图2 测试系统实测图Fig.2 Actual measurement diagram of the test system

2.2 阈值遴选方案

针对即将在云南天文台建设的L波段精密流量太阳射电望远镜，遴选介于北斗 B1(1 559.052～1 591.788 MHz)，B2(1 166.220～1 217.370 MHz)和B3(1 250.618～1 286.432 MHz)频点，全球定位系统L1(1 574.397～1 576.443 MHz)和L2(1 226.577～1 228.623 MHz)频点之间7个5 MHz干扰较少的无线电通带。已知全球定位系统中L1频点在北斗B1频点范围内，且B2与L2频点之间相差不超过10 MHz，无法遴选两处相对洁净的5 MHz通带，B2右边遴选的通带与L2左边遴选的通带重合，一共7个通道(经过无线电环境监测带内无干扰信号)，全球定位系统和北斗系统的所有通信频点都穿插在监测频点之间[17-20]。

本次测试目的是从L波段中遴选洁净率高于95%的可用信道。遴选方式是通过绘制频谱瀑布图，直观地选出较为干净的信道，通过遍历观测数据超过阈值的占比确定频谱洁净率，当遴选的信道的洁净率都大于95%，说明该信道可用。处理流程如下：

(1)选取一个短时 “宁静” 的观测频段作为阈值计算数据样本；

(2)在上述数据样本中通过本文提出的阈值算法计算得到阈值；

(3)以此阈值遍历其余观测时间段数据，对比阈值，超过阈值时间和低于阈值时间进行对比，求得该时间段的干扰信号占比；

(4)利用算法计算得到的阈值对本文遴选的7个5 MHz干扰较少的无线电通带进行评估。

2.2.1 阈值计算方案

阈值分析法因实现简单、检测结果精度较高而广泛应用。代表性阈值分析法有CUSUM法[21]和Simple Thresholding法[22]。在检测某行(列)观测数据时，CUSUM算法通过估算累积样本的方差与平均值得到阈值，Simple Thresholding算法使用该行(列)的中位数作为阈值[23]。

本文首先监测取对数之后的数据，取对数不会改变数据的性质和相关关系，且压缩了变量的尺度，数据更加平稳，也消弱了模型的共线性、异方差性等，所以计算阈值时需要将数据线性化，公式为

(1)

样本方差是表示射频干扰的统计异常的重要估计值。由于绝对中位差(Median Absolute Deviation, MAD)和中位数是剔除异常值的最有效方法，因此，本文提出一种基于Simple Thresholding算法与CUSUM算法的改进阈值算法对台址的L波段无线电测试结果和观测频段进行遴选，如果所选取信道阈值水平超过Th，视为射频干扰。Th计算公式为

(2)

MAD=1.4286×med1≤i≤n{|xi-med(xi)|}，

(3)

(3)式中，n=2m时，med=0.5(xm+xm+1)；n=2m+1时，med=x(m+1)；MAD是Simple Thresholding算法使用的方差估计值，中位数是原始总体数组xi的中位数。然后用此中位数减去原始样本的每个元素，创建与原始元素大小相同的新修改样本。然后计算这个新样本的中位数，并用恒定的比例因子1.428 6相乘[24]，使该估计与预期的高斯分布一致。得到线性计算的阈值Th之后，需要将阈值取对数转换在监测的数据中计算通带的洁净率，通过

Th*=10-3log10(Th)

(4)

得到阈值Th*。整个阈值算法框图以及算法伪代码如图3。

图3 阈值算法框图及算法伪代码。(a)阈值算法框图；(b)阈值算法伪代码

2.2.2 信道选取与数据比较

首先选取一个短时 “宁静” 的观测频段作为阈值计算数据样本，从表2中选取此次检测组数最多的一次作为信道的初选，西向测试在第4次04月27日时共采集了2 219组数据，从西向第4次数据的频谱瀑布图中选取短时 “宁静” 的观测频段，频谱图如图4～图6； 其次在上述数据样本中通过阈值算法计算得到阈值，以此阈值遍历其余观测时段数据，对比阈值，超过阈值时间和低于阈值时间进行对比，求得该时段的干扰信号所占比例，对比结果如表3；最后利用算法计算得到的阈值对本文遴选的7个5 MHz无线电干扰较少的无线电通带进行评估。

图4 L波段频谱瀑布图Fig.4 L-band spectral waterfall

图4中横坐标为L波段1～2 GHz，纵坐标为数据采样的组数，竖轴为采样的信号强度，因为能较为直观地从图中看出本文所需频段的初步干扰情况，故此次L波段的所有测试均采用此种绘制方法。

图5 L波段频谱右视图Fig.5 L-band spectrum right view

图5为L波段频谱瀑布图的右视图，图中能清晰看到所有测试组数的信号强度，且能直观看到本文所需的频段的干扰情况，能为本文寻找短时 “宁静” 的观测频段确定初步筛选工作。图6为L波段频谱瀑布图的俯视图，图中标记北斗 B1, B2和B3频点及全球定位系统 L1和L2频点位置，为后续通带的遴选起辅助作用。

图6 L波段频谱俯视图Fig.6 L-band spectrum vertical view

通过阈值算法对数据阈值的确定，实现对整个L波段监测数据遴选的多个5 MHz无线电干扰较少的无线电通带中，最终确定介于北斗系统及全球定位系统频点附近的7个无线电通带。信道选取流程如图7。

本次测试主要是为了遴选介于北斗 B1, B2和B3频点，全球定位系统 L1和L2频点之间7个5 MHz无线电干扰较少的无线电通带，首先从图4～图6中L波段频谱图中初选出5 MHz干扰较少的无线电通带1 161～1 165 MHz频段作为初选阈值确定的可用信道，再根据上述算法确定阈值后，通过阈值遍历对比程序计算考察信道的洁净率，得到各组信道监测数据的洁净率如表3。从表3可以看出，在16次L波段测试中所遴选的5 MHz干扰较少的无线电通带1 161～1 165 MHz频段信道洁净率都在95%以上，符合太阳观测需求，所以本文选用此阈值为所遴选的7个5 MHz干扰较少的无线电通带的阈值。

本文结合图5～图7以及表3初步遴选介于北斗 B1, B2和B3频点以及全球定位系统 L1和L2频点之间7个5 MHz干扰较少的无线电通带如图8，从图8看到本文所选的通带具体的干扰情况：图中B2a(1 161～1 165 MHz)和B2b/L2a(1 221～1 225 MHz)为B2频点附近通带；B2b/L2a(1 221～1 225 MHz)和L2b(1 231～1 235 MHz)为L2频点附近通带；B3a(1 246～1 250 MHz)和B3b(1 291～1 295 MHz)为B3频点附近通带；L1a(1 551～1 555 MHz)和L1b(1 596～1 600 MHz)为B1/L1频段附近通带。利用上述所求阈值确定遴选的7个5 MHz干扰较少的无线电通带的洁净率如表4。

图7 遴选信道流程框图Fig.7 Block diagram of the selection channel flow

表3 1 161～1 165 MHz频段洁净率Table 3 Cleanliness rates in the 1 161-1 165 MHz band

图8 L波段频谱遴选通带Fig.8 L-band spectrum selection channel

表4 遴选通带洁净率Table 4 Selection belt cleanliness rate

由表4可知，北斗系统B1、全球定位系统L1频点遴选通带为1 551～1 555 MHz和1 596～1 600 MHz，洁净率为98.329%和98.301%；北斗系统B2频点遴选通带为1 161～1 165 MHz和1 221～1 225 MHz，洁净率为98.315%和98.335%；北斗系统B3频点遴选通带为1 246～1 250 MHz和1 291～1 295 MHz，洁净率为98.224%和97.650%；全球定位系统L2频点遴选通带为1 231～1 235 MHz，洁净率为98.260%，且所遴选通带的洁净率都大于95%以上，符合太阳观测需求。

3 总结与展望

针对即将在云南天文台建设的L波段精密流量太阳射电望远镜，本文提出一种基于Simple Thresholding算法和CUSUM算法的改进阈值算法，并对台址的L波段无线电测试和观测频段遴选出介于北斗 B1, B2和B3频点，全球定位系统L1和L2频点之间7个5 MHz干扰较少的无线电通带，分别为1 551～1 555 MHz, 1 596～1 600 MHz, 1 161～1 165 MHz, 1 221～1 225 MHz, 1 246～1 250 MHz, 1 291～1 295 MHz和1 231～1 235 MHz，其洁净率分别为98.329%, 98.301%, 98.315%, 98.335%, 98.224%, 97.650%和98.260%，均符合太阳观测需求，为下一步接收机的设计和信号处理提供了依据。