一种射电天文信号的自适应射频干扰消除方法与流程

本发明属于射电天文领域，尤其涉及一种射电天文信号的自适应射频干扰消除方法，其适用于脉冲星、快速射电暴等天文信号的干扰处理及相关的科学研究。

背景技术：

随着信息技术与制造技术的快速发展，射电天文装备不断得到升级或换代，以使其具有更好的观测性能，包括提高数据记录的时间和频率分辨率以及获得更高的接收和记录带宽等。然而与之形成矛盾的是：国际电信联盟仅为射电天文分配了非常有限的频谱资源，导致的后果是射电观测设备不可避免地受到日益增强的非天文信号的影响，这些干扰信号主要来自于人类的通信活动和日常生活，这些人为的电磁信号即是射电天文中所说的射频干扰(radiofrequencyinterference，rfi)。

rfi(radiofrequencyinterference，射频干扰)作为观测目标之外的信号，是影响天文观测设备灵敏度的主要因素之一。射频干扰会降低数据质量甚至导致数据无效，将严重影响科学研究结果的准确性与有效性。

来自宇宙的天体辐射及其微弱，因此要求射电望远镜具备很高的灵敏度。由于射电望远镜灵敏度高、信号收集能力强，在接收到有用天文信息的同时人为的干扰也被接收。射频干扰作为观测目标之外的信号，是影响天文观测设备灵敏度的主要因素之一。如何快速准确地识别并抑制射频干扰已成为射电观测过程中急需解决的问题。

目前，世界各地的大口径射电望远镜均严重受到射频干扰的影响，天文观测面临商业、国防及其他利益方面使用的无线电频谱带来的巨大压力。射电天文学家不能再依靠无线电管理机构来实现无干扰的观测环境，需要研究和探索共存的可能性。现在，我国的fast500米、tianma65米、nsrt25米等望远镜在低频段天文观测中射频干扰比较严重，尚未找到良好的解决办法。

射电望远镜虽然建设在比较偏远、无线电干扰相对较小的区域，但是其天线接收面积大、天线增益高，灵敏度远远超过普通的通信设备，因此容易受到其它信号的干扰，包括电视信号、调频无线电传输、gps、手机和飞机导航通讯等。这些强干扰导致重要频段天文数据的丢失，大大降低天文信号的信噪比，对脉冲星、frb等天文现象研究产生很大的障碍。面对rfi环境越来越复杂的趋势，我们需要进一步研究新的射频干扰(rfi)消除算法和消干扰技术，以提高天文观测设备灵的敏度和数据质量。

目前使用的很多射频干扰(rfi)处理方法是在时域或频域里设立不同的阈值，如果信号强度大于阈值，则它被认为是射频干扰并剔除，rfi的识别和消除相关算法无法恢复原数据。

国外许多专家和学者对射电天文rfi处理理论和技术进行了深入研究，并实现了很多rfi抑制算法。一些著名天文研究机构也逐步关注并采取一系列抗干扰措施减少干扰对射电天文观测的影响，例如：美国国家天文台(nrao)、澳大利亚国家望远镜机构(atnf)、荷兰天文学研究机构(nfra)、柏克莱大学(byu)等天文研究机构均提出了rfi应对策略。虽然国内国家天文台、上海天文台等天文研究机构，围绕fast、天籁计划等项目做了相关的研究，但针对有效rfi消除技术的研究工作还不够深入。现有国内外常用的rfi消除方法包括阈值法、空间滤波、中值滤波、kurtosis滤波、mad滤波等，其缺点在于天文信号受到一定的损失，只得到比较接近的值，无法恢复原数据，如果干扰信号很强且带宽比较宽，该频段数据无法使用。

观测数据的好坏关系到科学成果的质量甚至结论的真伪,射电频率干扰抑制和消减技术也在科学需求的指引下得以不断进步和发展。但是，越来越复杂的rfi环境对消干扰技术提出了很大的挑战，无法使用之前使用的简单消干扰技术，需要研究新的rfi消除算法，进一步完善消干扰策略。

技术实现要素：

本发明的目的在于提供一种射电天文信号的自适应射频干扰消除方法，以解决观测设备灵敏度和信号信噪比下降问题，提升射电天文系统的观测及搜寻的工作效率。

为了实现上述目的，本发明提供一种射电天文信号的自适应射频干扰消除方法，包括：

步骤s1：获取参考数据和观测数据，所述观测数据是通过射电望远镜观测并记录射电天文信号获得的，所述参考数据是由一参考天线系统观测并记录射频干扰参考信号而获得的；

