CN109581422A - 一种基于子空间误差修正的线性扫频干扰抑制方法及系统 - Google Patents
一种基于子空间误差修正的线性扫频干扰抑制方法及系统 Download PDFInfo
- Publication number
- CN109581422A CN109581422A CN201811540877.4A CN201811540877A CN109581422A CN 109581422 A CN109581422 A CN 109581422A CN 201811540877 A CN201811540877 A CN 201811540877A CN 109581422 A CN109581422 A CN 109581422A
- Authority
- CN
- China
- Prior art keywords
- frequency
- interference
- time
- signal
- subspace
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/21—Interference related issues ; Issues related to cross-correlation, spoofing or other methods of denial of service
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
- Noise Elimination (AREA)
Abstract
本发明公开了一种基于子空间误差修正的线性扫频干扰抑制方法,包括步骤1:对用阵列天线接收到的信号白化处理并进行时频分析得到优化后的接收信号时频分布图;步骤2:从优化后的时频分布图中遴选自源点;步骤3:通过自源点估计线性扫频频率参数即初始频率fi和调频斜率gi,i为第i个线性扫频干扰信号,i=2,…,P,P为干扰信号的数量;步骤4:对所估计的频率参数进行修正;步骤5:使用修正后的频率参数构建干扰子空间,用子空间投影消除干扰。本发明通过对估计的频率参数进行修正,构建更准确的干扰子空间,然后通过子空间投影的方法消除干扰取得了较好的抑制效果,避免了在构建子空间基矢量序列长度过长和过短都损伤子空间投影抗干扰性能。
Description
技术领域
本发明属于信号处理领域,尤其涉及一种子空间误差修正的线性扫频干扰抑制方法及系统。
背景技术
线性扫频(LFM:Linear Frequency Modulation)干扰是战场环境下普遍的一种干扰,属于非平稳干扰,使用处理平稳干扰的方法处理效果不好,文献1“D C Ricks,P GCifuentes,J S Goldstein.What is Optimal Processing for Nonstationary Data[C].Pacific Grove,CA:34th Asilomar Conference on Signals,Systems,and Computers,2000,11,656-661”对此进行了分析和仿真验证。采用子空间投影的方法,在子空间基矢量构建精准的情况下,可以取得很好的效果,如文献2“Moeness G Amin,Zhao L,Lindsey AR.Subspace Array Processing for the Suppression of FM Jamming in GPSReceivers[J].IEEE Trans on Aerospace and Electronic Systems,2004,40(1):80-92”中的仿真结果。构建LFM干扰子空间基矢量的关键问题在于其参数的估计。文献2和文献3“周柱,卢树军,张尔扬等.GPS接收信号中线性扫频干扰的抑制[J].国防科技大学学报,2013,35(1):70-76”首先估计LFM干扰的参数,据此构建干扰子空间,在此基础上用子空间投影法消除干扰。文献4“李利,司锡才,张雯雯等.改进的多分量LFM信号参数估计算法及其快速实现[J].系统工程与电子技术,2009,31(11):2560-2562”通过二阶非线性变换在信号参数空间形成最大值来估计LFM信号参数,在多分量的情况下,讨论了信号自项和交叉项与时间的关系,发现自项和交叉项对时间有不同的依赖性,提出了加权平均的方法来改进算法,克服了交叉项的影响。
LFM参数估计必然存在误差,如果采用现有的子空间投影法构造干扰子空间,并用子空间投影抑制干扰,微小的估计偏差会带来输出信干噪比极大的下降。研究发现用于构造子空间的序列长度对于基于子空间投影抗干扰方法的性能有很大影响,具体体现为:1、用较短的序列缩短干扰子空间的基矢量,能够极大提高子空间投影方法对于参数估计的容错性,提高了抗干扰系统的鲁棒性;2、用较短序列相当于给信号加窗,窗口越短,信号频谱泄漏越严重。因为构建子空间基矢量序列长度过长和过短都损伤子空间投影抗干扰性能,所以需要确定最佳的序列长度,而最佳序列长度难于用显式求得。本发明将从误差修正的角度解决这一问题,即通过修正LFM干扰信号的频率参数误差,从而精确构建干扰子空间,完成干扰抑制。
本发明主要用到的背景技术是非平稳信号矢量的时频分析。
以M元天线阵接收GPS信号为例进行分析,采用何种导航信号对于抗干扰而言,处理方法相同。接收机接收到的信号由信号、白噪声和干扰组成。对于用阵列天线接收到的信号,其数学形式可以表示为
x(t)=as·s(t)+Ay(t)+v(t) (1)
其中s(t),as分别表示GPS信号及其天线阵响应;A为混合矩阵,A=[A1a(θ1),A2a(θ2),…,APa(θP)]T,其中Ai、a(θi)、θi分别表示第i个干扰的振幅、在阵列上的响应即导向矢量、干扰入射角,设有P个干扰从不同方向入射,y(t)表示LFM干扰的合集,只包含信号的相位,幅度已经归一化,可写为y(t)=[y1(t),y2(t),…,yP(t)]T;v(t)是均值为零的加性高斯白噪声。
取接收信号矢量的任意一路信号作时频分析,即可完成LFM信号的参数估计。时频分析能够反映信号的时频局部性质,时频分析具有代表性的两种方法是Wigner-Ville变换和短时傅里叶变换。信号的Wigner-Ville变换生成的时频分布称为Wigner-Ville分布(WVD),WVD具有良好的局部性质,是一种最基本、也是应用最多的时频分布,取接收信号矢量的任意一路,设为x(t),其Wigner-Ville变换如下
WVD具有理想的时频聚集性,但有多分量时,交叉项严重。干扰信号的本身在时频分布图上的体现称之为自源点,不同干扰信号在时频分布图上的体现惩治为互源点。用Wigner-Ville变换处理多个的非平稳信号,在时频分布图中,自源点和互源点均存在,而且互源点大量存在。通过对接收信号矢量进行白化可以消除时频分布的交叉项,体现在时频分布图上也就是消除互源点。
当接收信号为矢量形式时,矢量各分量的自Wigner-Ville变换和分量之间的互Wigner-Ville变换形成了一个M×M维的矩阵,这就是文献5“Adel Belouchrani,MoenessG.Amin.Time-frequency MUSIC[J].IEEE Signal Processing Letters,1999,6:109-110.”中描述的空时频分布(STFD:Spatial Temporal Frequency Distribution)矩阵。STFD矩阵的离散表示式为
假设采样间隔为Ts,式中n代表时间nTs,ω为归一化频率。
按照(3)式求取接收信号矢量对应的STFD矩阵,(3)式中所取的是无限长度,实际只需要取有限长序列,如下式所示。
由上式可知,STFD矩阵是一个矩阵序列,每一个时间和频率点对(n,ω)对应一个M×M矩阵。
发明内容
本发明要解决的技术问题是针对现有技术中使用子空间投影方法消除非平稳干扰时,因为线性扫频参数估计存在误差,微小的估计偏差会带来输出信干噪比极大的下降的问题。对此提出一种可精确构建干扰子空间,从而有效抑制干扰的一种基于子空间误差修正的线性扫频干扰抑制方法。
为解决该问题,本发明所采用的技术方案是:
一种基于子空间误差修正的线性扫频干扰抑制方法,包括以下步骤:
步骤1:对用阵列天线接收到的信号白化处理并进行时频分析得到优化后的接收信号时频分布图;
步骤2:从优化后的时频分布图中遴选自源点;
步骤3:通过自源点估计线性扫频频率参数即初始频率fi和调频斜率gi,i为第i个线性扫频干扰信号,i=2,…,P,P为干扰信号的数量;
步骤4:对所估计的频率参数进行修正;
步骤5:使用修正后的频率参数构建干扰子空间,用子空间投影消除干扰。
进一步地,步骤1中得到时频分布图的具体方法为:
步骤1.1:求取接收信号矢量的白化矩阵W;
步骤1.2:用白化矩阵将接收信号矢量白化并用Wigner-Ville变换求取时频分布,得到优化后的接收信号时频分布图。
进一步地,步骤2中遴选自源点的方法为:
步骤2.1:对白化后的接收信号矢量z(n)求取空时频分布矩阵序列,设时频分析中所取的数据段长度为L,如下式
n表示时间点,ω表示归一化频率,l表示偏移量,Dzz(n,ω)表示白化后的接收信号矢量z(n)的空时频分布矩阵序列。
步骤2.2:求得时频分布图上所有时频点对应的空时频分布STFD矩阵的迹并归一化,即Trace(E{Dzz(n,ω)})/L,设ζ=Trace(E{Dzz(n,ω)})/L;
步骤2.3:设置门限值ε,当STFD矩阵的迹ζ>ε时,则此矩阵对应的时频点(n,ω)属于自源点,0.1≤ε≤0.3。
进一步地,步骤3通过自源点估计线性扫频频率调制参数即初始频率fi和调频斜率gi的方法为;
步骤3.1:寻找第i个线性扫频干扰对应时频脊线的邻域;在第i个线性扫频干扰信号对应的时频分布脊线上任取一个自源点,过该点作一条直线,使倾角α在0°~180°中变化,设定直线的邻域半径rline,所述邻域半径rline不超过干扰在时频分布图的时频脊线在任意时间点上频宽的一半,将时频分布图上各自源点到该直线的距离小于rline的自源点集合做为该直线的邻域Ωi,统计不同倾角下直线的邻域Ωi中自源点的数量,当某一倾角下直线的邻域Ωi中自源点的数量比该直线处于其它倾角时自源点的数量都大时,则将该直线倾角作为第i个线性扫频干扰信号在时频分布面上的倾角,根据该倾角计算出第i个线性扫频干扰信号的调频斜率;寻找第i个线性扫频干扰对应时频脊线的邻域,i=2,…,P,P为干扰信号的数量;在第i个线性扫频干扰信号对应的时频分布脊线上任取一个自源点,过该点作一条直线,使倾角α在0°~180°中变化,设定直线的邻域半径rline,将时频分布图上各自源点到该直线的距离小于rline的自源点集合做为该直线的邻域Ωi,统计不同倾角下直线的邻域Ωi中自源点的数量,当某一倾角下直线的邻域Ωi中自源点的数量比该直线处于其它倾角时自源点的数量大时,则将该直线倾角作为在第i个线性扫频干扰信号在时频分布面上的倾角,根据该倾角计算出第i个线性扫频干扰信号的调频斜率;
步骤3.2:去除第i个线性扫频干扰对应时频脊线的邻域与其它线性扫频干扰对应时频脊线的邻域的交集,得到去除交集后的区域Ω′i,用最小二乘法拟合Ω′i邻域中的自源点,得到第i个干扰信号的拟合直线方程则为估计出的第i个线性扫频干扰信号的初始频率,为估计出的调频斜率。
进一步地,所述步骤4中对所估计的频率调制参数进行修正的具体方法为:
步骤4.1:求取根据所估计的频率参数构建的LFM干扰信号和接收信号的相关值;;
所述相关值是指:对用估计参数构建的LFM干扰信号与接收信号的某一路x(t)求内积
式中T为求内积的时间长度,其离散形式为:
L为基矢量序列长度;
步骤4.2:以估计值为中心,选取周边点;
设定初始频率和调频斜率的步进值δf和δg,在以初始频率f为横轴,调频斜率g为纵轴的平面上,以为中心,选取其左右上下四个点,即
步骤4.3:求取各个点对应的内积
将周围的四个点对应的频率参数值代入(1.2)式当中,计算相应的内积值
步骤4.4:更新中心点
比较中心点和周围四个点对应内积的大小,用五个点中对应内积最大的点取代中心点,返回步骤4.2,如果中心点对应的内积已经是最大,则停止运算,输出该内积最大值所对应的点坐标,该坐标为修正后的初始频率和调频斜率。
进一步地,步骤5中使用修正后的频率参数构建干扰子空间,用子空间投影消除干扰的具体方法为:
步骤5.1:根据修正后的频率参数构建干扰信号序列;
用修正后的频率参数构建的信号为 取xi(n)序列的L长度数据作为子空间基矢量,如下
步骤5.2:构建干扰子空间及其正交子空间。
用时域扩展后的干扰信号矢量构建成一个矩阵,如下
B=[X1 X2 … XP] (9)
该组矢量张成了干扰的子空间,其正交子空间表示如下
步骤5.3:进行子空间投影消除干扰。
将接收信号矢量投影到干扰子空间的正交子空间,如下式
Y=PX=PXs+PV (11)
一种基于子空间误差修正的线性扫频干扰抑制系统,包括处理器和存储器,所述存储器中存储有基于子空间误差修正的线性扫频干扰抑制方法的程序,所述处理器在运行所述基于子空间误差修正的线性扫频干扰抑制方法的程序时执行抑制方法的各步骤。
与现有技术相比,本发明取得的有益效果为:
本发明一种基于子空间误差修正的线性扫频干扰抑制方法及系统,研究发现,通过构造干扰子空间进行子空间投影消除干扰时,用于构造子空间的基矢量长度较大时,由于干扰频率参数估计的微小误差将导致干扰不能被抑制,而在基矢量长度取值较小,虽然干扰频率参数估计有微小误差时,干扰也能被消除,但是子空间基矢量长度过短相当于对信号进行加窗,窗口很窄,频谱泄露越严重,所以子空间基矢量长度不是越短越好,因此本发明通过对估计的频率参数进行修正,构建更准确的干扰子空间,然后通过子空间投影的方法消除干扰取得了较好的抑制效果,避免了在构建子空间基矢量序列长度过长和过短都损伤子空间投影抗干扰性能。
附图说明
图1为时频分布的优化效果图;
图2为不同倾角的直线邻域内自源点数目;
图3为频率参数偏差对于所构建信号与接收信号的相关值的影响;
图4为三种情况下抗干扰效果对照。
具体实施方式
为了更好的说明本发明,本发明设置一个实验场景,接收信号中含有3个非平稳干扰,不妨设为LFM干扰信号,设为xi(t)=Aiexp(j(2πfit+gi/2t2)),i=1,2,3,其中f1=0.1,g1=0.1,f2=0.3,g2=-0.15,f3=0.45,g3=-0.1。干扰入射方向分别为40°、80°、120°。图1至图4示出了本发明一种基于子空间误差修正的线性扫频干扰抑制方法及系统的一种具体实施例。包括以下步骤:
步骤1:对用阵列天线接收到的信号白化处理并进行时频分析得到优化后的接收信号时频分布图;
步骤1.1:求取接收信号矢量的白化矩阵W;
若将接收信号矢量白化,接收信号采用(1)式的数学形式,则不含白噪声的接收信号矢量记为y(t),设其对应白化矩阵为W,按定义,白化矩阵应写为
上式中为纯粹非平稳信号的自相关矩阵。
上式表明如果W为白化矩阵,则WA为一个酉矩阵,令该酉矩阵为U,则有
A=W#U (13)
式中W#表示W的Moore-Penrose逆矩阵,白化过程中将阵列流形矩阵转化为了酉矩阵,矩阵维数从阵列的维数M降到了子空间维数P。
下面推导白化矩阵W的求解方式,首先对含有白噪声的接收信号求自相关矩阵,有
对于纯粹信号,不考虑白化,则有
结合(14)(15)两式,有
AAH=Rx-σ2I (16)
式中σ2为白噪声方差。因为A=W#U,则有
AAH=(W#U)(W#U)H=W#(W#)H (17)
将(17)代入(16)式,有
式中Λ是由Rx的特征值组成的对角矩阵,即Λ=diag[λ1,λ2,,λP].令上式中的对角矩阵可以推出
W#=UK (19)
根据定义,白化矩阵W与其Moore-Penrose逆矩阵W#存在关系:W#WW#=W#。将(19)式代入该关系有
(UK)W(UK)=UK (20)
因为U为酉矩阵,上式最终可以化为
W=K-1UH (21)
设U=[u1,u2,…,uP],则有
可见求解白化矩阵关键是接收信号自相关矩阵的特征值和特征矢量构成的酉矩阵。据以上分析将白化矩阵的求解过程总结如下:第一步,协方差矩阵估计
第二步,计算矩阵的特征值及对应特征矢量,将特征值按降序排列可写为:λ1≥λ2≥…≥λP;相应的特征矢量为:u1,u2,…,uP。
第三步,估计噪声方差
第四步,构造如下的对角矩阵
第五步,构造酉矩阵
第六步,计算白化矩阵。
至此求得白化矩阵,可以用该白化矩阵对接收信号矢量进行白化。
步骤1.2:用白化矩阵将接收信号矢量白化并用Wigner-Ville变换求取时频分布,得到优化后的接收信号时频分布图;
求得白化矩阵之后,将接收信号矢量白化,将白化后的接收信号矢量记为z(n),则有
z(n)=Wx(n)=W(Ay(n)+v(n))=Φy(n)+Wv(n) (26)
将接收信号矢量白化之后,再用Wigner-Ville变换求取时频分布,可以得到优化后的时频分布图。可以观察白化的效果,采用场景一中的设置,结果如附图1所示,Wigner-Ville分布的时频聚集性好,但是存在严重的交叉项,所以时频分布图中有大量的互源点,通过白化可以消除时频分布的交叉项,也即消除了时频分布图中的互源点,突出自源点。
非平稳干扰的时频特性体现在时频分布图中,有几种方法可以得到时频分布图,典型的有Wigner-Ville变换和短时傅里叶变换,其中Wigner-Ville变换时频聚集性更好,但缺点在于交叉项严重。对于这一问题,文献7“Adel Belouchrani,Moeness G.Amin.BlindSource Separation Based on Time–Frequency Signal Representations[J].IEEETransactions on Signal Processing,1998,46(11):2888-2897”提出一种方法,推导出白化矩阵,并用于对接收信号矢量进行白化,对白化后的接收信号矢量用Wigner-Ville变换进行时频转换,将消除掉原本在Wigner-Ville变换之后在时频分布图中存在的互源点,此时干扰信号的自源点清晰地显现出来。
步骤2:从优化后的接收信号时频分布图中遴选自源点;
将接收信号矢量白化之后,也就消除了其时频分布的交叉项,虽然消除了交叉项,但是在时频分布图中选出自源点还需要做进一步处理,选出干扰信号的自源点,步骤如下。
步骤2.1:对白化后的接收信号矢量z(n)求取STFD矩阵序列,设时频分析中所取的数据段长度为L,如下式
步骤2.2:求得时频分布图中所有时频点对应的STFD矩阵的迹并归一化,即Trace(E{Dzz(n,ω)})/L,设ζ=Trace(E{Dzz(n,ω)})/L
步骤1.3.3:设置门限值ε,上述空时频分布矩阵的迹均已归一化,其值域为(0,1),当空时频分布矩阵的迹ζ>ε时,则此矩阵对应的时频点(n,ω)属于自源点,0.1≤ε≤0.3,本实施例中ε取为0.3;
当STFD矩阵的迹ζ>ε时,认为此矩阵对应的时频点属于干扰信号的自项,即得到自源点。部分矩阵的迹小于ε的自项对应的时频点未选入对结果影响很小,可以忽略不计。通过本步骤可以清晰的从时频分布图中取得自源点。
步骤3:通过自源点估计线性扫频干扰的频率调制参数即初始频率fi和调频斜率gi,i为第i个线性扫频干扰信号;
寻找第i个线性扫频干扰对应时频脊线的邻域;在第i个线性扫频干扰信号对应的时频分布脊线上任取一个自源点,过该点作一条直线,使倾角α在0°~180°中变化,设定直线的邻域半径rline,所述邻域半径rline不超过干扰在时频分布图的时频脊线在任意时间点上频宽的一半,将时频分布图上各自源点到该直线的距离小于rline的自源点集合做为该直线的邻域Ωi,统计不同倾角下直线的邻域Ωi中自源点的数量,当某一倾角下直线的邻域Ωi中自源点的数量比该直线处于其它倾角时自源点的数量都大时,则将该直线倾角作为第i个线性扫频干扰信号在时频分布面上的倾角,根据该倾角计算出第i个线性扫频干扰信号的调频斜率;对于本实施例中的LFM干扰,邻域半径rline根据归一化的频率设定,干扰在时频分布图的时频脊线占有一定宽度,即在任意时间点上,某个干扰占有一段很窄的频谱,邻域半径rline不超过任意时间点上频宽的一半。对于本发明中的LFM干扰,时频脊线在任意时间点上频宽用归一化频率衡量不超过0.01,因此邻域半径rline用归一化频率衡量不超过频宽的一半即0.005,而0.005即0.5的1/100,在本研究中归一化的半频带宽0.5对应于数值1024,所以邻域半径rline对应于1024/100≈10。因此,对于LFM干扰的情况,取rline=10足够。需要指出的是,rline的设定不需要很严格,因为rline即使不取到理想值,对估计结果影响微弱。
步骤3.2:去除第i个线性扫频干扰对应时频脊线的邻域与其它线性扫频干扰对应时频脊线的邻域的交集,得到去除交集后的区域Ω′i,用最小二乘法拟合Ω′i邻域中的自源点,得到第i个干扰信号的拟合直线方程则为估计出的第i个线性扫频干扰信号的初始频率,为估计出的调频斜率。这里需要指出的是,在步骤2.1中,用旋转的方法可以粗略估计出LFM干扰的频率参数,而这一步在估计每个LFM干扰频率参数时,排除掉含有的别的LFM干扰自源点,所以更为精确。表1给出了场景1中按照步骤2所估计出的频率和调频斜率参数值。
表1 LFM干扰参数估计值
步骤4:对所估计的频率参数进行修正。
首先分析必须对频率参数估计值修正的必要性。
如上文所述,子空间基矢量通过截取某一个LFM干扰序列的一段构成,设对某一路接收信号x(t)截取长为L的数据序列,如下式
X=[x(1) x(2) … x(L)]T=Xs+Xu+V (27)
式中V为白噪声矢量,Xs、Xu分别表示有用信号和干扰矢量。
不妨设接收信号中只含有一个非平稳干扰,在构建干扰子空间时,子空间基矢量的幅度影响很小,可以采用归一化的幅度,因此子空间基矢量可以只考虑相位,前文已假设接收信号矢量中,幅度归一化的单纯非平稳干扰部分为y(t),下文中y(t)表示同样的意义,只是假设只含一个非平稳干扰,可写成标量形式。
故该干扰信号数学表达式可设为
选取该干扰信号中L长度的序列构成一列矢量并进行归一化,设为Y,如下式
干扰子空间构建的必要条件是估计干扰参数,参数估计很可能存在微小的误差,参数估计误差导致的结果是根据估得参数构建的信号在任意时刻与真实的信号存在一个相位差,可设为据此估计出的干扰信号可表示为
用上式表示的非平稳信号构建干扰子空间,由于只有一个干扰信号,故干扰子空间只有一个基矢量,将干扰子空间基矢量写为
上式中矢量形式的Y′表示选取长度为L的干扰信号序列,排成矢量形式。得出干扰子空间基矢量后可将(10)式表示的干扰子空间正交矩阵写为
抑制接收信号中干扰的过程也就是将接收信号投影到干扰子空间的正交子空间,用数学形式表示如下
Y=PX=PXs+PY+PV (32)
式(11)没有列出干扰本身投影到干扰正交子空间的项,这是因为在干扰子空间估计准确的情况下,干扰投影到其正交子空间的结果为零,当干扰参数估计存在偏差的时候,这一项就不能忽略不计了。由于PXs和PV在前文已经进行分析,所以下面着重分析PY项,如下式
将(29)式代入(30)式可以得到干扰信号矢量的表达式
同理,存在误差时的干扰信号矢量可表示为
将(34)(35)两式代入Y′HY的计算,可得
将(36)式代入(33)式,有
对于式(37),当L取很大的值时,比如说GPS信号C/A码码长1023,则有
当L取值较小时,因为对于i=1,2,…,L,相位差均很微小,因此可以推出
由此可以推出
据此可以得出结论:1、当按照文献2中的方法以GPS信号C/A码码长1023为长度构建子空间的基矢量时,对于干扰频率参数估计的微小误差将导致干扰不能被抑制;2、基矢量长度取值较小,虽然干扰频率参数估计有微小误差时,干扰也能被消除,但是子空间基矢量长度过短相当于对信号进行加窗,窗口很窄,频谱泄露越严重,所以子空间基矢量长度不是越短越好。正由于以上原因,所以有必要对频率参数误差进行修正。
步骤4.1、求取根据估计的初始频率参数构建的LFM干扰信号和接收信号的相关值。
所构建的信号与接收信号相关度越大,就说明用于构建信号的参数估计得越准确。所需要估计的频率参数有两个,所以这个问题就转化为了在以两个频率变量为横纵坐标的平面上寻找最大值的问题。显然,在平面上寻找最大值需要一个起点,这个起点在平面上的坐标即是前一步估计到的初始频率和调频斜率,这两个值都是估计值,记为和计算起点位置的相关值,即求估计参数构建的LFM干扰信号与接收信号的某一路x(t)的内积,即
式中T为求内积的时间长度。写成离散形式,即
至此就求得了初始时刻构建信号和接收信号的相关值。
步骤4.2、以估计值为中心,选取周边点。
设定初始频率和调频斜率的步进值δf和δg,在以初始频率f为横轴,调频斜率g为纵轴的平面上,以为中心,选取其左右上下四个点,即
步骤4.3、求取各个点对应的内积
将周围的四个点对应的频率参数值代入(1.2)式当中,计算相应的内积值
步骤4.4、更新中心点
求取各个点对应的内积,比较中心点和周围四个点对应内积的大小,用五个点中对应内积最大的点取代中心点,返回S2步,如果中心点对应的内积已经是最大,则停止运算,输出该内积最大值所对应的点坐标,该坐标为修正后的初始频率和调频斜率。
步骤5:用修正后的频率参数构建干扰子空间并用子空间投影消除干扰。
步骤5.1、根据修正后的LFM干扰参数构建干扰信号序列。
用修正后的频率参数构建的信号为 取xi(n)序列的L长度取一段数据作为子空间基矢量,如下
步骤5.2、构建干扰子空间及其正交子空间。
用时域扩展后的干扰信号矢量构建成一个矩阵,如下
B=[X1 X2 … XP] (8)
该组矢量张成了干扰的子空间,其正交子空间表示如下
步骤5.3、进行子空间投影消除干扰。
将接收信号矢量投影到干扰子空间的正交子空间,如下式
Y=PX=PXs+PV (10)
式中X为接收信号矢量,通过投影过程消除了接收信号中的干扰分量,剩余有用信号和白噪声分量Xs和V。
本发明同时还给出了一种基于子空间误差修正的线性扫频干扰抑制系统,包括处理器和存储器,存储器中存储有基于子空间误差修正的线性扫频干扰抑制方法的程序,处理器在运行所述基于子空间误差修正的线性扫频干扰抑制方法的程序时执行本发明所提供的方法各步骤。
下面通过两个实验来说明本发明的关键步骤的执行和所取得的有益效果。这两个实验分别为:1、频率参数修正;2、通过频率参数修正以修正干扰子空间这种方法的抗干扰效果。
1、频率参数修正
使用场景一设置,以第一个LFM干扰为例,其数学表达式为x1(t)=exp(j(2π*0.1t+0.1/2t2)),即初始频率f1=0.1,调频斜率g1=0.1。估计得到的在(f1,g1)邻域内,根据以上分析,当估计得到的初始频率和调频斜率与其真实值一致时,接收信号序列与用估计参数得到的信号序列做内积,获得的内积值最大。现以(f1,g1)为中心,取其邻域,绘出邻域内的点对应的内积值如附图3所示。
附图3(a)用等高线表示(f1,g1)邻域内对应内积值的大小;附图3(b)则用将内积值表现在Z轴上。由该图可以看出:1、当估计值等于真实值(f1,g1)时,内积值最大;2、同等调频斜率误差对内积值的影响大于初始频率。由此可以推出,当估计了LFM干扰的频率参数之后,用本发明方法逐步寻找内积最大值对应的点,也就可以准确估计出干扰信号的初始频率和调频斜率,即完成了对初始估计初始频率和调频斜率的修正。
2、本发明方法的抗干扰效果
仍然采用场景一,现对照三种情况下的抗干扰效果,这三种情况为:一、用估计得到频率参数构建干扰子空间;二、对估计的频率参数进行修正,再构建干扰子空间;三、使用精确频率参数构建干扰子空间。第一种情况即用步骤2方法估计得到频率参数后,即构建干扰子空间,然后进行投影。第二种情况即在步骤2的基础上用步骤3方法对频率参数进行修正,用修正后的频率参数构建干扰子空间,然后进行子空间投影抗干扰。第三种情况是使用精确频率参数构建子空间进行抗干扰。列举第三种情况是因为使用本发明方法也不一定能够将频率参数估计得十分精确。
由文献2可知,当构建子空间的参数完全准确时,即使在构建子空间的序列长度很大情况下(如取L=1023),子空间投影的效果都很理想。但在实际环境中,即使用本发明方法对估计的参数进行修正,也不一定能得到完全精确的值,但能够确定的是,误差至少减少了一个数量级。现假设干扰参数经修正后还存在极小的误差,经多次仿真得知,参数最多也就差一个步进值的误差。经多次仿真,参数修正准确的情况占绝大多数,但应该考虑较苛刻的情况,即不妨对x3(t)中的初始频率引入一个步进值的误差δf=0.0001。抗干扰仿真结果如附图4。
如附图4,图中实线、点划线、星划线分别对应于频率参数取精确值、步骤2得出的估计值、经步骤3得到的修正值三种情况。图中横坐标为构建子空间的序列长度,纵坐标为输出信干噪比(SINR)。由图可知:1、当用于构建干扰子空间的频率参数取精确值时,序列长度增加对于输出SINR没有影响,始终保持很高的输出SINR;2、步骤2方法估计干扰的频率参数,并据此构建干扰子空间,在此情况下,当用于构建干扰子空间序列较短时,能够获得较高的输出SINR,随着序列长度增加,SINR下降很快,一般序列长度增加到32,使得SINR下降到-25dB,有用信号就难以从噪声中恢复了;3、在步骤2基础上,通过步骤3修正,能够将误差较大程度地缩小,甚至完全消除,此处假设频率参数仍有微小误差,在此情况下,输出SINR也随着序列长度增加下降,但是相比情况2,SINR下降要缓慢得多,直到子空间序列长度取120,仍能保持较高的输出SINR,约为-21dB。
至此可以得出结论:1、用本发明可以有效减小干扰频率估计误差,从而很大程度上提高输出SINR;2、尽管构建子空间序列长度取短可以提高输出SINR,但是截取序列相当于加窗,窗口越短,频谱泄露越严重;3、能够通过缩减频率误差,从而避免使用短序列构建子空间,是相对更优化的方式。
Claims (7)
1.一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:包括以下步骤:
步骤1:对用阵列天线接收到的信号白化处理并进行时频分析得到优化后的接收信号时频分布图;
步骤2:从优化后的时频分布图中遴选自源点;
步骤3:通过自源点估计线性扫频频率参数即初始频率fi和调频斜率gi,i为第i个线性扫频干扰信号,i=2,…,P,P为干扰信号的数量;
步骤4:对所估计的频率参数进行修正;
步骤5:使用修正后的频率参数构建干扰子空间,用子空间投影消除干扰。
2.根据权利要求1所述的一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:步骤1得到时频分布图的具体方法为:
步骤1.1:求取接收信号矢量的白化矩阵W;
步骤1.2:用白化矩阵将接收信号矢量白化并用Wigner-Ville变换求取时频分布,得到优化后的接收信号时频分布图。
3.根据权利要求2所述的一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:所述步骤2中遴选自源点的方法为:
步骤2.1:对白化后的接收信号矢量z(n)求取空时频分布矩阵序列,设时频分析中所取的数据段长度为L,如下式
n表示时间点,ω表示归一化频率,l表示偏移量,Dzz(n,ω)表示白化后的接收信号矢量z(n)的空时频分布矩阵序列。
步骤2.2:求得时频分布图上所有时频点对应的空时频分布矩阵的迹并归一化,即Trace(E{Dzz(n,ω)})/L,设ζ=Trace(E{Dzz(n,ω)})/L;
步骤2.3:设置门限值ε,上述空时频分布矩阵的迹均已归一化,其值域为(0,1),0.1≤ε≤0.3,当空时频分布矩阵的迹ζ>ε时,则此矩阵对应的时频点(n,ω)属于自源点。
4.根据权利要求3所述的一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:步骤3通过自源点估计线性扫频频率调制参数初始频率fi和调频斜率gi的具体方法为:
步骤3.1:寻找第i个线性扫频干扰对应时频脊线的邻域;在第i个线性扫频干扰信号对应的时频分布脊线上任取一个自源点,过该点作一条直线,使倾角α在0°~180°中变化,设定直线的邻域半径rline,所述邻域半径rline不超过干扰在时频分布图的时频脊线在任意时间点上频宽的一半,将时频分布图上各自源点到该直线的距离小于rline的自源点集合做为该直线的邻域Ωi,统计不同倾角下直线的邻域Ωi中自源点的数量,当某一倾角下直线的邻域Ωi中自源点的数量比该直线处于其它倾角时自源点的数量都大时,则将该直线倾角作为第i个线性扫频干扰信号在时频分布面上的倾角,根据该倾角计算出第i个线性扫频干扰信号的调频斜率;
步骤3.2:去除第i个线性扫频干扰对应时频脊线的邻域与其它线性扫频干扰对应时频脊线的邻域的交集,得到去除交集后的区域Ωi′,用最小二乘法拟合Ωi′邻域中的自源点,得到第i个干扰信号的拟合直线方程则为估计出的第i个线性扫频干扰信号的初始频率,为估计出的调频斜率。
5.根据权利要求4所述的一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:所述步骤4中对所估计的频率调制参数进行修正的具体方法为:
步骤4.1:求取根据所估计的频率参数构建的LFM干扰信号和接收信号的相关值;
所述相关值是指:对用估计参数构建的LFM干扰信号与接收信号的某一路x(t)求内积,
式中T为求内积的时间长度,其离散形式为:
L为基矢量序列长度;
步骤4.2:以估计值为中心,选取周边点;
设定初始频率和调频斜率的步进值δf和δg,在以初始频率f为横轴,调频斜率g为纵轴的平面上,以为中心,选取其左右上下四个点,即
步骤4.3:求取各个点对应的内积
将周围的四个点对应的频率参数值代入(1.2)式当中,计算相应的内积值
步骤4.4:更新中心点
比较中心点和周围四个点对应内积的大小,用五个点中对应内积最大的点取代中心点,返回步骤4.2,如果中心点对应的内积已经是最大,则停止运算,输出该内积最大值所对应的点坐标,该坐标为修正后的初始频率和调频斜率。
6.根据权利要求5所述的一种基于子空间误差修正的线性扫频干扰抑制方法,其特征在于:步骤5中所述构建干扰子空间,用子空间投影消除干扰的具体方法为:
步骤5.1:根据修正后的频率参数构建干扰信号序列;
用修正后的频率参数构建的信号为 取xi(n)序列的L长度数据作为子空间基矢量,如下
步骤5.2:构建干扰子空间及其正交子空间
用时域扩展后的干扰信号矢量构建成一个矩阵,如下
B=[X1 X2 … XP] (8)
该组矢量张成了干扰的子空间,其正交子空间表示如下
步骤5.3:进行子空间投影消除干扰;
将接收信号矢量投影到干扰子空间的正交子空间,如下式
Y=PX=PXs+PV (10)
式中X为接收信号矢量,通过投影过程消除了接收信号中的干扰分量,剩余有用信号和白噪声分量Xs和V。
7.一种基于子空间误差修正的线性扫频干扰抑制系统,其特征在于:包括处理器和存储器,所述存储器中存储有基于子空间误差修正的线性扫频干扰抑制方法的程序,所述处理器在运行所述基于子空间误差修正的线性扫频干扰抑制方法的程序时执行1-6任一项所述方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811540877.4A CN109581422B (zh) | 2018-12-17 | 2018-12-17 | 一种基于子空间误差修正的线性扫频干扰抑制方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811540877.4A CN109581422B (zh) | 2018-12-17 | 2018-12-17 | 一种基于子空间误差修正的线性扫频干扰抑制方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109581422A true CN109581422A (zh) | 2019-04-05 |
CN109581422B CN109581422B (zh) | 2020-07-17 |
Family
ID=65929734
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811540877.4A Active CN109581422B (zh) | 2018-12-17 | 2018-12-17 | 一种基于子空间误差修正的线性扫频干扰抑制方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109581422B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113419228A (zh) * | 2021-06-02 | 2021-09-21 | 中国人民解放军海军航空大学航空作战勤务学院 | 基于时频脊-Radon变换的海面小目标检测方法和装置 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR20030094545A (ko) * | 2002-06-04 | 2003-12-18 | 한국항공우주연구원 | 디지피에스 보정신호에 선형주파수변조를 이용하여지피에스와 독립된 항법해 및 측위를 구하는 방법 |
CN102999473A (zh) * | 2012-10-18 | 2013-03-27 | 中国人民解放军电子工程学院 | 一种线性调频信号的检测与参数估计方法 |
CN104698475A (zh) * | 2015-04-02 | 2015-06-10 | 芜湖航飞科技股份有限公司 | 一种卫星导航接收机仿真抗干扰测试方法 |
CN106019237A (zh) * | 2016-06-23 | 2016-10-12 | 哈尔滨工业大学(威海) | 雷达lfm复合波形设计方法 |
CN108107452A (zh) * | 2018-01-16 | 2018-06-01 | 北京理工大学 | 一种基于自适应分段子空间投影的干扰抑制方法 |
-
2018
- 2018-12-17 CN CN201811540877.4A patent/CN109581422B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR20030094545A (ko) * | 2002-06-04 | 2003-12-18 | 한국항공우주연구원 | 디지피에스 보정신호에 선형주파수변조를 이용하여지피에스와 독립된 항법해 및 측위를 구하는 방법 |
CN102999473A (zh) * | 2012-10-18 | 2013-03-27 | 中国人民解放军电子工程学院 | 一种线性调频信号的检测与参数估计方法 |
CN104698475A (zh) * | 2015-04-02 | 2015-06-10 | 芜湖航飞科技股份有限公司 | 一种卫星导航接收机仿真抗干扰测试方法 |
CN106019237A (zh) * | 2016-06-23 | 2016-10-12 | 哈尔滨工业大学(威海) | 雷达lfm复合波形设计方法 |
CN108107452A (zh) * | 2018-01-16 | 2018-06-01 | 北京理工大学 | 一种基于自适应分段子空间投影的干扰抑制方法 |
Non-Patent Citations (2)
Title |
---|
Z. ZHU 等: ""Research on interference suppression of GPS receiver"", 《NATIONAL UNIVERSITY OF DEFENSE TECHNOLOGY》 * |
周柱: ""GPS接收机抗干扰研究"", 《中国博士学位论文全文数据库 基础科学辑》 * |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113419228A (zh) * | 2021-06-02 | 2021-09-21 | 中国人民解放军海军航空大学航空作战勤务学院 | 基于时频脊-Radon变换的海面小目标检测方法和装置 |
Also Published As
Publication number | Publication date |
---|---|
CN109581422B (zh) | 2020-07-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104537249B (zh) | 基于稀疏贝叶斯学习的波达方向角估计方法 | |
CN109752710B (zh) | 一种基于稀疏贝叶斯学习的快速目标角度估计方法 | |
CN107870314B (zh) | 基于极化敏感阵列的完备电磁分量加权融合测向优化方法 | |
CN107290742B (zh) | 一种非线性目标跟踪系统中平方根容积卡尔曼滤波方法 | |
CN109471082A (zh) | 基于信号子空间重构的阵元缺损mimo雷达角度估计方法 | |
CN105445718B (zh) | 一种基于阵列重构的分布式多载舰超视距雷达的doa估计方法 | |
CN109116293A (zh) | 一种基于离格稀疏贝叶斯的波达方向估计方法 | |
CN110196410A (zh) | 一种阵列天线主瓣干扰抑制方法及系统 | |
Xiaoli et al. | Weighted least squares state estimation based on the optimal weight | |
CN109597046A (zh) | 基于一维卷积神经网络的米波雷达doa估计方法 | |
CN110231620B (zh) | 一种噪声相关系统跟踪滤波方法 | |
CN109507635A (zh) | 利用两个未知方位辅助源的阵列幅相误差估算方法 | |
CN108535698A (zh) | 基于波束空间的米波雷达低仰角估计方法 | |
CN107229032A (zh) | 一种构建四阵元立体阵列的方法和装置 | |
CN103616661A (zh) | 一种稳健的远场窄带信号源个数估计方法 | |
CN110161489A (zh) | 一种基于伪框架的强弱信号测向方法 | |
CN104375121A (zh) | 基于先验信息的mimo雷达波形与有偏估计器的联合优化方法 | |
CN104215939A (zh) | 一种融合广义对称结构信息的知识辅助空时自适应处理方法 | |
CN113534199A (zh) | 一种自适应的广义累计和gps欺骗攻击检测方法 | |
CN109541306A (zh) | 一种基于tls-esprit的谐波间谐波检测方法 | |
CN110244260B (zh) | 基于声能流矢量补偿的水下目标高精度doa估计方法 | |
CN109581422A (zh) | 一种基于子空间误差修正的线性扫频干扰抑制方法及系统 | |
CN105572629B (zh) | 一种适用于任意阵列结构的低运算复杂度的二维测向方法 | |
CN104298850A (zh) | 信源数未知的相干信号测向方法及系统 | |
CN105652256B (zh) | 一种基于极化信息的高频地波雷达tbd方法 |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |