CN109557513A - 一种动态环境下非平稳干扰的抑制方法及系统 - Google Patents

Abstract

本发明公开了一种动态环境下非平稳干扰的抑制方法，包括以下步骤；1、对用阵列天线接收到的信号进行时频分析，估计干扰的波达方向DOA；2、在动态环境下抑制干扰。本发明同时公开了一种动态环境下非平稳干扰的抑制系统，包括以下模块：干扰波达方向估计模块；时频分析子模块，自源点遴选子模块，动态环境干扰抑制模块；本发明可以清晰的从信号的时频分布图中选出干扰信号的自源点，用自源点对应的时频分布矩阵作干扰信号DOA估计，依据DOA估计结果进行零陷扩展，可以实现消除动态环境下的非平稳干扰的效果。

Description

一种动态环境下非平稳干扰的抑制方法及系统

技术领域

本发明主要涉及到信号处理领域，特指一种动态环境下非平稳干扰的抑制方法及系统。

背景技术

战场环境下，接收机容易受到非平稳干扰和载体振动的影响。对于非平稳干扰，使用处理平稳干扰的方法处理效果不好，文献1“D C Ricks,P G Cifuentes,J SGoldstein.What is Optimal Processing for Nonstationary Data[C].Pacific Grove,CA:34th Asilomar Conference on Signals,Systems,and Computers,2000,11,656-661”对此进行了分析和仿真。接收机载体振动引起的结果是：在权值更新的时间间隔内，干扰入射方向容易移出原零陷位置，致使抗干扰处理之后的干扰残余较大。

运用现有技术处理以上问题主要步骤可归纳为两大步骤：一、对接收信号进行时频分析，在时频分布图上找出干扰信号对应的时频点，称之为自源点，算出自源点对应的空时频分布(STFD:Spatial Temporal Frequency Distribution)矩阵序列，该矩阵见于文献2“Adel Belouchrani,Moeness G.Amin.Blind Source Separation Based on Time–Frequency Signal Representations[J].IEEE Transactions on Signal Processing,1998,46(11):2888-2897”之中，用空间谱方法如MUSIC算法分析该矩阵序列，即可估计出干扰的波达方向(DOA:Direction Of Arrival)；二、在获取干扰信号DOA的条件下，可在干扰方向形成零陷并且扩展零陷宽度，从而应对载体的振动，根据最小方差准则，可用文献3“Meng Hua Er.Linear Antenna Array Pattern Synthesis with Prescribed BroadNulls[J].IEEE Transactions on Antennas and Propagation,1990,38(9):1496-1498.”的方法将零陷展宽下抗干扰转化为基于接收信号矢量协方差矩阵的二次约束问题。通过以上处理即可在动态环境下消除非平稳干扰。具体步骤为：

一、干扰方向估计

现以GPS信号接收为例进行分析，采用何种导航信号对于抗干扰而言，处理方法相同。接收机接收到的信号由信号、白噪声和干扰组成。对于用阵列天线接收到的信号，其数学形式可以表示为

其中s(t)，S表示GPS信号的天线阵响应；y_i(t)，a_i分别表示第i个干扰及其天线阵列响应。设有P个干扰从不同方向入射。干扰带宽相对于载频来说可以看作窄带信号，不同阵元接收到信号的相位差仅仅取决于信号入射方向，认为不同阵元间接收信号的时差内信号频率不会发生变化。v(t)是均值为零的加性高斯白噪声。

为方便分析，按照文献2的方法，将干扰幅度提取到空间加权矢量之前，使干扰能量归一化，如下式所示

干扰带宽相对于载频可以看作窄带。y(t)表示多个干扰将振幅归一化后的合集，记为y(t)＝[y₁(t),y₂(t),…,y_P(t)]^T,各干扰之间不相关，其中y_i(t)的振幅为b_i。A由各干扰的导向矢量与对应干扰振幅的乘积组成，可写为A＝[b₁a(θ₁),b₂a(θ₂),…,b_Pa(θ_P)]^T,其中a(θ_i)为导向矢量，表示第i个信号在各个阵元上相对于参考阵元的相位差，θ_i为该信号的到达角。GPS有用信号很微弱，在对非平稳干扰的分析过程中不产生影响，可以忽略不计。

对于接收信号用Wigner-Ville变换进行时频分析可以得到良好的时频聚集性，信号的Wigner-Ville变换生成的时频分布称为Wigner-Ville分布(WVD)，WVD具有良好的局部性质，是一种最基本、也是应用最多的时频分布，取接收信号矢量的任意一路，设为x(t)，其Wigner-Ville变换如下

WVD具有理想的时频聚集性，但有多分量时，交叉项严重，交叉项在时频分布图中的体现是大量的互源点。

对于矢量形式的接收信号x(t)，矢量各分量的自Wigner-Ville变换和分量之间的互Wigner-Ville变换形成了一个M×M维的矩阵，即STFD矩阵，信号的时频分布图中每一个时频点均对应一个STFD矩阵。根据文献4“Adel Belouchrani,Moeness G.Amin.Time-frequency MUSIC[J].IEEE Signal Processing Letters,1999,6:109-110.”，STFD矩阵的离散表示式为