步骤s2：对参考数据和观测数据进行通道划分，生成多个子带通道中的参考数据和观测数据的数据块；

步骤s3：初始化一自适应滤波器的增益函数g(t)和该自适应滤波器的输出；

步骤s4：计算各个所述参考数据和观测数据的每个数据块的平均功率；

步骤s5：根据所述平均功率计算标准化参数；

步骤s6：使用标准化参数计算自适应滤波器的输出和参考数据的互相关函数；

所述互相关函数c(i)为：

式中，vast(i)-g(n)×vref(i)为自适应滤波器的输出，vref为参考数据，f(n)为标准化参数；

步骤s7：对各个互相关函数进行积分并取平均值，以获取互相关平均值；

步骤s8：采用所述互相关平均值计算所述自适应滤波器的增益函数；

步骤s9：使用所述自适应滤波器的增益函数来修改参考数据的幅值和相位，并获取自适应滤波器的输出；

所述自适应滤波器的输出为：

vfilt＝va-gvr，

式中，vfilt为所述自适应滤波器的输出，va为观测数据，g为自适应滤波器的增益函数，vr为参考数据；

步骤s10：循环执行所述步骤s4-s9，直到数据处理完成为止。

在所述步骤s1中，所述射电望远镜系统和参考天线系统的记录同时进行，或者通过硬件或软件方式进行数据同步。

在所述步骤s2中，所述对基带数据进行通道划分是通过采用多相滤波器组进行的，所述多相滤波器组包括多路fir滤波器和与其相连的离散傅里叶变换装置，所述多路fir滤波器与一接合器相连。

在所述步骤s2中，每个子带通道的带宽为0～1mhz。

在所述步骤s3中，所述自适应滤波器的增益函数的初始值为0。所述自适应滤波器的输出的初始值为所述观测数据。

在所述步骤s4中，所述计算为并行计算，

且所述参考数据和观测数据的每个数据块的平均功率为：

其中，pref(n)，past(n)分别为参考数据和观测数据的每个数据块的平均功率，n为数据块索引，i为采样索引，k为数据块长度，vref，vast分别为参考数据和观测数据。

在所述步骤s5中，所述标准化参数f(n)为：

f(n)＝past(n)+pref(n)，

其中，u是更新因子，pref(n)，past(n)分别为参考数据和观测数据的每个数据块的平均功率，n为数据块索引。

在所述步骤s7中，所述互相关平均值为：

其中，c为互相关函数，cav为互相关平均值，n为数据块索引，i为采样索引，k为数据块长度。

在所述步骤s8中，所述自适应滤波器的增益函数为：

g(n)＝g(n-1)+ε×cav(n)，

式中，g滤波器增益，n是整数，n＝1,2,3,….,cav为互相关平均值，ε为自适应滤波器的固定常数。

所述自适应滤波器的固定常数ε为小于1的复数。

本发明的射电天文信号的自适应射频干扰消除方法按照参考数据与观测数据之间的相关性特性设计互相关器，根据参考数据灵活实现了自适应滤波器的增益函数的自动调整，使参考数据更接近观测数据的干扰成分并从观测信号中减去，有效消除了混入天文信号中的射频干扰，具有很高的抗干扰性能，从而实现观测数据中射电天文信号的无损恢复，解决观测设备灵敏度和信号信噪比下降问题，提升脉冲星观测及搜寻的工作效率，提高了观测系统灵敏度。此外，本发明的射电天文信号的自适应射频干扰消除方法通过计算标准化参数，解决了自适应滤波器的增益函数对输入信号强度变化的依赖性问题，使其不会受到滤波器输入信号幅值和相位变化的影响，算法性能更加稳定。此外，本发明采用多相滤波器组方法，首先把基带信号分成多路子带信号，然后在每一路信号进行独立的rfi消除，使得算法数据处理速度很大提高，还控制信号分通道过程中产生的频谱泄露，在子带通道上有效消除rfi，易于gpu和fpga并行实现，rfi处理速度快，算法稳定性高。

附图说明

图1为根据本发明的一个实施例的射电天文信号的自适应射频干扰消除方法的流程图。

图2为如图1所示的射电天文信号的自适应射频干扰消除方法的所采用的多相滤波器组的实现结构图。

图3为如图1所示的射电天文信号的自适应射频干扰消除方法的所采用的fir滤波器的结构示意图。

图4(a)-图4(b)为如图1所示的射电天文信号的自适应射频干扰消除方法的自适应rfi消除结果图，其中图4(a)示出了滤波之前的观测数据，图4(b)示出了滤波之后的观测数据。

图5为未采用本发明的射电天文信号的自适应射频干扰消除方法进行滤波的脉冲星轮廓图。

图6为采用本发明的射电天文信号的自适应射频干扰消除方法滤波之后的脉冲星轮廓图。

图7为本发明的射电天文信号的自适应射频干扰消除方法的自适应滤波器的增益函数的演化示意图。

具体实施方式

如图1所示为根据本发明的一个实施例的射电天文信号的自适应射频干扰消除方法的流程图。本发明的自适应射频干扰消除方法可适用于各种射电天文信号，在本实施例中，所述射电天文信号为脉冲星信号。

如图1所示，所述射电天文信号的自适应射频干扰消除方法具体包括以下步骤：

步骤s1：获取基带数据，该基带数据包括参考数据和观测数据。

其中，所述基带数据采用自适应消干扰系统来获取，所述自适应消干扰系统包括参考天线系统(即，参考系统)和射电望远镜(即，观测系统)，在实际应用中，所述观测数据是通过射电望远镜观测并记录射电天文信号获得的。参考数据是由参考天线系统观测并记录射频干扰参考信号而获得的。射电望远镜系统和参考天线系统的记录需要同时进行，如果存在延迟，可以通过硬件或软件方式进行数据同步。

观测数据为：

va(t)＝va_rx(t)+va_rfi(t)+va_ast(t)(1)

式中，va_rx(t)是射电望远镜的噪声，va_rfi(t)、va_ast(t)分别为射电望远镜接收的射频干扰及射电天文信号。

参考数据为：

vr(t)＝vr_rx(t)+vr_rfi(t)+vr_ast(t)(2)

式中，vr_rx(t)是参考天线系统的噪声，vr_rfi(t)、vr_ast(t)分别为参考天线系统接收的射频干扰及射电天文信号。

步骤s2：对参考数据和观测数据进行通道划分，生成多个子带通道中的参考数据和观测数据的数据块。

其中，每个子带通道的带宽越细越好，但是需要根据实际应用需求而设定，因此，每个子带通道的带宽为0～1mhz比较合适。所述对基带数据进行通道划分是通过采用多相滤波器组(polyphasefilterbank，pfb)进行的，信号处理过程包含混频、抽取、低通滤波、fft变换等一系列信号处理流程。

所述多相滤波器组的算法实现结构如图2所示，所述多相滤波器组包括多路用于多相分解的fir(有限长单位冲激响应)滤波器和与其相连的dft(离散傅里叶变换)装置，fir滤波器的具体结构如图2中虚线框以及图3所示，其由混频、抽取、低通滤波的信号处理流程组成，分别用于对其中一路信号进行子滤波，获取其子带通道中的参考数据和观测数据的数据块。

此外，所述多路fir滤波器与接合器相连，并通过所述接合器分配每一路fir滤波器的输入数据的采样，实现多相分解。因此在所述多相滤波器组(pfb)的算法设计中，首先设计基于窗函数的fir滤波器，然后对fir滤波器的系数进行多相分解。其中，对fir滤波器的系数进行多相分解，是指将fir滤波器的传输函数h(z)分解,并形成不同相位的m个分支，然后在每个分支上实现数字滤波。所形成的多相滤波结构可以降低原本阶数(即，tap数)较高的滤波器，有效提高fir滤波器速度，为高速率实时信号处理提供了实现途径。

所述对fir滤波器的系数进行多相分解，具体参见以下公式：

如图3所示，fir滤波器是信号处理系统中最基本的元件，是一种非递归型滤波器。fir滤波器可以做到严格的线性相移，系统性能始终非常稳定，其单位冲击响应有限长的，可以设计满足特殊频率特性形状要求的滤波器。

其中，fir滤波器的计算公式为：

式中，h(k)为滤波器的固定系数，x(n-k)为经过k个延迟单元的输入数据，n是冲击函数长度，n＝1,2，…，n。

由此，fir滤波器的传输函数h(z)可表示成：

式中，h(n)为滤波器的固定系数，n是冲击函数长度，n＝1,2，…，n，z表示z变换，是信号处理中最常用的概念。

fir滤波器的传输函数h(z)分解所形成的不同相位的m个分支为：

如果设，

k＝0，1，2，….，m-1

则，该不同相位的m个分支表示为：

其中，pk(z^m)z^-k是fir滤波器的传输函数h(z)的多相表示。

将z变换可转换成傅里叶变换，则：

式中，k＝0，1，2，…，m-1，因此分解的fir滤波器具有不同的相位。

由此，形成了m阶的fir滤波器，每一路采样信号和多相分解的fir滤波器系数卷积运算，最终用于进行离散傅里叶变换。其中，fir滤波器的阶数m根据实际应用设定，目前常用的fir滤波器包括4，8，16，32，64或128阶的fir滤波器，即m通常为4，8，16，32，64或128。

一般记录的信号带宽比较宽，为了有效地消除干扰，需要进行通道划分。在现有技术中，对基带数据进行通道划分一般通过采用脉冲星观测设备并使用快速傅里叶变换(fastfouriertransform,fft)进行，从而生成多个子带通道的数据作为观测数据文件。但是，fft产生的频谱泄露会降低脉冲星终端系统的观测分辨率，对微弱的射电信号研究影响较大。由此，本发明的射电天文信号的自适应射频干扰消除方法通过采用多相滤波器组(polyphasefilterbank，pfb)进行通道划分，在基于窗函数的有限长单位冲激响应(finiteimpulseresponse，fir)滤波器基础上使用数字滤波技术，从而能够有效降低信号的频谱泄露。

步骤s3：初始化自适应滤波器的增益函数g(t)和该自适应滤波器的输出。

其中，自适应滤波器是整个消干扰系统的核心部分，它修改来自参考数据的幅值和相位，并从天文信号中减去干扰信号，从而获得干净的天文数据。该方法的最大优势是它不会影响天文数据，只消除天文信号的rfi成分，不会导致数据的丢失，从而有效提高观测系统信噪比，因此，自适应滤波是一种可以优先选择的消干扰方案。自适应消干扰系统为天文学家直接进入脉冲星消色散，折叠和脉冲星计时计算等数据后处理阶段提供高质量的观测数据。

自适应滤波器的增益函数g(t)用于在下文中修改所述参考数据的幅值和相位。自适应滤波器的增益函数g(t)的初始值可以设置为0。所述自适应滤波器的输出的初始值为所述观测数据。

步骤s4：计算各个所述参考数据和观测数据的每个数据块的平均功率。该平均功率用于自适应滤波器的增益函数g(t)的更新。其中所述步骤s4中的计算为并行计算，使得每一个数据块的数据处理相互独立，以实现并行计算。

其中，所述参考数据和观测数据的每个数据块的平均功率为：

其中，pref(n)，past(n)分别为参考数据和观测数据的每个数据块的平均功率，n为数据块索引，i为采样索引，k为数据块长度，vref，vast分别为参考数据和观测数据。数据块的大小根据实际应用需求设置。

步骤s5：根据所述参考数据和观测数据的每个数据块的平均功率计算标准化参数。所述标准化参数用于解决自适应滤波器的增益函数g(t)对输入信号功率强度变化的依赖性问题。

所述标准化参数f(n)为：

f(n)＝past(n)+pref(n)(7)

其中，u是更新因子，可根据实际应用需求设计u的大小，pref(n)，past(n)分别为参考数据和观测数据的每个数据块的平均功率，n为数据块索引。

步骤s6：使用标准化参数计算自适应滤波器的输出和参考数据的互相关函数c(i)。其中，自适应滤波器的输出和参考数据之间的相关性越大，则说明算法rfi消除性能越好。

所述互相关函数c(i)为：

式中，vast(i)-g(n)×vref(i)为自适应滤波器的输出，vref，vast分别为参考数据和观测数据，f(n)为标准化参数，i为采样索引，n为数据块索引，g(n)为自适应滤波器的增益函数。

由此，构成了一个反馈回路，通过自适应滤波器的增益函数g(t)直到相关性为零，使得互相关函数c(i)＝0的时候，输出信号中的rfi完全被消除，得到干净的脉冲星信号。

步骤s7：对各个互相关函数进行积分并取平均值，以获取互相关平均值。由此，使用积分方法，提高算法性能稳定性，尽量保持自适应滤波器的增益函数的线性状态。

所述互相关平均值为：

其中，c为互相关函数，cav为互相关平均值，n为数据块索引，i为采样索引，k为数据块长度。

步骤s8：采用所述互相关平均值计算所述自适应滤波器的增益函数。其中，自适应滤波器的增益函数直观展示了自适应滤波器的消干扰性能，随着rfi信号的变化，滤波器增益在稳定点的附近波动变化。

