CN117335772A - 一种天文瞬时干扰抑制方法 - Google Patents
一种天文瞬时干扰抑制方法 Download PDFInfo
- Publication number
- CN117335772A CN117335772A CN202311106072.XA CN202311106072A CN117335772A CN 117335772 A CN117335772 A CN 117335772A CN 202311106072 A CN202311106072 A CN 202311106072A CN 117335772 A CN117335772 A CN 117335772A
- Authority
- CN
- China
- Prior art keywords
- signal
- threshold value
- astronomical
- noise signal
- value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 230000001629 suppression Effects 0.000 title claims abstract description 23
- 230000003044 adaptive effect Effects 0.000 claims abstract description 93
- 238000001914 filtration Methods 0.000 claims abstract description 16
- 238000000926 separation method Methods 0.000 claims abstract description 13
- 238000005070 sampling Methods 0.000 claims description 31
- 238000000605 extraction Methods 0.000 claims description 5
- 230000014509 gene expression Effects 0.000 claims description 3
- 230000001052 transient effect Effects 0.000 claims description 3
- 230000003111 delayed effect Effects 0.000 claims description 2
- 238000004422 calculation algorithm Methods 0.000 abstract description 15
- 238000005516 engineering process Methods 0.000 abstract description 6
- 238000012545 processing Methods 0.000 abstract description 4
- 230000006870 function Effects 0.000 description 55
- 238000004364 calculation method Methods 0.000 description 15
- 230000001360 synchronised effect Effects 0.000 description 7
- 238000005311 autocorrelation function Methods 0.000 description 5
- 238000005314 correlation function Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 239000011159 matrix material Substances 0.000 description 5
- 230000000875 corresponding effect Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000006978 adaptation Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000007781 pre-processing Methods 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 230000001627 detrimental effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000007667 floating Methods 0.000 description 1
- 229910052739 hydrogen Inorganic materials 0.000 description 1
- 239000001257 hydrogen Substances 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007935 neutral effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
Classifications
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
- H03H21/00—Adaptive networks
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
- H03H21/00—Adaptive networks
- H03H21/0012—Digital adaptive filters
- H03H21/0043—Adaptive algorithms
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04B—TRANSMISSION
- H04B15/00—Suppression or limitation of noise or interference
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Signal Processing (AREA)
- Filters That Use Time-Delay Elements (AREA)
Abstract
本发明提供一种天文瞬时干扰抑制方法包括：将望远镜观测信号通过阈值进行信号分离，提取出噪声信号；将通过阈值提取出来的噪声信号经过LMS自适应滤波器滤波，得到待去除的噪声信号的估计值，最终将望远镜观测信号减去待去除的噪声信号的估计值，得到误差函数，将误差函数作为经过干扰抑制的天文信号。本发明的方法基于RFI粗分离技术和LMS自适应滤波器，无需使用参考天线数据就可以进行LMS自适应噪声抵消操作，粗提取的RFI信号完全可以满足算法要求；改进后的RFI处理流程能够直接使用射电望远镜输入信号进行LMS自适应滤波操作，而不需参考天线的信号输入。
Description
技术领域
本发明属于射电天文领域，具体涉及一种天文瞬时干扰抑制方法，其用于在没有参考天线的情况下使用最小均方误差自适应滤波器消除射电天文信号的射频干扰。
背景技术
RFI(Radio Frequency Interference)，射频干扰是射电望远镜在观测过程中接收到的人为、非人为信号干扰的总称。RFI会降低射电望远镜测得天文信号的信噪比，使得谱线模糊、分辨率降低，对脉冲星、中性氢等天文信号的测量有着非常大的负面影响。
随着无线通信的不断发展与射电望远镜口径的不断增加，射电望远镜越来越容易受到其他无线电信号源的干扰。此前的学者提出了许多卓有成效的方法规避RFI对射电天文信号的影响。FAST建立了电磁波宁静区，通过立法的方式划出一片无线电管制区，可以有效地避免来自通信基站和其他地面RFI源。但因为电磁波宁静区的面积受到种种条件的约束，FAST灵敏度又非常高，所以电磁波宁静区外的地面RFI还是会对射电望远镜产生影响。除此之外，电磁波宁静区也无法规避飞跃望远镜上空的人造地球卫星的射频干扰。
此外，还有学者使用SumThreshold、深度学习等方法在时频域图像中标记RFI。此类标记方法准确率较高，但并没有进一步实现RFI与有用信号的分离，没有从干扰中提取出期望测量到的天文信号。架设参考天线、建立射电望远镜所指天区的卫星RFI数据库，通过卫星轨道计算排列观测计划，从而规避卫星通过望远镜上空的时间或频率波段是一种较为成熟的被动避免RFI的方式。但随着近地轨道卫星数量的不断增加，在可预计的将来，采用被动的方式避免卫星RFI一定会使得射电望远镜的时间观察窗口和频率观察窗口大幅减少。
在科学观测过程中避免RFI是一种最直接的提高信号信噪比的方式。目前天文领域的消除RFI的方法都来自于雷达信号处理领域的自适应旁瓣对消算法，其往往需要利用到多个方向的辅助信号。
使用LMS(LeastMean Square，最小均方误差算法)自适应滤波器消除RFI是近期天文学家进行的另一种更为积极的尝试。LMS是一种最常用的自适应滤波算法，迭代滤波器的抽头系数使代价函数的均方值达到最小。与旁瓣对消算法类似，现有的LMS自适应噪声消除算法需要与天文信号同步的参考信号作为自适应滤波器的输入。具体来说，如图1所示，基于LMS自适应滤波器的RFI消除算法需要架设与望远镜同步的参考天线，没有参考信号无法使用自适应滤波器，因此对观测站硬件条件要求高。其中，a表示射电望远镜测得的天文信号，v表示干扰信号，v1和v2分别表示被射电望远镜和参考天线测得的干扰信号，因为v1和v2有共同的源v，所以v1和v2高度相关。通常可以认为天文信号与噪声信号v、v1和v2不相关。
但考虑到许多天文台站并没有架设参考天线的条件，此外，目前大多数天文信号并没有与之同步的参考信号，因此受到RFI影响而被舍弃的观测数据也无法通过这种方式处理，导致不能从现有的大量受射频干扰的信号中提取出所需要的天文信号，因此，急需提出一种不需要额外硬件设备就可以实现的自适应滤波技术。
发明内容
本发明的目的在于提供一种天文瞬时干扰抑制方法，使得大量没有同步参考信号的天文信号也可以使用最小均方误差算法进行自适应噪声抵消。
为了实现上述目的，本发明提供一种天文瞬时干扰抑制方法包括：
S1：将望远镜观测信号d通过阈值进行信号分离，得到通过阈值提取出来的噪声信号v′；
S2：将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值得到待去除的噪声信号v的估计值/>最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e，将误差函数e作为经过干扰抑制的天文信号a。
在所述步骤S1之前，还包括通过以下步骤来获取阈值：
A1：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线不再有明显减小的最大阈值作为所需阈值的第一估计值；
A2：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线不再有明显增加的最小阈值即为所需阈值的第二估计值；
A3：计算得到各个阈值候选值到所述阈值的第一估计值、阈值的第二估计值的距离之和，将距离之和最小的阈值候选值作为最终的阈值。
在所述步骤A1中，所述天文信号a的估计值等于望远镜观测信号d减去通过阈值提取出来的噪声信号v′。
在所述步骤A1中，通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2；
在所述步骤A2中，通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2。
所述阈值候选值是2的整数次幂。
所述阈值候选值包括16、32和64。
所述LMS自适应滤波器的输出信号y为：
其中，M为LMS自适应滤波器的总阶数，k表示LMS自适应滤波器的阶序数，k可取0到M-1间的整数，n表示信号采样点的序数，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，其在收敛时的值为滤波器的维纳最优抽头系数wopt，v′(n-k)表示把第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，为待去除的噪声信号v的估计值；
所述误差函数e为：
其中，e表示误差函数，为待去除的噪声信号v的估计值；
在LMS自适应滤波器中，损失函数J(n)表示为：
J(n)＝E{e(n)e*(n)}，
其中，J表示损失函数，上标*表示共轭运算，E表示对大括号内的表达式取均值，e*表示e的共轭，e是误差函数，n为信号采样点的序数；
LMS自适应滤波器的系数为：
wk(n)＝wk(n-1)+μ·v′(n一k)e*(n)，
其中，μ代表步长，是一个介于0到1之间的值，μ值的选取对自适应滤波器收敛性能至关重要，通常由经验给出；v′(n-k)表示将第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，wk(n-1)表示LMS自适应滤波器在第n个信号采样点的前一时刻的第k阶的抽头系数，e表示误差函数，上标*表示共轭运算。
所述步骤S2具体包括：
S21：选择一步长μ；
S22：在当前的步长μ下，将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值得到待去除的噪声信号v的估计值/>最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e；
S23：判断误差函数e的收敛速度和稳态误差是否均满足要求；
S24：重新选择步长μ，并回到步骤S22，直到收敛速度和稳态误差均满足要求，则输出误差函数e作为经过干扰抑制的天文信号a。
所述步长在1/400000到1/800000之间。
在所述步骤S24中，输出误差函数e在收敛后的部分作为经过干扰抑制的天文信号a。
本发明的天文瞬时干扰抑制方法将RFI与天文信号初步分离，利用初步剥离出的RFI信号与真实RFI信号的强相关性、与天文信号的弱相关性，使用最小均方误差算法实现RFI与天文信号的进一步分离。本发明的方法与现有的基于经验的分离方法和基于RFI近场假设的分离方法相比，在算法上更为简单，可以更好的消除距离遥远的人造地球卫星RFI。
(1)本发明的天文瞬时干扰抑制方法基于阈值的RFI信号粗提取方法。使用阈值法从天文信号与RFI信号的混合输入中粗提取出RFI信号。这里的RFI信号不需要非常准确，也可以包含部分天文信号。这是因为在粗提取之后还会进一步使用自适应滤波器对RFI信号进行滤波，此处提取的信号仅相当于预处理。
(2)本发明的天文瞬时干扰抑制方法基于新的RFI粗分离技术和LMS自适应滤波器，因为使用了新的RFI粗分离技术，所以无需使用参考天线数据就可以进行LMS自适应噪声抵消操作，因为对自适应滤波的输入信号只有相关性上的要求，所以粗提取的RFI信号完全可以满足算法要求。改进后的RFI处理流程能够直接使用射电望远镜输入信号进行LMS自适应滤波操作，而不需参考天线的信号输入。
(3)因为不需要与射电望远镜同步的参考天线输入信号，所以该方法更加灵活、适用性更广、具备更高的可行性。传统的自适应噪声抵消算法无法处理现存的大量天文数据，而本方法可以用于射电天文领域几乎所有的噪声抵消问题。此方法扩大了自适应噪声抵消的应用范围，大大降低了该方法对硬件设备的要求，能够处理大量已经存在的具有强射频干扰的天文数据。
(4)自适应滤波器抽头系数更新公式中的步长参数选择。一般地，给定一个平稳数据可以计算出该数据下的最优步长μ。本实验中用到的数据明显不符合平稳性假设，而且存在偏置，所以实际使用到的步长与计算出来的存在较大差异，本次实验中的步长参数选择更多是一种经验性的选择。
附图说明
图1是现有的自适应RFI消除器的原理框图。
图2是本发明的天文瞬时干扰抑制方法的原理框图。
图3是本发明的天文瞬时干扰抑制方法的流程图。
图4是新疆南山26米射电望远镜L波段的观测数据的示意图。
图5是图4的观测数据中选取第190个信道数据的示意图。
图6A和图6B是原观测数据和16阶LMS滤波器的输出数据的对比图，其中图6A示出了原观测数据，图6B示出了16阶LMS滤波器的输出数据。
图7是16阶LMS自适应滤波器最小均方误差曲线图。
图8是最优权系数误差滤波器的误差曲线图。
具体实施方式
以下结合具体实施例，对本发明做进一步说明。应理解，以下实施例仅用于说明本发明而非用于限制本发明的范围。
传统的自适应噪声消除器由雷达领域的自适应旁瓣消除器发展而来，需要参考天线提供与射电望远镜RFI相关的近场射频干扰信号作为自适应算法的迭代依据。在图1中，我们假设：
d＝a+v (1)
x＝v2 (2)
其中，d为望远镜观测信号，x为滤波器的输入信号，a为天文信号，v为待去除的噪声信号，v2为参考信号，即噪声信号。
根据此前的假设，天文信号a与待去除的噪声信号v不相关，可以计算出滤波器的输入信号x与望远镜观测信号d的互相关函数rxd、滤波器的输入信号x的自相关函数rxx：
其中，rxd为滤波器的输入信号x与望远镜观测信号d的互相关函数，rxx为滤波器的输入信号x的自相关函数，为待去除的噪声信号和参考信号的互相关函数，/>为参考信号的自相关函数。
从本质上说，待去除的噪声信号v和参考信号v2都源自噪声源的噪声信号，由于射电望远镜和参考天线在硬件上存在差异，使得接收到的信号存在差异，将这个差异抽象为一个噪声转换信道，理想情况下这个噪声转换信道可以使用FIR滤波器建模得到，进而得到噪声转换信道所对应的FIR滤波器的维纳最优抽头系数wopt(即维纳最优解)。
由此，可以计算得到，滤波器的维纳最优抽头系数wopt为：
wopt＝rxx -1·rxd (5)
在图1中可以看出，当滤波器达到维纳最优解时，参考信号v2作为滤波器的输入信号，在经过滤波后输出得到待去除的噪声信号v，滤波器的输入信号x与待去除的噪声信号v相消得到天文信号a。
综上，在传统的自适应噪声消除器中，在上述消除干扰信号的过程中需要用到参考天线测量得到滤波器的输入信号，即需要得到参考信号v2。
但在实际工作中我们常常并没有与射电望远镜信号同步的滤波器输入信号，因此，本发明提出了一种天文瞬时干扰抑制方法，其使用阈值从原信号中提取出RFI信号来代替参考信号作为滤波器的输入信号x，从而不需要参考天线来测量信号。
本发明的天文瞬时干扰抑制方法的原理框图如图2所示，其中∑符号表示求和运算，在这里表示信号与信号a+v的和。>>符号对应硬件中的逻辑右移，即，将数据的低比特位抛弃，抛弃的位数与阈值的选取有关，而阈值的选取与具体的信号有关。如图2所示，本发明的天文瞬时干扰抑制方法首先需要将望远镜观测信号a+v通过阈值提取出噪声信号v′，v′是通过阈值提取出来的噪声信号，a是天文信号，v是待去除的噪声信号；随后将通过阈值提取出来的噪声信号v′经过一个LMS(最小均方)自适应滤波器滤波，得到待去除的噪声信号v的估计值/>最终将望远镜的信号a+v减去待去除的噪声信号的估计值/>得到误差函数e，误差函数e的值被认为等于经过干扰抑制的天文信号a的值。
与上文类似，这些信号满足以下公式：
d＝a+v (6)
x＝v′ (7)
rxd＝rvv′ (8)
rxx＝rv′v′ (9)
wopt＝rxx -1·rxd (10)
其中，d为望远镜观测信号，x为滤波器的输入信号，a为天文信号，v为待去除的噪声信号，v′是通过阈值提取出来的噪声信号，rxd为滤波器的输入信号x与望远镜观测信号d的互相关函数，rxx为滤波器的输入信号x的自相关函数，wopt为LMS自适应滤波器的维纳最优抽头系数。
相应地，LMS自适应滤波器的输出信号y为：
其中，M为LMS自适应滤波器的总阶数，k表示LMS自适应滤波器的阶序数，k可取0到M-1间的整数，n表示信号采样点的序数，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，其在收敛时的值为滤波器的维纳最优抽头系数wopt，v′(n-k)中表示把第n个信号采样点的通过阈值提取出来的噪声信号延时k点，为待去除的噪声信号v的估计值。
由此，LMS自适应滤波器的输出信号y是由通过阈值提取出来的噪声信号v′作为滤波器的输入信号卷积LMS自适应滤波器的抽头系数wk(n)得到的，LMS自适应滤波器的维纳最优抽头系数wopt是通过自适应过程所期望达到的理论最优解。LMS自适应滤波器的本质是调整FIR滤波器的抽头系数，使得滤波器输出尽可能地满足期望要求。
需要注意的是，在计算滤波器的维纳最优抽头系数wopt时，需要用到所有时刻的数据，在无法做到在硬件设备上实时消除噪声，所以才需要自适应迭代过程(即LMS自适应滤波器的系数的更新过程)通过计算逼近自适应滤波器的最优解。
本发明的天文瞬时干扰抑制方法所输出的误差函数e为：
其中，e表示误差函数，是望远镜观测信号d与待去除的噪声信号v的估计值的差值，/>为待去除的噪声信号v的估计值。
在wk(n)等于LMS自适应滤波器的维纳最优抽头系数wopt时待去除的噪声信号v的估计值一般情况下可视同待去除的噪声信号v，所以经过抵消之后***输出的误差函数e可认为只包含天文信号a。
基于上述原理，如图3所示，本发明的天文瞬时干扰抑制方法包括：
步骤S1：将望远镜观测信号d通过阈值进行信号分离，得到通过阈值提取出来的噪声信号v′；
其中，在所述步骤S1中，所述阈值是经过选取的。
其中，阈值的选取标准为：
1)尽可能多地抛弃天文信号a(即，使得阈值尽可能大)，使通过阈值提取出来的噪声信号v′尽可能不包含天文信号成分，从而尽可能使公式(8)互相关计算所满足的前提，即通过阈值提取出来的噪声信号v′与天文信号a的不相关假设成立；
如上文所述，公式(8)为：rxd＝rvv′，
为了尽可能使通过阈值提取出来的噪声信号v′与天文信号a的不相关假设成立时，需要计算通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′，并且通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′越接近0越好。
天文信号a的估计值是通过以下方式来得到的：天文信号a是通过d-v^得到，由于此处其实并不需要得到精确的天文信号a，只需要将望远镜观测信号d与通过阈值提取出来的噪声信号v′的差值d-v′作为天文信号a的估计值，计算通过阈值提取出来的噪声信号v′与d-v′的相关性并作为通过阈值提取出来的噪声信号v′与天文信号a的估计值的大致的相关性即可。
在本实施例中，通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线不再有明显减小的最大阈值即为所需阈值的第一估计值。
其中，使得通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的相关性的函数曲线不再有明显减小说明通过阈值提取出来的噪声信号v′(即疑似RFI信号)不含有疑似天文信号的成分。
需要注意的是，参考信号v2与通过阈值提取出来的噪声信号v′实际上是通过不同的方法提取出的疑似RFI信号。
2)尽可能多的保留射频干扰信号(即，使得阈值尽可能小)，保留下的信号与原信号中的RFI信号相关度越大越好；也就是说，需要计算通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性，并且通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性越大越好。
在本实施例中，通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线不再有明显增加的最小阈值即为所需阈值的第二估计值。
3)因为在硬件实现过程中选取阈值对应于移位操作，所以为方便起见选取的阈值最好为2的整数次幂，如16、32、64等。
在本实施例中，原始的望远镜观测信号d是1×65536的double型数组，在硬件实现过程中，RFI为8位浮点数。
基于上述的阈值的选取标准，在所述步骤S1之前，还包括通过以下步骤来获取阈值：
步骤A1：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线收敛时的最大阈值作为所需阈值的第一估计值；
在本实施例中，通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2。函数曲线的斜率的绝对值小于0.2时达到收敛。
由此，步骤A1确保了疑似RFI信号中不含有疑似天文信号。
步骤A2：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线收敛时的最小阈值即为所需阈值的第二估计值；
通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2。
由此，步骤A2确保了疑似天文信号中不含有疑似RFI信号。
步骤A3：计算得到各个阈值候选值到所述阈值的第一估计值、阈值的第二估计值的距离之和，将距离之和最小的阈值候选值作为最终的阈值。
其中，阈值候选值是2的整数次幂，如16、32、64等。
步骤S2：将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值得到待去除的噪声信号v的估计值/>最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e，将误差函数e作为经过干扰抑制的天文信号a。
其中，在LMS自适应滤波器的抽头系数wk(n)等于LMS自适应滤波器的维纳最优抽头系数wopt时待去除的噪声信号v的估计值一般情况下可视同待去除的噪声信号v，所以经过抵消之后***输出的误差函数e可认为只包含天文信号a。
具体的计算公式可以参见上述的公式(6)至公式(12)，即，
d＝a+v (6)
x＝v′ (7)
rxd＝rvv′ (8)
rxx＝rv′v′ (9)
wopt＝rxx -1·rxd (10)
其中，d为望远镜观测信号，x为滤波器的输入信号，a为天文信号，v为待去除的噪声信号，v′是通过阈值提取出来的噪声信号，rxd为滤波器的输入信号x与望远镜观测信号d的互相关函数，rxx为滤波器的输入信号x的自相关函数，wopt为LMS自适应滤波器的维纳最优抽头系数。
LMS自适应滤波器的输出信号y为：
其中，M为LMS自适应滤波器的总阶数，k表示LMS自适应滤波器的阶序数，k可取0到M-1间的整数，n表示信号采样点的序数，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，其在收敛时的值为滤波器的维纳最优抽头系数wopt，v′(n-k)表示把第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，为待去除的噪声信号v的估计值。
误差函数e为：
其中，e表示误差函数，为待去除的噪声信号v的估计值。
在LMS自适应滤波器中，损失函数J(n)可表示为：
J(n)＝E{e(n)e*(n)} (13)
其中，J表示损失函数，上标*表示共轭运算，E表示对大括号内的表达式取均值，e*表示e的共轭，e是误差函数，e(n)代表第n个采样点的误差，n为信号采样点的序数。
需要说明的是，d是向量，也是向量，数据长度等于采样点数，所以误差函数e也是一个向量，e(n)代表第n个采样点的误差，即/>
LMS自适应滤波器的系数的更新过程由以下公式表示：
wk(n)＝wk(n-1)+μ·v′(n-k)e*(n) (14)
其中，μ代表步长，是一个介于0到1之间的值，μ值的选取对自适应滤波器收敛性能至关重要，通常由经验给出；v′(n-k)表示将第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，wk(n-1)表示LMS自适应滤波器在第n个信号采样点的前一时刻的第k阶的抽头系数，e表示误差函数，e(n)代表第n个采样点的误差，上标*表示共轭运算。
如上文所述，滤波器的维纳最优抽头系数wopt是由理想维纳解的滤波器系数计算公式，即公式(10)得到，LMS滤波器的计算公式由公式(14)计算得到。而因为维纳解的计算需要用到整个观测的数据，无法做到实时消除，而LMS滤波器计算公式只需要用到当前与过去一段时间的数据，因此本发明必须采用公式(14)来获取LMS自适应滤波器的抽头系数，进而逐渐逼近滤波器的维纳最优抽头系数wopt。滤波器的维纳最优抽头系数wopt是LMS滤波器计算的滤波器抽头系数的比较基准，是LMS滤波器系数与wopt越相似越好。
在本实施例中，误差函数e是逐渐收敛的函数，因此，所述步骤S2具体包括：
步骤S21：选择一步长μ；
步骤S22：在当前的步长μ下，将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值得到待去除的噪声信号v的估计值/>最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e；
步骤S23：判断误差函数e的收敛速度和稳态误差是否均满足要求；
其中，可认为权系数误差曲线斜率的绝对值小于0.2时达到收敛，或观察MSE(最小均方误差)曲线开始围绕某一固定的值上下波动而不是有固定的方向时达到收敛。
收敛速度可以根据MSE与最优权系数误差曲线的下降所需的迭代次数(或者达到稳态的迭代次数)看出，也可以观察自适应***对输入中的瞬时脉冲干扰信号的抑制情况得到大概的收敛速度。稳态误差可以根据MSE与最优权系数误差的收敛值得到，也可以观察***输出信号的波形判断。在自适应滤波算法中不存在特定的收敛速度计算公式，因为收敛速度可以从图像中非常直观的看出来。应用最广的稳态误差的计算公式就是此前介绍的MSE与最优权系数误差。收敛速度和稳态误差没有严格的数值标准，一般都是在具体的情况中具体地分析。比如，对于瞬时脉冲干扰来说如果收敛速度不够快，可能不等***收敛干扰就已结束，所以这时的收敛速度需要短于脉冲信号的持续时间。稳态误差就更不存在数值范围了，每种情况对数据的要求可能也会不太一样，需要根据情况具体确定，在本发明中，误差函数e是作为输出信号的。
步骤S24：重新选择步长μ，并回到步骤S22，直到收敛速度和稳态误差均满足要求，则输出误差函数e作为经过干扰抑制的天文信号a。
优选地，输出误差函数e在收敛后的部分作为经过干扰抑制的天文信号a。
实验结果：
实验数据使用2019年9月3日新疆南山26米射电望远镜L波段的观测数据，如图4所示。
望远镜观测数据为65536×512矩阵，其中数据采样频率为共65536个采样点，频率范围从1300Mhz到1811Mhz共512Mhz，512个信道每个信道带宽为1Mhz，第一个信道对应最高频率成分，最后一个信道对应最低频率成分。此处选择第190个信道进行最小均方误差自适应滤波处理。第190个信道数据如图5所示。
图5中持续时间较短、幅度较大的为RFI，天文信号蕴含在底部30至40附近波动的值中，目的是去除样本中的瞬时脉冲干扰信号。
此时的输入信号相当于图2中的a+v，选取阈值为32、步长μ为2.5×10-6、FIR滤波器的阶数为16，构建LMS自适应RFI抵消***。
经过matlab程序仿真之后得到使用16阶LMS自适应滤波器的滤波前后对比图如图6A和图6B所示。
图6A和图6B是原观测数据和16阶LMS滤波器的输出数据的对比图，其中图6A示出了原观测数据，图6B示出了16阶LMS滤波器的输出数据，图中横坐标表示迭代次数，纵坐标表示信号幅值。
通过图6A和图6B对比可以看出，在滤波之后绝大部分瞬时脉冲信号被很好的抑制了。但是根据图6A和图6B可以明显地看出，在自适迭代的起始位置有许多干扰没有被抑制，这是因为自适应滤波器达到最优效果需要一个收敛过程，这个收敛过程与步长μ的选取有关。一般地，μ越大收敛速率越快，但随着μ的增大稳态误差也会增大。而且μ值过大会使得自适应滤波器不收敛，μ值存在上限，这个上限与信号本身的特征有关。此处的μ值是综合考虑各个性能指标后选择的。同时需要特别注意的是，在原信号中，天文信号存在一个偏置，天文信号均值在30附近。这个偏置不利于自适应滤波器的迭代优化，而且天文信号本身蕴含在信号波形中而不是均值中，所以这里的滤波器输出信号去除了这个偏置，信号在0值附近波动。
本发明所采用的最小化均方误差(Mean Square Error，MSE)的LMS自适应滤波器和作为对比的权系数误差滤波器两种方法量化滤波器的滤波结果。
LMS自适应滤波器的MSE(均方误差)的计算公式由公式(15)给出：
其中，MSE(n)是LMS自适应滤波器的均方误差，M为滤波器阶数，n为表示信号采样点的序数，由于每次自适应迭代过程输入一个采样点信号，所以n也可以表示迭代次数，最大值为采样点个数，e(n)代表第n个采样点的误差，由期望信号d与滤波器输出相减后第n个采样点的值得到。
图7是16阶LMS自适应滤波器最小均方误差曲线，横坐标表示迭代次数，纵坐标表示均方误差，单位是分贝(db)。图7所展现出来的最小均方误差图像显示，误差收敛趋势十分明显，稳态误差较小。这表明该自适应过程能够较快地达到所期望的输出响应信号。曲线收敛后的部分波动是瞬时脉冲干扰信号造成的，现阶段还无法完全去除，误差波动总体处于一个可以接受的水平上。
使用权系数误差是通过当前时刻自适应滤波器抽头系数与滤波器的维纳最优抽头系数差值的平方来描述误差的一种方法，见公式(16)：
其中，k表示LMS自适应滤波器的阶序数，k可取0到M-1间的整数，wk(n)表示LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，由式(16)计算得出，wopt，k为自适应滤波器的第k阶的维纳最优抽头系数，其由公式(10)计算得出。
使用最优权系数误差滤波器生成的误差曲线图如图8所示，横坐标表示迭代次数，纵坐标表示与最优权系数误差滤波器的权系数误差，单位是分贝(db)。从图8中可以看出，最优权系数误差曲线收敛较快，稳态误差较小。这表明该自适应滤波器能够通过自适应迭代过程较快地达到最优维纳滤波器对应的系数。然而，需要指出的是，在式(10)中wopt的计算过程中需要计算自相关矩阵rxx的逆矩阵。矩阵求逆运算是一个非常复杂的运算，尤其当矩阵的维度非常大时，对其求逆几乎是一个不可能完成的任务。此处只选择前1024个数据用于最优滤波器系数的计算。
本发明使用基于阈值分离RFI信号的LMS自适应噪声消除器能在没有额外参考天线的输入下达到良好的滤波效果，同时保留最多的天文信号。在相同情况下，现有的其他RFI处理方法均无法达到同样效果。
(1)本发明的天文瞬时干扰抑制方法基于阈值的RFI信号粗提取方法。使用阈值法从天文信号与RFI信号的混合输入中粗提取出RFI信号。这里的RFI信号不需要非常准确，也可以包含部分天文信号。这是因为在粗提取之后还会进一步使用自适应滤波器对RFI信号进行滤波，此处提取的信号仅相当于预处理。
(2)本发明的天文瞬时干扰抑制方法基于新的RFI粗分离技术和LMS自适应滤波器，因为使用了新的RFI粗分离技术，所以无需使用参考天线数据就可以进行LMS自适应噪声抵消操作，因为对自适应滤波的输入信号只有相关性上的要求，所以粗提取的RFI信号完全可以满足算法要求。改进后的RFI处理流程能够直接使用射电望远镜输入信号进行LMS自适应滤波操作，而不需参考天线的信号输入。
(3)因为不需要与射电望远镜同步的参考天线输入信号，所以该方法更加灵活、适用性更广、具备更高的可行性。传统的自适应噪声抵消算法无法处理现存的大量天文数据，而本方法可以用于射电天文领域几乎所有的噪声抵消问题。此方法扩大了自适应噪声抵消的应用范围，大大降低了该方法对硬件设备的要求，能够处理大量已经存在的具有强射频干扰的天文数据。
(4)自适应滤波器抽头系数更新公式中的步长参数选择。一般地，给定一个平稳数据可以计算出该数据下的最优步长μ。本实验中用到的数据明显不符合平稳性假设，而且存在偏置，所以实际使用到的步长与计算出来的存在较大差异，本次实验中的步长参数选择更多是一种经验性的选择。步长根据信号以及信号预处理的方式的不同会有很大的差异，差异有时甚至能够达到五六个数量级，在本次实验中所尝试的步长在1/400000到1/800000之间。
以上所述的，仅为本发明的较佳实施例，并非用以限定本发明的范围，本发明的上述实施例还可以做出各种变化。凡是依据本发明申请的权利要求书及说明书内容所作的简单、等效变化与修饰，皆落入本发明专利的权利要求保护范围。本发明未详尽描述的均为常规技术内容。
Claims (10)
1.一种天文瞬时干扰抑制方法，其特征在于，包括：
步骤S1：将望远镜观测信号d通过阈值进行信号分离，得到通过阈值提取出来的噪声信号v′；
步骤S2：将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e，将误差函数e作为经过干扰抑制的天文信号a。
2.根据权利要求1所述的天文瞬时干扰抑制方法，其特征在于，在所述步骤S1之前，还包括通过以下步骤来获取阈值：
步骤A1：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线收敛时的最大阈值作为所需阈值的第一估计值；
步骤A2：通过多次改变阈值，计算通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线并观察该函数曲线，使通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线收敛时的最小阈值即为所需阈值的第二估计值；
步骤A3：计算得到各个阈值候选值到所述阈值的第一估计值、阈值的第二估计值的距离之和，将距离之和最小的阈值候选值作为最终的阈值。
3.根据权利要求2所述的天文瞬时干扰抑制方法，其特征在于，在所述步骤A1中，所述天文信号a的估计值等于望远镜观测信号d减去通过阈值提取出来的噪声信号v′。
4.根据权利要求2所述的天文瞬时干扰抑制方法，其特征在于，在所述步骤A1中，通过阈值提取出来的噪声信号v′与天文信号a的估计值的相关性rav′的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2；
在所述步骤A2中，通过阈值提取出来的噪声信号v′与望远镜观测信号d的相关性的函数曲线收敛是指函数曲线的斜率的绝对值小于0.2。
5.根据权利要求2所述的天文瞬时干扰抑制方法，其特征在于，所述阈值候选值是2的整数次幂。
6.根据权利要求5所述的天文瞬时干扰抑制方法，其特征在于，所述阈值候选值包括16、32和64。
7.根据权利要求1所述的天文瞬时干扰抑制方法，其特征在于，所述LMS自适应滤波器的输出信号y为：
其中，M为LMS自适应滤波器的总阶数，k表示LMS自适应滤波器的阶序数，k可取0到M-1间的整数，n表示信号采样点的序数，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，其在收敛时的值为滤波器的维纳最优抽头系数wopt，v′(n-k)表示把第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，为待去除的噪声信号v的估计值；
所述误差函数e为：
其中，e表示误差函数，为待去除的噪声信号v的估计值；
在LMS自适应滤波器中，损失函数J(n)表示为：
J(n)＝E{e(n)e*(n)}，
其中，J表示损失函数，上标*表示共轭运算，E表示对大括号内的表达式取均值，e*表示e的共轭，e是误差函数，e(n)代表第n个采样点的误差，n为信号采样点的序数；
LMS自适应滤波器的系数为：
wk(n)＝wk(n-1)+μ·v′(n-k)e*(n)，
其中，μ代表步长，是一个介于0到1之间的值；v′(n-k)表示将第n个信号采样点的通过阈值提取出来的噪声信号v′延时k点，wk(n)为LMS自适应滤波器在第n个信号采样点的第k阶的抽头系数，wk(n-1)表示LMS自适应滤波器在第n个信号采样点的前一时刻的第k阶的抽头系数，e表示误差函数，e(n)代表第n个采样点的误差，上标*表示共轭运算。
8.根据权利要求1所述的天文瞬时干扰抑制方法，其特征在于，所述步骤S2具体包括：
步骤S21：选择一步长μ；
步骤S22：在当前的步长μ下，将通过阈值提取出来的噪声信号v′经过一个LMS自适应滤波器滤波，得到待去除的噪声信号v的估计值得到待去除的噪声信号v的估计值/>最终将望远镜观测信号d减去待去除的噪声信号的估计值/>得到误差函数e；
步骤S23：判断误差函数e的收敛速度和稳态误差是否均满足要求；
步骤S24：重新选择步长μ，并回到步骤S22，直到收敛速度和稳态误差均满足要求，则输出误差函数e作为经过干扰抑制的天文信号a。
9.根据权利要求8所述的天文瞬时干扰抑制方法，其特征在于，所述步长在1/400000到1/800000之间。
10.根据权利要求8所述的天文瞬时干扰抑制方法，其特征在于，在所述步骤S24中，输出误差函数e在收敛后的部分作为经过干扰抑制的天文信号a。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311106072.XA CN117335772A (zh) | 2023-08-30 | 2023-08-30 | 一种天文瞬时干扰抑制方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311106072.XA CN117335772A (zh) | 2023-08-30 | 2023-08-30 | 一种天文瞬时干扰抑制方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117335772A true CN117335772A (zh) | 2024-01-02 |
Family
ID=89289224
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311106072.XA Pending CN117335772A (zh) | 2023-08-30 | 2023-08-30 | 一种天文瞬时干扰抑制方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117335772A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118033549A (zh) * | 2024-04-15 | 2024-05-14 | 华中科技大学 | 一种自适应盲源分离的遥感图像射频干扰去除方法及*** |
-
2023
- 2023-08-30 CN CN202311106072.XA patent/CN117335772A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118033549A (zh) * | 2024-04-15 | 2024-05-14 | 华中科技大学 | 一种自适应盲源分离的遥感图像射频干扰去除方法及*** |
| CN118033549B (zh) * | 2024-04-15 | 2024-06-28 | 华中科技大学 | 一种自适应盲源分离的遥感图像射频干扰去除方法及*** |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110634500B (zh) | 一种先验信噪比的计算方法、电子设备及存储介质 | |
| US6157909A (en) | Process and device for blind equalization of the effects of a transmission channel on a digital speech signal | |
| CN101729461B (zh) | 消除单频干扰及多频干扰的***及方法 | |
| CN109977914B (zh) | 基于vmd的自适应降噪方法 | |
| CN104835503A (zh) | 一种改进gsc自适应语音增强方法 | |
| CN117335772A (zh) | 一种天文瞬时干扰抑制方法 | |
| CN104808219A (zh) | 一种新型的空时联合抗干扰方法 | |
| Meller | Cheap cancellation of strong echoes for digital passive and noise radars | |
| CN110926455B (zh) | 一种射电天文信号的自适应射频干扰消除方法 | |
| Lobov et al. | Dispersion distortion tracking compensator based on the sigma-point Kalman | |
| CN114938232B (zh) | 基于lstm的同时同频全双工数字域自干扰抑制方法 | |
| CN116155306A (zh) | 一种基于维纳滤波和自适应滤波算法的磁电联合低频信号接收机及信号接收方法 | |
| EP3712626B1 (en) | High-rate dft-based data manipulator and data manipulation method for high performance and robust signal processing | |
| Banchhor et al. | GUI based performance analysis of speech enhancement techniques | |
| CN110444222B (zh) | 一种基于信息熵加权的话音降噪方法 | |
| US11888538B2 (en) | Receiving apparatus | |
| CN103929150B (zh) | 一种子带自适应滤波器的权值向量更新方法 | |
| CN113824488A (zh) | 基于判决反馈自适应对消的卫星通信非恶意干扰抑制方法 | |
| Wang | Direct signal recovery and masking effect removal exploiting sparsity for passive bistatic radar | |
| Lampl | Implementation of adaptive filtering algorithms for noise cancellation | |
| CN108375760B (zh) | 雷达杂波抑制的fir滤波器设计方法 | |
| George et al. | Multiple signal detection and measurement using a configurable wideband digital receiver | |
| Soni et al. | Design, performance and cost analysis of various band pass IIR filters for Myriametre band applications | |
| KR101784607B1 (ko) | 다중 레이더 시스템에서의 간섭 신호 제거 장치 및 방법 | |
| CN112769414B (zh) | 一种基于自适应滤波器的干扰信号处理方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |