基于自适应滤波的干扰消除方法研究∗

马 琴 裴 鑫 陈卯蒸 托乎提努尔 王兆军 李 健段雪峰

(1新疆大学物理科学与技术学院 乌鲁木齐 830046)

(2中国科学院新疆天文台 乌鲁木齐 830011)

(3中国科学院大学 北京 100049)

(4新疆微波技术重点实验室 乌鲁木齐 830011)

1 引言

随着科技的发展,射电望远镜的灵敏度、观测视场以及空间分辨率等诸多性能都有着比较明显的提升.但与此同时,无线电业务领域的应用范围也在不断扩大,使得射电天文望远镜不可避免地受到非天文信号的影响,即射频干扰(Radio Frequency Interference,RFI),其主要来源包括无线通信、无线网络、雷达、广播电视、卫星和人们使用的各类电子设备[1].RFI不但使得观测天文信号质量降低,而且将耗费大量的信号处理资源和存储空间.因此如何避免或减少RFI被观测系统接收和记录,已经成为射电望远镜在观测中需要解决的突出问题.

传统的单通道射电天文干扰消除方法[2–5]在信噪比较低时只能去除部分干扰,且容易造成观测信号失真.基于自适应滤波器的干扰消除方法可以在信号和噪声的统计特性未知的前提下,从噪声中提取信号,当信号与噪声的统计特性发生变化时,能调节自身的参数以适应不同状况,从而实现最优滤波的特性[6].它通过将射电天文信号自身特点与雷达及无线电通信等领域微弱信号处理方法相结合,提出了自适应干扰对消的策略.该方法借助一架低增益参考天线来实时监测射频干扰,然后由自适应滤波器跟踪两个输入信号并去除RFI,保留所需的天文信号.

自适应干扰对消早在20世纪60年代就被引入,作为减少低频系统中干扰噪声的一种方法,成功应用于语音信号噪声消除、心电图周期噪声信号抵消等领域[7].其最主要的部件是自适应滤波器,而自适应滤波器在普通滤波器的基础上加入了调整滤波器系数的自适应算法,根据自适应算法优化准则可选取不同的算法,但这些算法的计算量差异较大,故而在选择算法时既要考虑性能也要考虑计算量的大小.自适应滤波器两种常用的算法是最小均方(Least Mean Square,LMS)和递推最小二乘法(Recursive Least Square,RLS)[8],RLS算法具有较快的收敛速度,但其运算量较大且具有发散性,LMS算法收敛速度一般,但收敛性好.本文将对这两种方法进行分析和比较,重点采用LMS进行建模,利用仿真软件对自适应滤波器进行设计和分析,以系统误差和收敛性为评判标准,通过改变步长与阶数对滤波效果进行优化.

2 研究内容

2.1 自适应干扰对消模型

自适应干扰对消系统是基于自适应滤波器原理的一种扩展[9],其原理如图1所示,它由主信道(射电望远镜)和一个单独的参考通道组成.主信道接收被RFI(无线电频率干扰)污染的天文信号,即进入射电望远镜的天文信号s(n)和干扰iP(n),其中iP(n)也包含了白噪声,n为时间序列.一架低增益参考天线被用来接收干扰iR(n),其中由于参考天线增益有限,因此在滤波器自适应时间尺度上的参考输入中基本上没有天文信号[10].

图1 自适应干扰对消原理图Fig.1 Schematic diagram of adaptive interference cancellation

参考信道iR(n)中的干扰与主信道iP(n)中的干扰相关,自适应滤波器工作的主要目的是将这种相关性作为时间的函数.自适应算法通过比较前一时刻和当前的信息,并发送更新后的加权系数到数字滤波器,数字滤波器利用这些系数改变iR(n),产生与主信道干扰相似的输出y(n);从一次输入中减去y(n),得到系统输出ε(n):

对于s(n)、iP(n)、iR(n)或它们之间的相互关系,不需要事先知道.自适应算法通过比较ε(n)和ε(n−1)找到新的系数,使用LMS算法最小化总功率:

由于s(n)与主路和参考路的干扰不相关,交叉项为零,因此系统输出的期望值(时间平均)为:

当滤波器对系数进行调整使E{ε2(n)}最小时,天文信号E{s2(n)}中的功率被吸收,因此E{[iP(n)−y(n)]2}达到功率最小,即:

由于天文信号是恒定的,使总输出功率最小就使输出的干扰功率最小,从而使输出信噪比最大.此时的ε(n)就是s(n)的最佳估计,y(n)实际上是iP(n)在最小均方意义下的估计.

2.2 自适应滤波器结构

自适应滤波器分为有限冲激响应(Finite Impulse Response,FIR)和无限冲激响应(Infinite Impulse Response,IIR)这两种结构[11].FIR滤波器由有限个脉冲响应函数离散值的滤波器组成,它在实际应用中最常用,可以实现较为严格的线性相移特性,并且可以保证滤波后波形不失真,但若要求频域过渡带快速衰减,就需要更多的阶数.IIR滤波器可以以较小的计算量获得陡降的过渡带,但其收敛慢、稳定性差,且较难保证线性相移,在FIR滤波器中则不存在类似问题.FIR滤波器结构又可以分为3种结构类型:格型结构、横向结构和对称横向结构,其中横向结构因其形式简单及易于实现最为常用,故本文采用FIR横向滤波器结构进行设计,滤波器的结构如图2所示.

滤波器的输出为:

其中XT(n)为滤波器输入矢量,即X(n)=[x(n)x(n−1)···x(n−N+1)]T,W(n)为权系数矢量(滤波器的冲激响应),即W(n)=[w0(n)w1(n)w2(n)w3(n)···w N−1(n)]T,N为滤波器的阶数[12].

图2 FIR滤波器结构Fig.2 The structure of FIR filter

2.3 自适应滤波器算法

本文采用LMS算法,因其具有计算量小、性能稳健、易于实现、不依赖于模型等优点,在实践中被广泛使用[13].

LMS算法是一种基于最小均方误差准则的随机梯度下降算法,核心思想就是利用均方误差代替平方误差,其通过调节权系数使得滤波器的均方误差最小,误差曲面的梯度为:

最陡下降法迭代计算权矢量的公式为:

式中µ是控制步长的参数称为自适应增益常数,是表征迭代速度快慢的物理量.W(n+1)和W(n)分别为迭代后和迭代前的系数值,从上式可以看出,自适应迭代下一时刻权系数矢量由上一时刻的权系数矢量加以误差函数为收敛参数的输入矢量得到.自适应滤波器收敛的条件是:

其中λmax是输入信号自相关矩阵的最大特征值;当µ值越大自适应过程越快且时间越短,引起的失调越大,µ值越小,自适应过程越慢,系统越稳定同时失调也越小,µ的选取必须在失调和收敛速度之间取得较好的折中,既要具有较快的收敛速度,又要使稳态误差最小.

算法步骤:

(1)初始化滤波器系数,W(0)=0,或者可以根据先验知识来确定初始权值.

(2)对时间序列每一时刻进行以下计算:

滤波:y(n)=XT(n)W(n),误差:ε(n)=s(n)+iP(n)−iR(n),权向量更新:W(n+1)=W(n)+2µε(n)X(n).

综上所述,该方法的主要步骤为:(1)读入主信道和参考信道的信号;(2)设置滤波器的阶数和步长;(3)初始化自适应滤波算法的参数;(4)自适应滤波处理;(5)滤波器权系数的更新.

3 仿真与分析

基于上述理论分析,由(5)式和(7)式可知,权系数矢量和滤波器输入矢量乘积的和即为滤波器的估计输出,加权系数矢量是通过步长、ε(n)以及与上一次迭代权值之和计算得到.根据图1本文设计了不同阶数和不同步长的自适应滤波器,对加干扰的正弦波进行滤波,输入信号如图3所示.

图中,主通道中输入目标信号、干扰以及白噪声,参考信道输入的只有干扰和白噪声,通过LMS算法可以自适应调节线性组合器权系数,主信道与参考信道内的噪声信号相互对消,输出的信号即为所需的目标信号.

图3 输入信号Fig.3 The input signal

为获取滤波器最优参数,本文首先分析了当阶数一定时改变步长对滤波效果的影响,选取了阶数为2、步长设置分别为0.0005、0.0002、0.00009的滤波器进行仿真;然后分析了当步长一定时改变阶数对滤波效果的影响,选取了步长为0.0001、阶数设置分别为2、6、12的滤波器进行仿真;同时用误差对滤波性能进行了表征.如图4和图5分别是2阶时不同步长下滤波后的信号和误差,图6和图7分别是步长为0.0001时不同阶数的滤波信号和误差.

当阶数一定,改变步长时,由图4可知,计算结果随步长增加收敛速度变快,步长越短的计算时间越长.滤波性能随参数的改变在图5中可以看到,误差的收敛速度与滤波信号的收敛速度同步,但滤波效果随步长增大而变差.

图5 2阶时不同步长的误差(上:0.0005、中:0.0002、下:0.00009)Fig.5 The error of different step sizes when the order is 2(top panel:0.0005,middle panel:0.0002,bottom panel:0.00009)

图5 续Fig.5 Continued

图6 步长为0.0001时不同阶数的滤波信号(上:12、中:6、下:2)Fig.6 Filtered signals of different orders when the step size is 0.0001(top panel:12,middle panel:6,bottom panel:2)

图7 步长为0.0001时不同阶数的误差(上:12、中:6、下:2)Fig.7 The error of different orders when the step size is 0.0001(top panel:12,middle panel:6,bottom panel:2)

当步长一定改变阶数时,如图6所示,滤波信号收敛时间随着阶数的减小而增大,滤波信号的收敛速度随着阶数的增大而变快.滤波效果如图7所示,随着阶数的减小,误差也更接近0.

在进行多次仿真实验后发现,当阶数不同时,取不同的步长也可以达到较好的滤波效果.当阶数增大时,为了得到较好的滤波效果,步长也将随之增大.通过仿真我们发现,自适应滤波器在保证快速滤波的前提下,能有效还原信号的轮廓,将误差降到最小.

4 观测数据验证

为了检验该算法的有效性,本文使用新疆天文台南山观测基地观测的脉冲星数据进行了滤波实验.由于南山26 m射电望远镜是单天线观测,无法获取实验中所需的实时参考信号,本文将双极化信号中的一路作为参考信号,另一路作为主望远镜信号,将两路信号进行对消,验证自适应滤波器的有效性.实验采用L波段接收机观测的脉冲星J0332+5434数据,数据记录为天文VDIF(VLBI(Very Long Baseline Interferometry)Data Interchange Format)格式,信号带宽为256 MHz,中心频率为1428 MHz.数据处理时首先将记录的两个文 件snap2input02020-11-0710-46-36(极 化1)和snap2input12020-11-0710-46-34(极化2)分别作为主通道和参考通道的信号输入自适应滤波器,选择优化后的滤波参数对其进行滤波,阶数为2,步长为0.00015.然后利用VLBI处理软件DIFX(Distributed FX Correlator)对数据进行折叠,处理结果如图8所示.从上至下,第1、2幅图分别为未经滤波的极化1和极化2脉冲星轮廓,第3幅为采用自适应滤波器处理后的轮廓,具有相关性的两个轮廓基本被消除.由此可推断,如将该方法应用至实际观测中,参考通道与主通道中具有相关轮廓的信号均会被消除,即干扰信号被消除,因为干扰信号将同时出现在参考通道和主通道中.而天文信号将被保留,因为微弱的天文信号仅存在于主通道中.

本文算法还使用Parkes 64 m射电望远镜超宽带接收机经系统RFI处理后的测试数据进行了算法性能验证.射电望远镜观测脉冲星J0332+5434,3 m口径的小天线指向固定干扰源(4G塔),同时记录数据.观测系统将主信道和参考信道的数据交替记录在一个PSRDada格式的文件中,主信道和参考信道的数据都是16比特复数信号(实部和虚部各16比特),带宽为128 MHz.

首先,利用天文信号处理软件(DSPSR、PSRCHIVE等)进行消色散及折叠,判断记录的数据中是否存在脉冲星信号.折叠之后,在主信道中可以看到脉冲星J0332+5434的清晰轮廓,由于参考天线增益低、口径小、灵敏度低,在实际观测中无法收到微弱的天文信号,因此在参考信道中看不到脉冲星信号轮廓.然后,使用自适应滤波器进行RFI消除测试,数据处理时先将.dada文件中的主信道和参考信道分离出来,再将信号输入至本模型进行滤波处理,参数设置为阶数6、步长0.00006.图9显示了滤波前后的信号频域图,从上至下依次为:主天线信号(Parkes 64 m)、参考天线信号和滤波后的信号.从图中可以看出,位于0点位置的强干扰信号在经过本滤波模型后得到了明显的削弱.

5 讨论与总结

天文信号是一种极其微弱的天体辐射,为了观测到清晰可见的天文信号需要使用高灵敏度的观测设备.天线接收面积对观测灵敏度的影响很大,因此一般使用大口径射电望远镜进行天文观测.

图8 滤波前后的轮廓图Fig.8 The profile figure before and after filtered

图9 滤波前后的频域图Fig.9 Frequency domain figure before and after filtered

RFI主要来自于人类的通信活动和日常生活,例如无线电广播、电视信号、3G和4G通信信号等,这些干扰信号都很强,通过小天线接收.

主天线和参考天线都接收到天文信号的概率较低,主要原因如下:

(1)主天线与参考天线的指向不同.在观测过程中,射电望远镜跟踪天文信号源,而参考天线固定指向RFI源,二者指向重合的概率比较低;

(2)参考天线一般采用低增益小口径天线,灵敏度低,而绝大多数天文信号都非常微弱,能接收到天文信号的概率非常低.

以Parkes 64 m射电望远镜为例,在测量中采用了3 m口径的参考天线和常温接收机,信号带宽为128 MHz,系统噪声温度为120 K.对于目前天文上探测到的信号强度极高的快速射电暴(FRB)来说,参考天线的探测灵敏度为:

其中,Smin是在给定阈值γ0下的可探测流量密度,β是数字化因子,Tsys是系统温度,G是望远镜增益,BW是带宽,NP是偏振数[14].当天线效率取60%时,GDPFU(Degree Per Flux Unit)计算后为0.0015 K·Jy−1,τ为FRB信号的爆发持续时间,取3 ms计算,探测阈值取7时的最小可探测流量密度Smin约为667 Jy,通过查询FRB目录网站1https://frbcat.org/,目前探测到的FRB很少有超过这个数值的.此外,FRB信号出现在特定天区的特定时刻,而参考天线固定指向某一特定方向,因此,参考天线接收到天文信号的概率非常低.

本文分析了自适应滤波器的结构和算法,根据自适应干扰对消原理,设计并建立了适用于射电天文的自适应滤波器仿真模型,通过大量仿真实验得出了算法参数与滤波器性能之间的关系,获取了较优的参数模型,并将南山26 m射电望远镜和Parkes 64 m射电望远镜的观测数据输入至该模型,验证了该自适应滤波器的有效性.自适应滤波器RFI消除方法能滤除干扰信号而保留天文信号,是其他的RFI消除方法所不具备的.由于缺乏参考天线以及相关信号采集设备,本文仅对该方法进行了仿真和部分验证,作为启发性工作,未来将持续开展相关研究,以期将该方法应用至南山26 m射电望远镜、即将建设的新疆110 m大口径全可动射电望远镜以及国内其他大型射电望远镜.