CN108111240A - 一种射频干扰的数据处理方法及装置 - Google Patents

Abstract

本发明公开了一种用于消减天文观测数据中射频干扰信号的数据处理方法及装置，该方法不仅大幅度提高了射频干扰信号探测的可靠性和准确性，还能降低处理大量数据的时间消耗。该数据处理方法主要包括：获取待处理数据；统计数据中的数值大小，画出其柱状图；求出数据的经验累计分布函数；定义拟合函数并进行拟合；标记RFI；选定最佳滤波函数；对每一个数据点进行加权,用滤波函数对待处理的数据进行卷积。本发明运算效率高，速度快，能快速实时处理海量数据，该方法基于二维高斯函数对射频干扰的数据进行卷积处理，可利用二维高斯函数对数据进行多次迭代运算，使最后图像更加平滑，降低噪底，对强弱射频干扰信号均能有效消除。

Description

一种射频干扰的数据处理方法及装置

技术领域

本发明涉及无线电干扰消减领域，尤其针对射电天文数据中射频干扰消减的数据处理方法及装置。

背景技术

随着射电天文设备持续不断地更新发展，其观测性能越来越高。然而，由于射电天文分配的频谱资源非常有限，导致射电设备不断受到日益加强的非天文信号的影响，即射频干扰(Radio Frequency Interference，简称RFI)，其主要来源是人类的日常活动和通讯行为。射频干扰对观测结果的影响非常严重，大大降低了数据质量甚至产生数据无效，因此需要行之有效的消除射频干扰的技术手段，然而现有的射频干扰的数据处理技术具有如下缺陷：

(1)由于观测设备的更新升级，带宽不断增加，分辨率也越来越高，导致产生的数据量非常庞大，通常是TB量级的，而现有的RFI消除技术非常耗时，无法快速实时的处理海量数据。例如，一般，常规射频干扰的数据处理采用标记阈值方法，该方法是一种局域方法，将数据按时域或频域处理，处理的运算量是N(频率点数)*M(时间点数)*L(基线数量)，不适用于大数据；其相对应的标记算法，需先产生一个RFI统计表，然后在后续操作中删除，识别和删除数据是两个独立的操作过程，因此无法实时完成。

(2)弱RFI信号的处理识别效果差，检测RFI信号的真实度和灵敏度都不够，往往只能有效的识别强RFI信号，最终获得的数据结果与真实情况有很大的误差；例如，常规的标记算法在处理弱RFI信号方面效果不好，仅仅只能标识出RFI在数据序列中的位置，这样带来的潜在后果是强RFI信号会遗留一些残余在噪声中，而这些遗漏的噪声与弱RFI信号一起构成了噪底，给后续数据处理增加了难度，因此很难再进一步处理剩余的弱RFI信号。

(3)目前消减射频干扰的方法大多采用非线性的处理方式，无法进行迭代运算，这很大程度上限制了RFI噪声的消除精度。

因此，本领域的技术人员致力于开发一种射频干扰的数据处理方法及装置。

发明内容

有鉴于现有技术的上述不足，本发明所要解决的技术问题是提供一种可以快速实时的处理射电设备采集的大量图像数据，标记强RFI信号，同时可多次迭代处理以探测并消减弱RFI信号，获得高质量天文数据的射频干扰的数据处理方法及装置。

为了解决上述问题，本发明提供一种射频干扰的数据处理方法，包括以下步骤：

获取待处理数据及其滤波函数，所述待处理数据是射电设备获取的经过处理后的时频域数据，所述滤波函数基于核函数获得；

标记获取数据中的强RFI信号，包括求得数据的经验累计分布函数(empiricalcumulative distribution function,ECDF)，以及数据拟合

消减数据中的弱RFI信号，包括对所述待处理数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。

进一步地，获取数据的ECDF，将函数值大于0.99的部分标记为强RFI。

进一步地，将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图范围为数据的最大值到最小值，平均分配以得到组距；

对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；

基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据。

进一步地，还包括以下步骤：输出所述强干扰数据至文本。

进一步地，所述核函数采用二维高斯函数，其公式是

其中，x、y表示x*y矩阵的核大小，x、y分别为横纵方向上与核中心的距离，G(x,y)为在坐标点(x,y)处的概率密度，σ为标准差，e是比例系数，其中σ与经过卷积处理后的优化图像数据的平滑程度相关。

进一步地，采用多次所述卷积处理。

进一步地，所述滤波函数采用快速傅里叶变换后的二维高斯函数。

进一步地，所述待处理图像数据被网格化处理成多个待处理子图像数据，对每个待处理子图像数据并行处理以获取优化数据。

为了解决上述问题，本发明还提供一种射频干扰的数据处理装置，包括存储器、处理器及存储在存储器上可以在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序实现以下步骤：

获取待处理数据及其滤波函数，所述待处理数据是射电设备获取的经过处理后的时频域数据，所述滤波函数基于核函数获得；

标记获取数据中的强RFI信号，包括求得数据的经验累计分布函数(empiricalcumulative distribution function,ECDF)，以及数据拟合

消减数据中的弱RFI信号，包括对所述待处理数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。进一步地，所述处理器执行所述程序还实现以下步骤：

进一步地，获取数据的ECDF，将函数值大于0.99的部分标记为强RFI。

进一步地，将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图范围为数据的最大值到最小值，平均分配以得到组距；

对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；

基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据。

进一步地，还包括以下步骤：输出所述强干扰数据至文本。

进一步地，所述核函数采用二维高斯函数，其公式是

其中，x、y表示x*y矩阵的核大小，x、y分别为横纵方向上与核中心的距离，G(x,y)为在坐标点(x,y)处的概率密度，σ为标准差，e是比例系数，其中σ与经过卷积处理后的优化图像数据的平滑程度相关。

进一步地，采用多次所述卷积处理。

进一步地，所述滤波函数采用快速傅里叶变换后的二维高斯函数。

进一步地，所述待处理图像数据被网格化处理成多个待处理子图像数据，对每个待处理子图像数据并行处理以获取优化数据。

通过实施上述本发明提供的射频干扰的数据处理方法及装置，具有如下技术效果：

(1)运算效率高，速度快，能快速实时的处理海量数据，应用范围更广。

(2)基于二维高斯函数对射频干扰的数据进行卷积处理，由于二维高斯函数相当于对两个独立的一维空间分别进行计算，也就是说，首先将图像数据与一维高斯函数进行卷积，所述函数为：

式中x为离核中心的距离，G(x)为在离核中心x处的概率密度，σ为标准差。

然后，将卷积结果与方向垂直的相同一维高斯函数卷积，因此，二维高斯滤波的计算量是随着其中一维长度的增加而线性增长的，并不是平方的关系，因而可以利用二维高斯函数对数据进行多次迭代运算，使最后图像更加平滑，降低噪底，对弱射频干扰信号的消除更加明显。

(3)该射频干扰的数据处理方法可以并行化，进一步加快运算速度。例如针对大规模阵列数据，不同天线产生的数据可以分别并行处理再整合，也可以将数据分为几块分别进行处理。

(4)高精度的消减弱射频干扰数据。

以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。

附图说明

图1是本发明总体流程示意图；

图2是本发明优选实施例一的射频干扰的数据处理方法示意图；

图3是本发明优选实施例二的射频干扰的数据处理方法示意图；

图4是本发明优选实施例二的拟合结果示意图；

图5是本发明优选实施例二的最终RFI抑制结果图。

具体实施方式

下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面采用几个具体实施例对本发明上述实施例的技术方案进行详细说明。

实施例一

图2是本发明优选实施例一的射频干扰的数据处理方法示意图，如图1、图2所示，本实施例的射频干扰的数据处理方法可以包括：

步骤101、获取待处理数据及其滤波函数，所述待处理数据是射电设备获取的经过处理后的时频域数据，所述滤波函数基于核函数获得；

在本实施例中，引用了核函数对射电设备获取的时频域数据进行滤波处理，一方面由于现在射电设备的不断升级换代，采集的图像越来越清晰，分辨率愈来愈高，其获取的图像数据量越来越大，通过核函数的引入避免了“维数灾难”,大大减小了计算量。另一方面核函数可融合性能非常好，可与不同的算法相结合，形成多种不同的基于核函数的技术方案。

举例说明，本实施例中，根据射电设备采集的时频域数据，采用的核函数是二维高斯函数，其公式是

其中，x、y表示x*y矩阵的核大小，x、y分别为横纵方向上与核中心的距离，G(x,y)为在坐标点(x,y)处的概率密度，σ为标准差，e是比例系数，其中σ与经过卷积处理后的优化图像数据的平滑程度相关，也就是说σ决定滤波后图像的平滑程度。

步骤102、对所述待处理时频域数据进行卷积处理，包括对所述待处理时频域数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。换言之，将核函数与待处理时频域数据进行卷积，以获取待处理图像数据的加权平均数据结果即获得消减弱射频干扰信号的数据结果，具体是分别重新计算待处理图像数据中每个点的数值。

举例来说，在本实施例中具体实施为，利用快速傅里叶变换，将核函数与射电设备采集的图像数据进行卷积，重新计算图像中每点的数值，包括：移动核的中心，使其位于图像待处理点的位置，然后计算该中心的优化数据值，是由其本身与核范围大小的邻域内的其它数值经过加权平均后得到，其中所有权重由核函数决定。

需要说明的是，在另一种较佳实施方式中，可以对待处理图像数据进行多次步骤102所述的卷积处理，即利用核函数对数据进行多次迭代运算，使最后图像更加平滑，对弱射频干扰数据消除更加明显。例如，反复执行步骤102操作，直至图像效果满意为止。

需要说明的是，在另一种较佳实施方式中，可以对待处理时频域数据进行预处理以加快运算速度，即将所述待处理时频域数据被网格化成多个待处理子图像数据，对每个待处理子图像数据并行处理以获取优化数据即对每个待处理子图像数据并行执行步骤102操作。例如，具体实施方案有两种，一是针对不同天线产生的数据分别进行卷积处理，然后汇总数据；二是把所有数据分为多段进行并行化卷积处理。

综上所述，以上两个步骤是利用高斯滤波消除弱射频干扰数据的方法，使噪底下降，图像平滑。

实施例二

图3是本发明优选实施例二的射频干扰的数据处理方法示意图，其大体可以分为两部分：一、识别并标记强射频干扰数据；二、识别弱射频干扰数据并进行消减处理，该部分与实施例一的方法类似，如图3所示，本实施例的方法可以包括：

步骤201、获取待处理时频域数据，所述待处理图像数据是射电设备获取的时间序列数据转换至时频域的数据，

步骤202、获取数据的ECDF，将函数值大于0.99的部分标记为强RFI。

步骤203、将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图范围为数据的最大值到最小值，平均分配以得到组距；

步骤204、对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；

步骤205、基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据。

步骤206、输出所述强干扰数据至文本

步骤207、选取滤波函数，所述滤波函数基于核函数获得；

步骤208、对所述待处理数据进行卷积处理，包括对所述待处理图像数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。

步骤207、208的执行过程与实施例一所示的步骤101、102的执行过程类似，此处不再赘述。

在本实施例中采用仿真数据进行模拟。仿真法主要读取含有RFI噪声信号(即射频干扰数据)的模拟数据，然后通过高斯滤波和设定阈值的方法对RFI噪声信号进行处理，从而优化数据。在本实例中采用Python来实现仿真，Python中的astropy库比较全面，可以满足大部分射电设备采集的数据(例如天文数据)处理需求，可以很容易地对仿真数据模拟滤波、拟合过程，获取的模拟结果与真实情况接近。

下面以Python软件实现本实施例的技术方案为例，对本实施例的射频干扰的数据处理方法进行进一步详细的说明：

(1)对应于本实施中步骤201，首先创建仿真数据文件，文件为500*128矩阵的时频域数据，其中RFI噪声信号为随机分布加上固定的行、列。

(2)对应于本实施中步骤202，获取数据的ECDF，将函数值大于0.99的部分在原数据中标记为强RFI。标记结果为图4中深黑线条右侧的部分。

(3)对应于本实施中步骤203，将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图范围为数据的最大值到最小值，平均分配以得到组距；

(4)对应于本实施中步骤204，对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；拟合结果如图4

(5)对应于本实施中步骤205，基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据。

(6)对应于本实施中步骤206，输出经过步骤202和205标记的强干扰数据至文本。

(7)对应于本实施中步骤207，在Python中定义二维高斯函数

函数中具体数值的定义经过多次试验，优选为：标准差σ为0.5，核尺寸x*y为4*4的矩阵，比例系数e为10。

(8)对应于本实施中步骤208,利用快速傅里叶变换，将前面定义的二维高斯函数与仿真图像数据进行卷积，对整幅仿真图像进行加权平均，重新计算仿真图像中每点的数值。具体步骤为：移动核的中心，使其位于图像待处理点的位置，然后该中心的修改值(即优化数据)由其本身与核范围大小(即4*4的矩阵)邻域内的其它数值经过加权平均后得到，其中所有权重数值由前面定义的二维高斯函数决定；然后再次移动核的中心，重复该步骤，直至求的所有点的加权平均值。

需要注意的是，若有必要，可以重复运行此步骤，即迭代运算，直至图像满意为止。

最后结果见图5，左图为原始时频域数据所成的图，中间为经过本方案消减过的结果，从右边的强度条可以看出，整个图像信号强度由13降至9，图像更加平滑，右图中的白点为标记为强RFI的部分。

本实施例采用Python软件来实现该射频干扰的数据处理方法，先用滤波消除弱RFI噪声，再标记并删除强RFI信号，该方式在能够快速获得运行结果的同时又极大地提高了对噪声识别的灵敏度和结果信号的真实度，并且该算法具有实时性和全局性，而且能够实现并行化和迭代运算，大大改善了当前针对RFI噪声识别困难、效率低的处境。该算法实现方法的提出能够有效地促进天文、通讯及其它行业相关技术的快速发展。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种射频干扰的数据处理方法，其特征在于，包括以下步骤：

获取待处理数据及其滤波函数，所述待处理数据是射电设备获取的经过处理后的时频域数据，所述滤波函数基于核函数获得；

标记获取数据中的强RFI信号，通过获取数据的经验累计分布函数以及数据拟合两方面来判断；

消减数据中的弱RFI信号，包括对所述待处理数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。

2.如权利要求1所述的射频干扰的数据处理方法，其特征在于，获取数据的ECDF，将函数值大于0.99的部分标记为强RFI。

将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图横轴范围为数据的最小值到最大值，平均分配以得到组距；

对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；

基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据。

3.如权利要求2所述的射频干扰的数据处理方法，其特征在于，还包括以下步骤：综合统计两种方法中的强干扰数据，并输出所述强干扰数据至文本。

4.如权利要求1所述的射频干扰的数据处理方法，其特征在于，所述核函数采用二维高斯函数，其公式是：

<mrow> <mi>G</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <msup> <mi>&amp;pi;&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mfrac> <mrow> <msup> <mi>x</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> </mrow> </msup> </mrow>

其中，x、y表示x*y矩阵的核大小，x、y分别为横纵方向上与核中心的距离，G(x,y)为在坐标点(x,y)处的概率密度，σ为标准差，e是比例系数，其中σ与经过卷积处理后的优化图像数据的平滑程度相关。

5.如权利要求1所述的射频干扰的数据处理方法，其特征在于，采用多次所述卷积处理。

6.如权利要求1所述的射频干扰的数据处理方法，其特征在于，所述滤波函数采用快速傅里叶变换后的二维高斯函数。

7.如权利要求1所述的射频干扰的数据处理方法，其特征在于，所述待处理图像数据被网格化处理成多个待处理子图像数据，对每个待处理子图像数据并行处理以获取优化数据。

8.一种射频干扰的数据处理装置，其特征在于，包括存储器、处理器及存储在存储器上可以在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序实现以下步骤：

获取待处理数据及其滤波函数，所述待处理数据是射电设备获取的时间序列数据，所述滤波函数基于核函数获得；

标记获取数据中的强RFI信号，包括求得数据的经验累计分布函数，以及数据拟合；

消减数据中的弱RFI信号，包括对所述待处理数据中的每个待处理点数据进行加权平均，以获取各待处理点数据的优化数据，其中所述加权的权重基于所述滤波函数确定。

9.如权利要求8所述的射频干扰的数据处理装置，其特征在于，所述处理器执行所述程序还实现以下步骤：

获取数据的ECDF，将函数值大于0.99的部分标记为强RFI；

将获取的所有待处理点的优化数据基于组距处理获得柱状图数据，其中组距定义方法为：首先总组数定义为数据点数量总和除以10，柱状图范围为数据的最大值到最小值，平均分配以得到组距；

对所述柱状图数据进行最小二乘法拟合处理以获取拟合函数，考虑到自然界中噪声为高斯分布，因此所述拟合函数公式为y＝A*exp(-a*x)，其中A和a为待拟合值；

基于所得拟合函数，求出函数误差在6σ的比较函数，将柱状图中超过6σ部分的数据也标记为强干扰数据；

所述核函数采用二维高斯函数，其公式是：

<mrow> <mi>G</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <msup> <mi>&amp;pi;&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mfrac> <mrow> <msup> <mi>x</mi> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> </mrow> </msup> </mrow>

其中，x、y表示x*y矩阵的核大小，x、y分别为横纵方向上与核中心的距离，G(x,y)为在坐标点(x,y)处的概率密度，σ为标准差，e是比例系数，其中σ与经过卷积处理后的优化图像数据的平滑程度相关；

采用多次所述卷积处理；

所述滤波函数采用快速傅里叶变换后的二维高斯函数；

所述待处理图像数据被网格化处理成多个待处理子图像数据，对每个待处理子图像数据并行处理以获取优化数据。

10.如权利要求9所述的射频干扰的数据处理装置，其特征在于，所述处理器执行所述程序还实现以下步骤：输出所述强干扰数据至文本。