假设采样间隔为T_s，式中n代表时间nT_s，ω为归一化频率，l₁、l₂表示偏移量，用于度量x(n)在不同位置时的相关度。

按照(4)式求取接收信号矢量对应的STFD矩阵，(4)式中所取的是无限长度，实际只需要取有限长序列，如下式所示。

由上式可知，STFD矩阵是一个矩阵序列，每一个时间和频率点对(n,ω)对应一个M×M矩阵，其中M为矢量x(n)的维度。L为所取数据段长度，需要指出的是，对上式而言，如果L为偶数，则求和上下限不再是整数，在这种情况下将积分上下限取整即可，积分限有一位的误差对结果影响微乎其微。

如上文所述，信号的Wigner-Ville变换时频聚集性良好，但是交叉项严重，难以准确遴选出干扰信号对应的时频点，故而难以准确选出相应的STFD矩阵作空间谱估计。对于这一问题，文献2通过将接收信号矢量白化，达到了在时频分布图中消除互源点的目的。不妨设用文献2中的方法推算得到的白化矩阵为W，然后用W将接收信号矢量白化，将白化后的接收信号矢量记为z(n)，则有

z(n)＝Wx(n)＝W(Ay(n)+v(n))＝Φy(n)+Wv(n) (6)

将接收信号矢量白化之后，再用Wigner-Ville变换求取时频分布，可以得到优化后时频分布图。

文献2在获取优化后的时频分布图之后，没有给出在时频分布图中遴选自源点的方法，且文献2在分析中忽略了白噪声项，不能很精准的选出自源点。文献5“Yimin Zhang,Moeness G.Amin.Blind Separation of Nonstationary Sources Based on SpatialTime-Frequency Distributions[J].EURASIP Journal on Applied Signal Processing,Vol 2006,P 1–13.”给出一种遴选自源点的方法，即以下式作为门限进行选择。

上式中，η为判决变量，Trace表示求矩阵迹的运算，Norm表示求矩阵范数的运算，一般采用2-范数，即矩阵最大特征值开方。D_zz(n,ω)即对白化后接收信号矢量z(n)进行空时频分析得到的STFD矩阵序列，每一个时频点(n,ω)对应一个STFD矩阵。设定一个门限ζ，只有当η≥ζ时，才认为D_zz(n,ω)对应的时频点为自源点。

文献5中将矩阵的迹除以矩阵范数的目的是归一化，这种方法存在一种问题。通过观察(7)式可知，不仅矩阵的迹随时频点(n,ω)变化，矩阵的范数也随(n,ω)变化，这就带来一个问题，当矩阵的迹大时，其范数也大，导致自源点无法凸显，所以用(7)式作为判决变量是不合适的。而且该方法忽视了一个重要因素，也就是既不是自源点也非互源点的时频点，其对应矩阵的迹和矩阵范数均很小，其比值不一定是一个接近于零的小数，因此用该方法不能可靠的从所有时频点中选出自源点。

不妨设用现有技术已经在时频分布图中找出干扰对应的自源点以及相应的STFD矩阵。运用这些矩阵进行联合对角化即可估计出干扰方向。运用文献3“Adel Belouchrani,Moeness G.Amin,Karim Abed-Meraim.Direction Finding in Correlated Noise FieldsBased on Joint Block-Diagonalization of Spatio-Temporal Correlation Matrices[J].IEEE Signal Processing Letters,1997,4(9):266-268”中联合对角法综合自源点对应的时频分布矩阵，从而估计非平稳干扰DOA。

将各自源点对应的STFD矩阵进行平均后所得的矩阵记为这里需要指出的是：因自源点已经遴选得到，就没有必要再用白化后的矢量生成的STFD矩阵，只需要用原接收信号矢量x(n)算得的STFD矩阵即可。现对STFD矩阵进行特征分解，如下

矩阵为M维矩阵，对其进行特征值分解得到M个特征值λ₁,λ₂,…,λ_M，及对应的特征向量t₁,t₂,…,t_M。运用最小描述长度(MDL)准则可以估计出信号子空间维数，应为P，则利用MUSIC算法构造空间谱估计如下

式中a(θ)＝exp(-j2πfdcosθ/c),f为入射信号频率，d为阵元间距，θ为入射信号与天线阵夹角，c为光速。将信号入射角从0度到180度变化，所得空间谱的谱峰即是信号入射方向。

二、零陷扩展

在估计出干扰方向之后，通过对干扰方向施加约束，即可抑制干扰。如考虑载体快速运动，则需要扩展静态或低速运动情况下天线方向图对准干扰的零陷。

首先求得接收信号矢量的协方差矩阵，如下

式中省去了有用信号项，θ_i、a_i(θ_i)、r_i分别为第i个干扰的入射方向、导向矢量、干扰功率，为接收机噪声方差，I为单位矩阵。

然后求取不考虑零陷扩展时的最优权值，记为w₀，在MVDR准则约束下的最优权值，如下式

式中h为有用信号的在阵列天线上的导向矢量，如采用MVDR准则的特殊情况简单约束，则h＝[1,0,L,0]^T，维度为M。

考虑干扰入射角变动，设某干扰入射角为零陷加宽的预期目的是使得方向图中零陷扩展到邻近区域，设引入零限扩展之后的权值为w，则零陷展宽的数学意义为：式中为第i个干扰的导向矢量，ε是一个极小值。