所述自适应滤波器的增益函数为：

g(n)＝g(n-1)+ε×cav(n)(11)

式中，g滤波器增益，n是整数，n＝n＝1,2,3,….,cav为互相关平均值，ε为自适应滤波器的固定常数。

自适应滤波器的ε值是一个固定的常数，它确定滤波器响应达到稳定状态的快慢，直接影响rfi消除算法的性能。如果ε取值过小，自适应滤波器需要更长的时间才能达到稳定状态；如果取的值太大，滤波器的消干扰性能失效，导致输出数据的迅速增大，得不到正确的输出结果。一般ε值取小于1的复数，且大约等于1/rfi。

步骤s9：使用所述自适应滤波器的增益函数来修改参考数据的幅值和相位，并获取自适应滤波器的输出。这里，经过上文的对输入信号进行的一系列处理之后，这里计算得到干净的射电天文信号，如脉冲星信号。

其中，修改后的参考数据为g×vr，由此，通过复数乘法运算和自适应滤波器的增益函数来修改参考数据的幅值和相位。

所述自适应滤波器的输出为：

vfilt＝va-gvr

＝va_rx+(va_rfi-gvr_rfi)+(va_ast-gvr_ast)-gvr_rx(12)

式中，vfilt为所述自适应滤波器的输出，va为观测数据，g为自适应滤波器的增益函数，vr为参考数据，va_rx为观测系统的噪声，vr_rx为参考系统的噪声，va_rfi为观测天线收到的rfi，vr_rfi为参考天线收到的rfi，va_ast为观测天线收到的脉冲星信号，vr_ast为参考天线收到的脉冲星信号。

步骤s10：循环执行所述步骤s4-s9，直到数据处理完成为止。

因为有多个数据块，不能一次性处理所有数据块，由此，本发明通过对每个数据块处理一次，如此循环处理，并且每个数据块的增益参数传输到下一个数据块，以通过不断的循环执行来降低自适应滤波器的增益函数的误差。

实验结果

模拟仿真两组数据，其中，天文信号(即观测数据)有3个信号分量：白噪声、rfi和脉冲星信号。参考数据只包含白噪声及干扰信号，参考数据中忽略脉冲星信号成分，因为参考天线口径小、灵敏度低，并且参考天线指向定到rfi源，因此无法收到极其微弱的脉冲星信号。观测数据和参考数据的噪声相对独立、不相关，并且rfi幅值和相位是随机的。

自适应rfi消除结果如图4(a)-图4(b)所示，图4(a)是使用自适应滤波器之前的脉冲星信号时域图，图4(b)是滤波后的脉冲星信号。

从图4(b)中可以看出，rfi滤波器消除了所有干扰信号，脉冲星信号通过滤波器之后能够清楚看见4个脉冲星信号。滤波器输出只包含脉冲星信号和白噪声，大大提升了天文信号信噪比。

采用本发明的射电天文信号的自适应射频干扰消除方法进行滤波前后的脉冲星号轮廓对比如图5、图6所示。

图5、图6分别表示了未采用本发明的射电天文信号的自适应射频干扰消除方法滤波和采用本发明的射电天文信号的自适应射频干扰消除方法滤波的脉冲轮廓图。图5、6中可以看到，采用射电天文信号的自适应射频干扰消除方法进行rfi消除之前，脉冲轮廓很大失真，无法判断脉冲星实际轮廓。消干扰之后的脉冲轮廓非常接近实际轮廓，有效消除rfi信号，信噪比很大提高，脉冲轮廓得到了很大改善。

实验结果显示，本发明具有非常好的消干扰性能，并且滤波器响应达到稳定状态不依赖于输入信号幅值的变化。

图7是自适应滤波器的增益函数的演化示意图，决定自适应rfi消除算法的性能。图7中可以看到，数据块索引值在30之后自适应滤波器达到稳定状态，滤波器稳定状态之后才开始正常工作。自适应滤波器增益越稳定消干扰算法的rfi消除性能越好。消干扰算法根据rfi的变化，不停地自动调整滤波器增益，保持滤波器的稳定状态。

以上所述的，仅为本发明的较佳实施例，并非用以限定本发明的范围，本发明的上述实施例还可以做出各种变化。凡是依据本发明申请的权利要求书及说明书内容所作的简单、等效变化与修饰，皆落入本发明专利的权利要求保护范围。本发明未详尽描述的均为常规技术内容。