CN109100749A - 基于mwf中运用滑窗判定的噪声子空间估计方法 - Google Patents

Abstract

本发明涉及基于MWF中运用滑窗判定的噪声子空间估计方法。本方法根据MWF的特性，结合小矩阵最大特征值运算，提出一种运用滑窗定阶的方法，可以鲁棒的估计噪声子空间，更稳定地获得最优信干噪比。运用本方法可以使得干扰子空间分解更为准确可靠从而获得更好的抗干扰效果。

Description

基于MWF中运用滑窗判定的噪声子空间估计方法

技术领域

本发明主要涉及到信号处理领域，特指一种基于MWF运用滑窗判定的噪声子空间分解方法。

背景技术

地球表面接收到的导航信号功率很小，以GPS为例，比接收机热噪声低约21dB。当受到人为干扰时，接收信号难以恢复。运用天线阵，并且进行时延扩展，形成空时二维模式，可以大幅增加自由度，增强了抗干扰能力。设阵元数为M，时间延迟数为N，则接收数据X为MN×1维向量，其协方差矩阵R_x＝E{XX^H}，为MN×MN维矩阵，其最优处理的运算量约为O((MN)³)，随着空时处理维数的增加，运算量成立方倍增长，大运算量对计算资源消耗和计算时长都是不可接受的。多级维纳滤波(MWF:Multistage Weiner Filter)有避免大矩阵分解和降维处理两个优势^[1]，适用于解决空时二维抗干扰处理的问题。

多级维纳滤波是维纳滤波器的一种多级等效形式，因此运用这一方法首先需要将抗干扰处理转化为维纳滤波的问题。通过广义旁瓣对消器可以从接收信号矢量分离出期望信号，形成维纳滤波的模式，如附图1所示。接收信号矢量经过上面的支路得到期望信号d₀(n)，d₀(n)包含有用信号和干扰，阻塞矩阵B₀＝null(h)，下支路通过B₀阻塞期望信号，则x₀(n)仅含干扰。

经典的多级维纳滤波器结构如附图2所示。按照信息论的思想，从每一步观测数据矢量x_i(n)中尽量提取与d_i(n)相关的信息，用x_i(n)与d_i(n)的归一化互相关矢量可以达到这个目的。x_i(n)与d_i(n)的互相关为其归一化矢量为阻塞矩阵为与h_i正交的[(MN-i)×(MN-i-1)]维矩阵，记为B_i，可表示为B_i＝null(h_i)，其列矢量构成h_i的零空间，即

经典算法中的综合滤波器是一个反推的过程。假设预定的降维维数为P，即分解到第P级停止，由图2可知，e_P(n)＝d_P(n)，第P级实际构成了经典的维纳滤波器，根据经典维纳滤波的理论可以推出这一步标量权值其中ξ_P＝E{|d_P(n)|²}，从P级往前推，第P-1级的误差为相应的可以算得w_P，以此类推可以算得标量形式的各阶维纳滤波器权值。

由各阶标量维纳滤波器可嵌套组成多级维纳滤波器的综合滤波器部分，求得各阶标量维纳滤波器权值之后，逐级反推可得到综合滤波器的权值

最终算得的权值为

w_MWF＝Tw (2)

经典MWF算法中的迭代次数P是通过每步迭代后均方误差ξ_i＝E{|e_i(n)|²}与一个门限值(通常设为接收机热噪声方差)相比较来确定，只要ξ_i小于门限值，就让迭代停止，使P＝i。运用经典MWF对空时二维模式抗干扰处理降维，会使最优子空间维数难以估计准确，有两个原因：1、噪声随环境变化，用一个固定的门限难以准确估计干扰子空间维数；2、因为MWF的每一步按照与接收信号矢量的最大相关提取干扰分量，导致每一步提取之后观测数据方差成L型变化，即先陡降，后缓降，致使后续接收信号残差方差与白噪声方差区分不明显，很难找准。在实际中运用经典MWF，迭代阶数很容易超过最优，将引起信噪比较大降幅，约有4～5dB^[2]，所以准确估计干扰的子空间是必要的。

参考文献：

[1]J.S.Goldstein,I.S.Reed.A Multistage Representation of the WienerFilter Based on Orthogonal Projection[J].IEEE Transactions on InformationTheory,1998,44(7):2943-2959.

[2]周柱,张尔扬,卢树军.基于噪声子空间估计的MWF实现[J].宇航学报,2012,33(5):661-668.

[3]D.C.Ricks,J.S.Goldstein.Efficient Architectures for ImplementingAdaptive Algorithms[J].Proceedings of the 2000Applications Symposium,2000:29-41.

[4]郭艺.GPS接收机空时抗干扰理论与实现关键技术研究[D].国防科技大学,2007。

发明内容

为了克服以上缺陷，本发明根据MWF的特性，结合小矩阵最大特征值运算，提出了一种运用滑窗定阶的方法。空时二维处理的主要问题是权值计算，图3即为本专利提出的一种运用滑窗定阶的噪声子空间估计方法的流程图。运用该方法完成最优权值的计算，如图1中“w”标识的部分。该方法创新点主要在于：将各阶期望信号组成矢量形式，求取该矢量数据的协方差矩阵，根据该矩阵高阶对角元素接近于白噪声方差的特性，提出一种滑窗的方法。该方法用一个小矩阵沿所得到的协方差矩阵对角线下滑，逐步比较小矩阵最大特征值与协方差矩阵的对角元素，显然，当特征值略大于白噪声对角线元素时，也就达到了干扰子空间和白噪声子空间的分界点。该方法利用了这一特性：各阶特征值区分度大大高于传统MWF方法中各阶最小均方误差。因此运用该方法可以鲁棒的估计噪声子空间，更稳定地获得最优信干噪比。

本发明包括如下步骤：

S1、将抗干扰处理问题转化为维纳滤波模式：将接收信号矢量通过广义旁瓣对消器，得到上支路的期望信号，记为d₀(n)，下支路x₀(n)为待处理信号矢量，如附图1，形成维纳滤波的模式；

S2、粗略估计干扰子空间维数：如图2所示，将x₀(n)通过MWF的分析滤波器，生成d₁(n)与x₁(n)，任取矢量x₁(n)中一路信号，可设其方差为再用同样的方法分解x₁(n)，以此类推可以得到d_i(n)，x_i(n)，i＝1,2,…,MN-1。白噪声方差为σ²，将序列与k(k一般取3～5)倍σ²进行比较，当时，将此时的i作为粗略估计的噪声子空间维数，设为P₁，P₁几乎总是大于准确的子空间维数P；

S3、估计干扰子空间维数：运用MWF分析滤波器结构将接收信号矢量映射为各阶期望信号d_i(n)，i＝1,2,…,MN-1，将各阶d_i(n)组成矢量形式，设为x_T(n)＝[d₁(n)：d₂(n)，…，d_MN-1(n)]^T，求取x_T(n)的协方差矩阵，设为沿的对角线，从第P₁阶开始截取一个小矩阵，设为R_j，用幂法算得R_j的最大特征值，设为λ_j，将λ_j与1.5σ²比较，当λ_j小于1.5σ²时，认为此时的j即是干扰子空间维数，也就是最佳的迭代次数，令j＝P；

S4、执行MWF的综合滤波部分：在获得最佳迭代次数P之后，从第P级开始逆推，由图2可知，e_P(n)＝d_P(n)，根据经典维纳滤波的理论可以推出这一步标量权值以此类推可以算得标量形式的各阶维纳滤波器权值，逐级反推可得到综合滤波器的权值，如(1)式，并得到最终的权值，如(2)式。用(2)式所示权值对接收信号进行滤波。

本发明的优点在于：本发明根据MWF的特性，结合小矩阵最大特征值运算，提出了一种运用滑窗定阶的方法，用该方法进行噪声子空间判断区分度更大，可以鲁棒的估计噪声子空间，更稳定地获得最优信干噪比。运用本方法可以使得干扰子空间分解更为准确可靠从而获得更好的抗干扰效果。文献[2]将接收信号矢量映射为各阶期望信号组成的矢量，生成相应的协方差矩阵，将该矩阵的高阶特征值带入MDL准则进行运算，得到了噪声子空间维数。而本发明不需要算得所有高阶特征值，而且避免在文献[2]中出现的对数运算，便于实现。

附图说明

图1为广义旁瓣对消器结构图；

图2为多级维纳滤波器结构图；

图3为运用滑窗判定的噪声子空间分解方法流程图；

图4为输出SINR与迭代次数关系图；

图5为经典MWF与本发明在最佳迭代次数判定上的比较。

具体实施方式

下面结合附图，对本发明进行详细描述：

空时二维处理的主要环节是权值计算，附图3为在MWF理论基础上提出的一种运用滑窗定阶的噪声子空间估计方法。运用该方法完成最优权值的计算，如图1中“w”标识的部分。运用该方法进行空时二维处理的步骤如下：

S1、将抗干扰处理问题转化为维纳滤波模式：将接收信号矢量通过广义旁瓣对消器，得到上支路的期望信号，记为d₀(n)，下支路x₀(n)为待处理信号矢量，如附图1，形成维纳滤波的模式。

S2、粗略估计干扰子空间维数：根据附图2，按照最大相关并归一化算得h₁，即第一阶的归一化相关矢量，如上文所述，阻塞矩阵B₁为h₁的零空间矩阵，将x₀(n)通过这一级滤波，得到用第一阶期望信号d₁(n)、提取d₁(n)之后，生成的数据矢量x₁(n)，及矢量x₁(n)中任意一路的方差依此类推，构建各级h_i和B_i，生成d_i(n)、x_i(n)、i＝1,2,…,MN-1，白噪声方差为σ²，将与k(k为常数，一般取3～5)倍σ²进行比较，当时，将此时的i作为粗略估计的噪声子空间维数，设为P₁，P₁几乎总是大于准确的子空间维数P；

阻塞矩阵B_i的选取不是唯一的，需要满足的条件是与MWF滤波每一级的互相关矢量h_i正交，运用文献[3]提出的方法，令阻塞矩阵为其中I为单位矩阵，可推导出满足阻塞矩阵的要求。因为所以有

由上式可知，可以分为两步：第一步为第二步将d_i(n)乘上h_i，得到接收信号矢量在第i维子空间的投影。减去这个投影，就完成了一次分解；

降维多级维纳滤波的过程相当于接收信号矢量经过降维矩阵的作用映射为一个低维数据矢量。降维矩阵的列矢量为根据阻塞矩阵B_i的形式可以把降维矩阵的列矢量化简；

已知t₁＝h₁，可见如果h_i之间相互正交，则有t₂＝h₂，以此类推可以得到t_i＝h_i。首先证明[h₁,h₂,…,h_P]为一组标准正交向量，已知分母为常量，在正交性分析中可以不考虑，令有

因为x_i-1(n)＝B_i-1x_i-1(n)，所以有

依此类推，可以得到

由此可得h_i相互正交，因此t₂＝h₂，依此类推，t_i＝h_i，所以用替代原阻塞矩阵B_i＝null(h_i)后，降维矩阵也可表示为T＝[h₁,h₂,…,h_P]；

S3、准确估计干扰子空间维数：如附图3所示，在MWF框架下，根据子空间理论，结合矩阵最大特征值的计算，运用滑窗逐级搜寻噪声子空间维数。现先给出其操作步骤，再给出其理论依据。这一步分解为以下四步执行：

S31、协方差矩阵计算：将接收信号矢量经过MWF生成的各阶期望信号d_i(n)组成矢量形式x_T(n)＝[d₁(n),d₂(n),…,d_MN-1(n)]^T，求取该矢量的协方差矩阵为三对角矩阵；

S32、截取小矩阵：在中用一个窗口截取一个小方阵，设为R_j，不妨设方阵宽度为width，j为方阵左上角元素在中的阶次(首次截取的位置j＝P₁，也就是粗略估计的噪声子空间维数)；

S33、求取最大特征值：用幂法求得矩阵R_j的最大特征值(R_j为三对角矩阵，且维度小，运算量小)，设为λ_j；

S34、判定噪声子空间维数：将λ_j与1.5倍σ²比较，如果λ_j大于1.5σ²，则j＝j+1，将窗口沿着待分析矩阵的对角线向右下滑动一阶，返回到S42。如果不满足λ_j大于1.5σ²，则认为此时的j为噪声子空间维数，即为P；

S3步的理论依据可概括如下：将的主对角线元素简记为副对角线元素简记为δ_i，可将协方差矩阵表示如下：

用MWF逐步提取相关干扰信息过程中，必然存在P＜MN-1，经过P步投影，观测数据中不再有干扰成分。文献[4]证明：如果干扰子空间维数为P，则多级维纳滤波器做P级分解后的(MN-P-1)组观测数据矢量为白噪声矢量，且为：

式中，(·)^(P)表示干扰子空间维数为P的情况，下标(·)_i表示MWF的第i级，v(n)为白噪声矢量；

现根据(7)式对矩阵进行分析。其副对角线元素可推导如下

当i≥P时，因此所有含有因子h_i的项均与正交，而与h_i正交，理论上可以推出δ_i＝0；

易知d_i(n)观测数据矢量的加权，如上文所述，P阶之后的观测数据矢量为白噪声，显然P阶之后观测数据矢量的加权d_i(n)也是白噪声。可知，当i＞P时，主对角元素约等于白噪声方差；

因此在理论上，可以按照下面形式进行分块：

左上部分仍为P阶三对角矩阵，右下部分为对角矩阵，其余位置元素为零。矩阵的特征值为的解，有下式

式中由此可知，高于P阶的的特征值等于白噪声方差；

可以分块，以噪声子空间维数P为界，左上部分为三对角矩阵，直到第P阶，右下部分为对角矩阵，对角元素约等于白噪声方差。因此，如果用一个窗口从的左上向右下滑动，当窗口内矩阵R_j的最大特征值接近于白噪声方差时，说明窗口刚好滑动到噪声子空间和白噪声分界的阶次，也就找到了准确的噪声子空间维数；

S4、权值计算：运用经典MWF综合滤波器反推算得最优权值，获得最佳迭代次数P之后，从第P级往前推算。由图2可知，e_P(n)＝d_P(n)，推出这一步标量权值其中ξ_P＝E{|d_P(n)|²}，从P级往前推，第P-1级的误差为相应的可以算得w_P-1，以此类推可以算得标量形式的各阶维纳滤波器权值。由各阶标量维纳滤波器可嵌套组成多级维纳滤波器的综合滤波器部分，如(1)式，最终算得的权值如(2)式。

仿真试验

采用5元加芯圆阵，成十字型排列，半径为d＝λ/2，λ为接收信号波长，时间延迟单元数N＝5。无干扰情况下的信噪比为-20dB。按照干噪比30dB设置两个功率相同的部分带宽干扰和两个单频干扰。进行两组仿真：一、试验输出信干噪比(SINR)与迭代次数的关系；二、检测本发明对于噪声子空间估计的作用。仿真中均重复多次取平均。

一、输出信干噪比(SINR)与迭代次数的关系

运用经典MWF方法进行抗干扰滤波，逐渐增加MWF的迭代次数，观测输出SINR的变化。结果如附图4所示，表明随着迭代次数增加，SINR在达到最大值以后将大幅下降，约有5dB，与文献[2]结果相符。

二、本发明对于噪声子空间估计的作用

仿真中将沿对角线截取的小方阵边长width定为3。图5的左图表示经典MWF中各步迭代运算之后接收信号矢量任一路方差的变化。因为白噪声方差用分贝表示为0dB，从图5(左)可以看出，在12次迭代之后，干扰信号提取完，进行第13次迭代之后，均方误差与白噪声比值降了3dB。

图5的右图表示运用本发明，在粗略判断最佳迭代次数(本场景中P₁＝10)之后，滑窗用幂法依次算得的小矩阵最大特征值与白噪声方差的比值。因为场景相同，且与运用经典MWF在同一次仿真中，所以无疑最佳迭代次数为12。而在图5(右)中，第12和13阶比值的差距有7dB。这就意味着区分度大大高于经典方法对最佳迭代次数的区分度。因此可以得出结论，运用本方法可以使得干扰子空间分解更为准确可靠从而获得更好的抗干扰效果。

Claims (2)

1.基于MWF中运用滑窗判定的噪声子空间估计方法，将各阶期望信号组成矢量形式，求取该矢量数据的协方差矩阵，根据该矩阵高阶对角元素接近于白噪声方差的特性，提出一种滑窗方法，该滑窗方法用一个小矩阵沿所得到的协方差矩阵对角线下滑，逐步比较小矩阵最大特征值与协方差矩阵的对角元素，当特征值略大于白噪声对角线元素时，也就达到了干扰子空间和白噪声子空间的分界点，其特征在于，包括如下步骤：

S1、将抗干扰处理问题转化为维纳滤波模式：将接收信号矢量通过广义旁瓣对消器，得到上支路的期望信号，记为d₀(n)，下支路x₀(n)为待处理信号矢量，形成维纳滤波的模式；

S2、粗略估计干扰子空间维数：按照最大相关并归一化算得h₁，阻塞矩阵B₁为h₁的零空间矩阵，将x₀(n)通过这一级滤波，得到d₁(n)、x₁(n)，及矢量x₁(n)中任意一路的方差依此类推，构建各级h_i和B_i，生成d_i(n)、x_i(n)、i＝1,2,…,MN-1，白噪声方差为σ²，将与k倍σ²进行比较，k取3～5，当时，将此时的i作为粗略估计的噪声子空间维数，设为P₁，P₁总是大于准确的子空间维数P；

S3、估计干扰子空间维数：运用MWF分析滤波器结构将接收信号矢量映射为各阶期望信号d_i(n)，i＝1,2,…,MN-1，将各阶d_i(n)组成矢量形式，设为x_T(n)＝[d₁(n),d₂(n),L,d_MN-1(n)]^T，求取x_T(n)的协方差矩阵沿的对角线，从第P₁阶开始截取一个小矩阵，设为R_j，用幂法算得R_j的最大特征值，设为λ_j，将λ_j与1.5σ²比较，当λ_j小于1.5σ²时，认为此时的j即是干扰子空间维数，也就是最佳的迭代次数，令j＝P；

S4、执行MWF的综合滤波部分：在获得最佳迭代次数P之后，从第P级开始逆推，e_P(n)＝d_P(n)，根据维纳滤波理论推出这一步标量权值以此类推得到标量形式的各阶维纳滤波器权值，逐级反推得到综合滤波器的权值。

2.根据权利要求1所述的基于MWF中运用滑窗判定的噪声子空间估计方法，其特征在于，所述阻塞矩阵B₁的选取不是唯一的，满足的条件是与MWF滤波每一级的互相关矢量h_i正交，令阻塞矩阵为推导出满足阻塞矩阵的要求；

因为所以有

由上式可知，可以分为两步：第一步为第二步将d_i(n)乘上h_i，得到接收信号矢量在第i维子空间的投影，减去这个投影，就完成了一次分解。