根据最小方差准则，可按照文献3的方法将零陷展宽下抗干扰转化为基于接收信号矢量协方差矩阵的二次约束问题：

minf(w)＝||w-w₀||²s.t.w^HQw＝ε (12)

式中Q为针对干扰入射角变量积分求平均获得的干扰入射方向导向矢量协方差矩阵，其表达式为式中为在变化区间内的概率密度函数，ε为预期的零陷深度。

式(12)描述的属于凸集最优化问题，即Covex Optimization。这个优化问题可以运用CVX软件求解，也可以用拉格朗日乘数法转化为显式的权值表达式。

以上两步即是运用现有技术处理接收机在动态环境下受到非平稳干扰的基本方法。但是因为不能可靠的从所有时频点中遴选出自源点导致不能准确的估计出干扰信号的入射方向，也就不能很好的消除非平稳干扰。

发明内容

本发明要解决的技术问题是提供了一种可准确估计干扰信号入射方向，从而可有效消除非平稳干扰的动态环境下非平稳干扰的抑制方法。

为解决上述问题，本发明采取的技术方案是：

一种动态环境下非平稳干扰的抑制方法，包括以下步骤：

步骤1：对用阵列天线接收到的信号进行时频分析，估计干扰的波达方向DOA；

步骤1.1：求取接收信号矢量的白化矩阵W；具体方法是：

步骤1.1.1求取接收信号矢量的自相关矩阵R_x，R_x＝E{x(n)x^H(n)}，其中，x(n)表示在时间点n时刻的接收信号矢量；

步骤1.1.2对自相关矩阵进行特征分解，

R_x＝ΦΛΦ^H＝Φdiag[λ₁,λ₂,…,λ_P]Φ^H

其中，Λ是由R_x的前P阶特征值组成的对角矩阵，即Λ＝diag[λ₁,λ₂,…,λ_P]，λ₁,λ₂,…,λ_P为自相关矩阵的前P阶特征值，Φ为特征值所对应的特征矢量组成的矩阵，

步骤1.1.3：根据特征分解所得到的特征值和特征向量求取接收信号矢量的白化矩阵；

式中为接收机噪声方差；

步骤1.2：用白化矩阵将接收信号矢量白化，z(n)＝Wx(n)，z(n)为白化后的接收信号矢量，用Wigner-Ville变换求取时频分布，得到优化后的时频分布图；

步骤1.3：从优化后的接收信号时频分布图中遴选自源点，具体为：

步骤1.3.1求取白化后接收信号矢量的STFD矩阵序列；

n为时间点，ω为归一化频率，l表示偏移量，L表示时频分析中所取的数据段长度；

步骤1.3.2在每一个时频点(n,ω)，计算在该时频点的STFD矩阵的迹，设为ζ，有

ζ＝Trace(E{D_zz(n,ω)})/L

其中，L为时频分析中所取的数据段长度；

步骤1.3.3设置门限ε_TH选出自源点；当时频点STFD矩阵的迹ζ＞ε_TH时，0.1≤ε_TH≤0.3，则此矩阵对应的时频点为自源点；

步骤1.4：估计干扰入射方向；

步骤1.4.1：将各自源点对应的STFD矩阵进行特征分解，获得特征向量，

矩阵为M维矩阵，对其进行特征值分解得到M个特征值λ₁，λ₂…λ_i…λ_M，及对应的特征向量t₁,t₂…t_i…t_M；

步骤1.4.2：利用MUSIC算法构造空间谱估计：

其中，a(θ)＝exp(-j2πfdcosθ/c),f为入射信号频率，d为阵元间距，θ为入射信号与天线阵夹角，c为光速；将信号入射角θ从0度到180度变化，所得空间谱的谱峰即是干扰信号入射方向；

步骤2：动态环境下抑制干扰；

步骤2.1：不考虑零陷扩展，求取静态或低速情况下针对接收信号矢量的协方差矩阵，

式中省去了有用信号项，x(n)为在时间点n时刻的接收信号矢量，θ_i、a(θ_i)、r_i分别为第i个干扰的入射方向、导向矢量、干扰功率，为接收机噪声方差，I为单位矩阵。；

步骤2.2：求取不考虑零陷扩展时协方差矩阵在MVDR准则约束下的最优权值

其中，h为有用信号的在阵列天线上的导向矢量，采用MVDR准则的特殊情况简单约束，则h＝[1,0,…,0]^T，维度为M；

步骤2.3：根据最小方差准则，将零陷展宽下抗干扰转化为基于接收信号矢量协方差矩阵的二次约束问题：

minf(w)＝||w-w₀||²s.t.w^HQw＝ε (12)

w为零限扩展之后的权值，为平均矩阵，表示对干扰入射角变量积分求平均获得的干扰入射方向导向矢量协方差矩阵,为在变化区间内的概率密度函数，为干扰入射角，为第i个干扰的导向矢量，ε为预期的零陷深度，是一个极小值；

步骤2.4：求解零陷加宽后的最优权值w；

步骤2.5：滤波，在获得最优权值w之后，对接收信号矢量进行矢量加权，即阵列滤波，实现在动态环境下消除干扰的目的。

一种动态环境下非平稳干扰的抑制系统，包括处理器，以及与所述处理器连接的存储器，所述存储器存储有动态环境下非平稳干扰的抑制程序，所述动态环境下非平稳干扰的抑制程序被所述处理器执行时实现上述方法的步骤，包括以下模块：

干扰波达方向估计模块；用于对用阵列天线接收到的信号进行时频分析，估计出干扰入射方向；包括以下子模块：

时频分析子模块：用于对阵列天线接收到的信号进行时频分析，求取接收信号矢量的白化矩阵，对接收信号矢量白化，得到优化后的时频分布图；

自源点遴选子模块，用于计算优化后的时频分布图中每一个时频点的STFD矩阵的迹，并判断每一个时频点的STFD矩阵的迹与设定的门限值的大小，获得自源点；

动态环境干扰抑制模块；用于根据遴选出的自源点估计出干扰入射方向，通过所估计出的干扰入射方向，扩展静态或低速运动情况下天线方向图对准干扰的零陷，求出零陷加宽后的最优权值，对接收信号矢量进行矢量加权，在动态环境下消除干扰。

与现有技术相比，本发明所取得的有益效果为：

本发明一种动态环境下非平稳干扰的抑制方法，在接收信号矢量白化后生成的优化时频分布图中，分析了每一时频点所对应STFD矩阵迹的取值范围，当该时频点所对应STFD矩阵迹大于设定的门限值时，则该时频点为自源点，本发明可以清晰的从信号的时频分布图中选出干扰信号的自源点，用自源点对应的时频分布矩阵作干扰信号DOA估计，依据DOA估计结果进行零陷扩展，可以实现消除动态环境下的非平稳干扰的效果。本发明一种动态环境下非平稳干扰的抑制系统，通过干扰波达方向估计模块中的自源点遴选子模块，清晰的遴选出自源点，用自源点对应的时频分布矩阵作干扰信号DOA估计

附图说明

图1为时频分布图的优化；

图2为本发明方法与现有方法在遴选时频点这一功能上的对比；

图3为用自源点对应STFD矩阵所作DOA估计；

图4为零陷扩展验证图；

图5为零陷扩展效果图。

具体实施方式

附图1至5示出了本发明一种动态环境下的非平稳干扰抑制方法的具体实施例。该方法包括以下步骤：

步骤1：对用阵列天线接收到的信号进行时频分析，估计干扰的波达方向DOA。

步骤1.1：求取接收信号矢量的白化矩阵W；具体方法是：

求取接收信号矢量的自相关矩阵，对自相关矩阵进行特征分解，根据特征分解所得到的特征值和特征向量求取接收信号矢量的白化矩阵；

按照公式(2)的设定，将所有干扰之和记作u(t)，采用离散化表示，记作u(n)，设干扰矢量u(n)对应白化矩阵为W，按照定义对白化过程进行分析，白化即是使得干扰相关矩阵转换为单位矩阵，将单位矩阵记为I，白化过程可分析如下

E{Wu(n)u^H(n)W^H}＝WR_uW^H＝WAE{y(n)y^H(n)}A^HW^H (13)

据前文所设，各干扰不相关，又干扰功率已经归一化，所以E{y(n)y^H(n)}＝I，可得

E{Wu(n)u^H(n)W^H}＝WAA^HW^H (14)

上式表明如果W为白化矩阵，则WA为一个酉矩阵，令该酉矩阵为Φ，则有

A＝W^#Φ (15)

式中W^#表示W的Moore-Penrose逆矩阵，白化过程中将阵列流形矩阵转化为了酉矩阵，矩阵维数从阵列的维数M降到了子空间维数P。

下面求解白化矩阵W，

步骤1.1.1对含有白噪声的接收信号求自相关矩阵R_x；

步骤1.1.2对自相关矩阵进行特征分解，

R_x＝ΦΛΦ^H＝Φdiag[λ₁,λ₂,…,λ_P]Φ^H

式中λ₁,λ₂,…,λ_P为矩阵特征值，不妨设干扰子空间的秩为P，对矩阵进行分解只取其P阶之前的主体部分，Φ为特征值所对应的特征矢量组成的矩阵，

根据(13)式，可知R_u＝AA^H，由此可以推出

式中为白噪声方差，将A＝W^#Φ代入(16)式，有

AA^H＝(W^#Φ)(W^#Φ)^H (18)

根据(16)(17)式，有

由此可以推导出W^#的表达式。设则有

令将(19)式简写为W^#＝ΦΛ′。

步骤1.1.3：根据特征分解所得到的特征值和特征向量求取接收信号矢量的白化矩阵；

根据Moore-Penrose逆矩阵的定义，白化矩阵W与其Moore-Penrose逆矩阵W^#存在关系：W^#WW^#＝W^#。将W^#＝ΦΛ′代入该关系有

(ΦΛ′)W(ΦΛ′)＝ΦΛ′ (21)

因为Φ为酉矩阵，上式可以转化为

W＝(Λ′)^-1Φ^H (22)

最终可化为

至此就求得了白化矩阵。

步骤1.2：用白化矩阵将接收信号矢量白化，使用Wigner-Ville变换求取时频分布图，得到优化后的时频分布图；

求得白化矩阵之后，将接收信号矢量白化，将白化后的接收信号矢量记为z(n)，则有

z(n)＝Wx(n)＝W(Ay(n)+v(n))＝Φy(n)+Wv(n) (24)

将接收信号矢量白化之后，再用Wigner-Ville变换求取时频分布，可以得到优化后的时频分布图。通过优化可以消除时频分布的交叉项，也即消除了时频分布图中的互源点，突出自源点。

步骤1.3：从优化后的接收信号空时频分布图中遴选自源点，具体为：

步骤1.3.1求取白化后接收信号矢量的STFD矩阵序列；

n为时间点，ω为归一化频率，l表示偏移量，l表示偏移量，L表示时频分析中所取的数据段长度；

步骤1.3.2，计算每一个时频点(n,ω)的STFD矩阵的迹；

由式(24)可知白化信号与接收信号矢量的STFD矩阵D_xx(n,ω)存在下面关系

D_zz(n,ω)＝WD_xx(n,ω)W^H (26)

再将(2)式代入D_xx(n,ω)，有

D_xx(n,ω)＝D_uu(n,ω)+D_uv(n,ω)+D_vu(n,ω)+D_vv(n,ω) (27)

因为噪声和信号不相关，所以D_uv(n,ω)＝D_vu(n,ω)＝0，因此有

其中，D_yy(n,ω)是各干扰的空时频分布矩阵，它的元素是干扰信号的自Wigner-Ville分布和互Wigner-Ville分布。

以上是接收信号矢量x(n)、经白化后的接收信号矢量z(n)、干扰之和矢量u(n)以及归一化后的干扰之和矢量y(n)所对应的STFD矩阵D_xx(n,ω)、D_zz(n,ω)、D_uu(n,ω)以及D_yy(n,ω)之间的关系。据此可以推导出D_zz(n,ω)的特点如下。

将白化信号矢量的表示式(24)代入(25)式，并引入(28)式，可以得到下式

将干扰信号矢量y(n)的STFD矩阵记为

将推得的白化矩阵表达式(23)代入(25)式，有

式中令式(31)为任意时频点上白化信号的STFD矩阵。

将自源点对应的STFD矩阵记为将互源点对应的STFD矩阵记为自源点对应STFD矩阵的迹可推算如下

式(32)中假设时频点处在第q个信号的自项上，因此除外，其余对角元素均为0。就是信号的Wigner变换，记为不妨设信号为其有限长度序列下的Wigner分布可表示如下

上式只有当

时不为0，此时由此可知，时频平面中任何非时频脊线范围之外的点其幅度理论值约为零。式(34)也决定了各信号时频分布脊线的方程。

设时频分析中所取的数据段长度为L，则只有当时间点n处于数据段的中央时，才可能向两端取到长度为L的数据段。容易推知，当n_i＝n±i时，在时频点满足(34)式的情况下当n_i＝n±n时，因此可以推出下式

按照同样的方法推算互源点对应STFD矩阵的迹，如下

协方差矩阵的前P个特征值远大于白噪声方差，因此有

由此可知，接收信号矢量经白化后，其STFD矩阵的迹对于信号自项来说约为一个0到1之间的数，而交叉项对应STFD矩阵的迹约等于0,因此在每一个时频点(n,ω)，该时频点的STFD矩阵的迹，设为ζ

ζ＝Trace(E{D_zz(n,ω)})/L (37)

其中，L为时频分析中所取的数据段长度。

步骤1.3.3：设置门限ε_TH选出自源点；当时频点STFD矩阵的迹ζ＞ε_TH时，0.1≤ε_TH≤0.3，则此矩阵对应的时频点为自源点；

当STFD矩阵的迹ζ＞ε_TH时，认为此矩阵对应的时频点属于干扰信号的自项。部分矩阵的迹小于ε的自项对应的时频点未选入对结果影响很小，可忽略不计，本实施例中取ε_TH＝0.1。至此就在时频分布图中选得了干扰信号对应的自源点。

由以上推导可知，在白化后的时频分布图中，只有时频点落在LFM干扰的时频脊线之上，其对应的STFD矩阵的迹的值才不为0，理论上恒等于L，即故此，将(37)式所示的ζ作为遴选自源点的门限是合理的。因为ζ的最大值显然为1，而最小值可以忽略不计，所以用ε_TH＝0.1进行区分既可以很好的遴选出自源点。

相比较而言，如果用文献5的门限表达式，即(7)式，所得到的判决门限在时频分布图上的分布如附图2(a)所示，非常杂乱，无法遴选出自源点。这是因为文献5没有对每个时频点对应的矩阵的迹进行分析。而文献2在讨论了白化之后，并未对如何选择自源点进行分析，且文献2在分析中未考虑加入白噪声。

步骤1.4：估计干扰入射方向；

将各自源点对应的STFD矩阵进行特征分解，获得特征向量，利用MUSIC算法构造空间谱，将信号入射角从0度到180度变化，所得空间谱的谱峰即是信号入射方向；

步骤2：动态环境下干扰抑制

步骤2.1：不考虑零陷扩展，求取静态或低速情况下针对接收信号矢量的协方差矩阵，

式中省去了有用信号项，x(n)为在时间点n时刻的接收信号矢量，θ_i、a(θ_i)、r_i分别为第i个干扰的入射方向、导向矢量、干扰功率，为接收机噪声方差，I为单位矩阵。；

步骤2.2：求取不考虑零陷扩展时协方差矩阵在MVDR准则约束下的最优权值

其中，h为有用信号的在阵列上的导向矢量，

步骤2.3：根据最小方差准则，将零陷展宽下抗干扰转化为基于接收信号矢量协方差矩阵的二次约束问题：

minf(w)＝||w-w₀||²s.t.w^HQw＝ε (12)

w为零限扩展之后的权值，为干扰入射角变量积分求平均获得的干扰入射方向导向矢量协方差矩阵,为在变化区间内的概率密度函数，为干扰入射角，为第i个干扰的导向矢量，ε为预期的零陷深度，是一个极小值；

步骤2.4：求解零陷加宽后的最优权值。

对于该凸集最优化问题，即Covex Optimization，运用CVX软件进行求解。

步骤2.5：滤波

在算得权值w之后，用来对接收信号矢量进行矢量加权，即阵列滤波，即可达到在在动态环境下消除干扰的目的。

步骤2.3中，平均矩阵Q为针对干扰入射角变量积分求平均获得的干扰入射方向导向矢量协方差矩阵，其表达式为式中为在变化区间内的概率密度函数，ε为预期的零陷深度。

针对干扰入射角变化，平均矩阵Q用积分的方式求解，如公式38所示

为第i个干扰的导向矢量，如下式

将(39)式代入(38)式，可算得矩阵Q的第(m,n)个元素为

以上积分形式难以求解，因而需要改变积分变量。可将积分变量改为因为在[θ_i-Δθ_i,θ_i+Δθ_i]内均匀分布，而Δθ_i是一个很小的值，因而在角度变化范围内，也可以看作服从均匀分布，令则y的概率密度函数为对积分限进行分析有

sin(θ_i+Δθ_i)＝sinθ_icosΔθ_i+cosθ_isinΔθ_i≈sinθ_i+Δθ_icosθ_i (41)

同理，sin(θ_i-Δθ_i)≈sinθ_i-Δθ_icosθ_i，有

令Q_i(m,n)＝exp(-jπ(m-n)sinθ_i)，可以推出

式中“ο”表示Hardmard积，式中

步骤2.3：计算权值

采用拉格朗日乘子法求解式(12)描述的非线性约束的问题，如下

f(w)＝(w-w₀)^H(w-w₀)+α(w^HQw-ε) (44)

式中α为拉格朗日乘数因子，对f(w)关于w求导，并令结果为零，可得

w＝(I+αQ)^-1w₀ (45)

至此就获得了在估计出干扰信号来向，并对干扰来向进行零陷扩展情况下的最优权值。

文献6“邹翔,朱然刚,史英春,蒋云霄.一种基于二次约束的零陷加宽算法[J].电路与系统学报,2012,17(5):134-138.”限定了拉格朗日乘数因子α的取值范围，可参照选择。而w₀已在前面步骤算得，至此也就算得了零限扩展情况下的权值w。

一种动态环境下非平稳干扰的抑制系统，包括处理器，以及与所述处理器连接的存储器，所述存储器存储有动态环境下非平稳干扰的抑制程序，所述动态环境下非平稳干扰的抑制程序被所述处理器执行时实现上述权利要求1～2任一项所述方法的步骤，包括以下模块：

干扰波达方向估计模块；用于对用阵列天线接收到的信号进行时频分析，估计出干扰入射方向；包括以下子模块：

时频分析子模块：用于对阵列天线接收到的信号进行时频分析，求取接收信号矢量的白化矩阵，对接收信号矢量白化，得到优化后的时频分布图；

自源点遴选子模块，用于计算优化后的时频分布图中每一个时频点的STFD矩阵的迹，并判断每一个时频点的STFD矩阵的迹与设定的门限值的大小，获得自源点；

动态环境干扰抑制模块；用于根据遴选出的自源点估计出干扰入射方向，通过所估计出的干扰入射方向，扩展静态或低速运动情况下天线方向图对准干扰的零陷，求出零陷加宽后的最优权值，对接收信号矢量进行矢量加权，在动态环境下消除干扰。

下面通过实验对本发明所提出的系统和方法进行验证。

现设置一个场景对接收信号矢量的白化以及在时频分布图中遴选自源点进行检验，不妨称此场景为场景一。场景一：接收信号中含有3个非平稳干扰，不妨设为LFM干扰信号，设为x_i(t)＝A_iexp(j(2πf_it+g_i/2t²)),i＝1,2,3,其中f₁＝0.1,g₁＝0.1,f₂＝0.3,g₂＝-0.15，f₃＝0.45,g₃＝-0.1。干扰入射方向分别为40°、80°、120°。

将接收信号矢量白化前和白化后进行Wigner变换得到的时频分布图进行对比，如附图1所示。由图1可知，当接收信号中含有多个干扰时，时频分布图的交叉项严重，白化处理可以消除时频变换后的交叉项。

再将文献5的提出的在时频分布图中选择自源点的方法与本发明方法进行对比，如附图2所示。在文献5与本发明选择自源点方法时，需要同时讨论图1(b)，现讨论如下。附图2(a)表示的是运用文献5方法将时频分布图上所有时频点对应的STFD矩阵之迹，用该矩阵范数进行归一化之后，得到的时频分布图。附图2(b)表示的是用本发明方法，针对图1(b)选得的干扰信号自源点散点图。即：接收信号矢量经过白化后，求取其时频分布，得到图1(b)所示的时频分布图，得到该图后，应遴选时频点，现有两种处理方式：第一种方式用文献5的方法处理，得到如附图2(a)所示的矩阵之迹归一化之后的时频分布图，由图可见很难在图中选得干扰信号的自源点；第二种方式是采用本发明方法选择自源点，如附图2(b)所示，可见用该方法可以清晰的选出干扰信号的自源点。

步骤1：根据干扰信号的自源点进行干扰信号来向估计。

将各自源点对应的STFD矩阵进行平均后所得的矩阵记为对其进行特征分解，如(8)式，对进行特征值分解得到M个特征值λ₁,λ₂,…,λ_M，及对应的特征向量t₁,t₂,…,t_M。运用最小描述长度(MDL)准则可以估计出信号子空间维数，应为P，则利用MUSIC算法构造空间谱估计如(9)式所示。在(9)式中，将信号入射角从0度到180度变化，所得空间谱的谱峰即是信号入射方向。

仍然使用场景一，用本发明方法将时频分布图优化之后，将自源点对应的STFD矩阵做平均，即可估计干扰信号来向。其结果如附图3所示。选用非平稳干扰在时频分布图上自源点对应的时频分布矩阵作DOA估计的结果。一共在相同条件下仿真5次，由图可知空间谱的谱峰很尖锐，DOA估计得也很准。

步骤2：动态环境下进行干扰抑制；

步骤2.1和2.2、求取静态情况下最优权值；

将不考虑零陷扩展时的最优权值记为w₀，在MVDR准则约束下的最优权值如式(11)所示。

步骤2.3动态情况下最优权值计算；

在估计出干扰方向之后，通过对干扰方向施加约束，即可抑制干扰。考虑载体振动，则需要扩展静态或低速运动情况下天线方向图对准干扰的零陷。“背景技术”中给出了现有的处理方法，即：将对于干扰来向邻域施加约束，然后用Covex Optimization(CVX)这一凸集最优化软件计算零陷扩展情况下的最优权值。用CVX软件处理虽然能够得到好的效果，但是算得的最优权值非显式表达，难于实现。因此，在本发明中，引入干扰入射方向情况下，显式表达最优权值，便于实现。

S5.1、计算平均矩阵

如上文所述，平均矩阵计算式为(42)式，即平均矩阵由R_i和T_i的积累加而成。R_i的计算式为T_i的计算式为其中坐标(m,n)表示矩阵T_i在该位置的元素。R_i和T_i的计算式中，θ_i为第i个干扰的入射角，通过前面步骤已经估计得到，Δθ_i为该干扰入射角可能变动的最大范围，可以根据经验设置，一般在1.5°之内。

步骤2.4、求解零陷加宽后的最优权值

零陷扩展后的权值表达式采用采用拉格朗日乘子法求取，求得权值如式(44)所示，即w＝(I+αQ)^-1w₀，式中α为拉格朗日乘数因子，文献6限定了拉格朗日乘数因子α的取值范围，可参照选择。平均矩阵Q在S5.1步算得，而w₀已在前面步骤算得，至此也就获得了零限扩展情况下的权值表达式。

仍然使用场景一对算法进行验证。运用求得的权值可以得到相应的天线方向图，如附图4所示：其中左图为不进行零陷扩展时的方向图，也就是权值矢量w₀对应的方向图；右图为零陷扩展以后所得权值w对应的方向图。由图可知，零陷扩展之后，零陷由6个减少为3个，均对准干扰方向，且零陷加宽效果较好，能容忍干扰入射角在一定范围内的变动。需要指出的是，附图4中，零陷设置较宽，是为了使得视觉效果明显。

零陷扩展需要解决的是动态环境下，自适应权值更新时间间隔内，干扰方向移出原零陷的问题。现将零陷扩展之前与之后，随着干扰入射角偏移量增加，用两种方式获得的输出信干噪比。结果如附图5所示。

由附图5可以看出：如果采用传统方法，不进行零陷扩展，则随着干扰入射角偏移加大，输出SINR下降很迅速；如果采用本文方法处理，随着干扰入射角偏移量增加，输出SINR能够基本保持平稳，偏移量超过1.5°之后缓慢下降。

由以上所有步骤以及结果可以得出结论：使用本发明方法可以清晰的从信号的时频分布图中选出干扰信号的自源点，用自源点对应的时频分布矩阵作干扰信号DOA估计，可以得到准确结果，依据DOA估计结果，并用本发明方法进行零陷扩展，就可以在芯片上实现消除动态环境下的非平稳干扰的功能。

Claims (3)

1.一种动态环境下非平稳干扰的抑制方法，其特征在于：包括以下步骤：

步骤1：对用阵列天线接收到的信号白化处理并进行时频分析，估计干扰的波达方向DOA；

步骤1.1：求取接收信号矢量的白化矩阵W；具体方法是：

步骤1.1.1求取接收信号矢量的自相关矩阵R_x，R_x＝E{x(n)x^H(n)}，其中，x(n)表示在时间点n时刻的接收信号矢量；

步骤1.1.2对自相关矩阵进行特征分解；

R_x＝ΦΛΦ^H＝Φdiag[λ₁,λ₂,…,λ_P]Φ^H

其中，Λ是由R_x的前P阶特征值组成的对角矩阵，即Λ＝diag[λ₁,λ₂,…,λ_P]，λ₁,λ₂,…,λ_P为自相关矩阵的前P阶特征值，Φ为特征值所对应的特征矢量组成的矩阵，

步骤1.1.3：根据特征分解所得到的特征值和特征向量求取接收信号矢量的白化矩阵；

式中为接收机噪声方差；

步骤1.2：用白化矩阵将接收信号矢量白化，z(n)＝Wx(n)，z(n)为白化后的接收信号矢量，用Wigner-Ville变换求取时频分布，得到优化后的时频分布图；

步骤1.3：从优化后的接收信号时频分布图中遴选自源点，具体为：

步骤1.3.1求取白化后接收信号矢量的STFD矩阵序列；

n为时间点，ω为归一化频率，l表示偏移量，L表示时频分析中所取的数据段长度；

步骤1.3.2在每一个时频点(n,ω)，计算在该时频点的STFD矩阵的迹，设为ζ，有

ζ＝Trace(E{D_zz(n,ω)})/L

其中，L为时频分析中所取的数据段长度；

步骤1.3.3设置门限ε_TH选出自源点；当时频点STFD矩阵的迹ζ＞ε_TH时，0.1≤ε_TH≤0.3，则此矩阵对应的时频点为自源点；

步骤1.4：估计干扰入射方向；

步骤1.4.1：将各自源点对应的STFD矩阵进行特征分解，获得特征向量，

矩阵为M维矩阵，对其进行特征值分解得到M个特征值λ₁，λ₂…λ_i…λ_M，及对应的特征向量t₁,t₂…t_i…t_M；

步骤1.4.2：利用MUSIC算法构造空间谱估计：

其中，a(θ)＝exp(-j2πfdcosθ/c),f为入射信号频率，d为阵元间距，θ为入射信号与天线阵夹角，c为光速；将信号入射角θ从0度到180度变化，所得空间谱的谱峰即是干扰信号入射方向；

步骤2：动态环境下抑制干扰；

步骤2.1：不考虑零陷扩展，求取静态或低速情况下针对接收信号矢量的协方差矩阵，

式中省去了有用信号项，x(n)为在时间点n时刻的接收信号矢量，θ_i、a_i(θ_i)、r_i分别为第i个干扰的入射方向、导向矢量、干扰功率，为接收机噪声方差，I为单位矩阵；

步骤2.2：求取不考虑零陷扩展时协方差矩阵在MVDR准则约束下的最优权值

其中，h为有用信号的在阵列天线上的导向矢量，采用MVDR准则的特殊情况简单约束，则h＝[1,0,L,0]^T，维度为M；

步骤2.3：根据最小方差准则，将零陷展宽下抗干扰转化为基于接收信号矢量协方差矩阵的二次约束问题：

minf(w)＝||w-w₀||²s.t.w^HQw＝ε (12)

w为零限扩展之后的权值，为平均矩阵，表示对干扰入射角变量积分求平均获得的干扰入射方向导向矢量协方差矩阵,为在变化区间内的概率密度函数，为干扰入射角， 为第i个干扰的导向矢量，

ε为预期的零陷深度，是一个极小值；

步骤2.4：求解零陷加宽后的最优权值w；

步骤2.5：滤波，在获得最优权值w之后，对接收信号矢量进行矢量加权，即阵列滤波，在动态环境下消除干扰。

2.根据权利要求1所述的一种动态环境下非平稳干扰的抑制方法，其特征在于：步骤2.4中对公式12的求解方法为运用CVX软件或拉格朗日乘子法进行求解。

3.一种动态环境下非平稳干扰的抑制系统，其特征在于：包括处理器，以及与所述处理器连接的存储器，所述存储器存储有动态环境下非平稳干扰的抑制程序，所述动态环境下非平稳干扰的抑制程序被所述处理器执行时实现上述权利要求1～2任一项所述方法的步骤，包括以下模块：

干扰波达方向估计模块；用于对用阵列天线接收到的信号进行时频分析，估计出干扰入射方向；包括以下子模块：

时频分析子模块：用于对阵列天线接收到的信号进行时频分析，求取接收信号矢量的白化矩阵，对接收信号矢量白化，得到优化后的时频分布图；

自源点遴选子模块，用于计算优化后的时频分布图中每一个时频点的STFD矩阵的迹，并判断每一个时频点的STFD矩阵的迹与设定的门限值的大小，获得自源点；

动态环境干扰抑制模块；用于根据遴选出的自源点估计出干扰入射方向，通过所估计出的干扰入射方向，扩展静态或低速运动情况下天线方向图对准干扰的零陷，求出零陷加宽后的最优权值，对接收信号矢量进行矢量加权，在动态环境下消除干扰。