CN107179521B - 一种低复杂度的鲁棒子空间估计方法 - Google Patents

Abstract

本发明属于信号处理与通信领域，涉及一种低复杂度的鲁棒子空间估计方法。本发明的方法为了减小计算复杂度，主要是通过将未知变量分解为一系列块变量并选择块更新规则，然后找到适当的对应于原始目标函数的上界函数对目前进行求解，本发明的方法基于不确切的块坐标梯度下降，对每一个子问题，通过找到一个可靠的逼近解且仅仅通过矢量的四则运算，避免了高复杂度的一维线搜索优化问题的计算开销。同时，通过引入一种新的可调参的鲁棒函数，既避免了代价函数的不平滑性而带来的计算问题，又可以提高对异常点或脉冲噪声的鲁棒性。

Description

一种低复杂度的鲁棒子空间估计方法

技术领域

本发明属于信号处理与通信领域，涉及一种低复杂度的鲁棒子空间估计方法。

背景技术

子空间方法在信号处理、通信以及计算机视觉中具有重要的应用，其主要作用是通过选取低维不相关的变量成分来表示高维数据，以便于达到降噪和降维的效果。传统子空间方法大多为基于欧式范数(l₂)空间下的低秩矩阵分解，例如基于截断奇异值或特征分解的PCA(principal component analysis)方法；《Multiple emitter location andsignal parameter estimation》。众所周知，这种传统基于l₂范数的方法不能够适用于脉冲噪声，因为l₂范数仅在高斯背景噪声下具有最优性。事实上，在许多实际应用中，噪声成分除了包含一些相对较小幅度的高斯扰动之外，还包含了一小部分的奇异点(Outliers)。这些奇异点通常是由于传感器的异常扰动、放大器饱和以及一些恶意的攻击或干扰等因素而造成的。

近年来，研究适用于脉冲噪声下的子空间估计受到了学者们的广泛关注。由于l₂范数对脉冲噪声不具有鲁棒性，学者们开始尝试利用其它对脉冲噪声不具有敏感性的代价函数。最常用的代价函数为l₁罚函数。首先，它具有比l₂函数更低罚函数值特性；另外，它的凸函数特性使得在算法求解上也比较容易实现。《Optimal algorithms for L₁-subspacesignal processing》利用基于l₁范数的投影最大化算法；然而，该算法的计算复杂度非常高(随着数据维数增加而呈指数型增长)，因此不实用于较大的数据模型。《Practical low-rank matrix approximation under robust L₁-norm》利用交替优化求解低秩矩阵分解的方法分离出子空间，主要通过线性或二阶锥规划求解每一个凸子优化问题。《L_p-MUSIC:Robust direction-of arrival estimator for impulsive noise environments》利用l_p(1＜p＜2)范数为代价函数求解复数数据条件下的低秩矩阵分解，但仍然使用交替凸优化方法求解最小化问题，不同的是子问题通过复牛顿迭代算法求解。然而，该方法的计算复杂度仍然相对较高，因为每一个子问题都需要搜索最优的步长和求解Hessian矩阵(包括矩阵求逆运算)。另外，与L₁范数和l_p(1＜p＜2)范数相比较，l_p(1＜p＜2)范数对脉冲噪声更具有鲁棒性。

发明内容

本发明的目的在于克服现有技术的不足，从低秩矩阵分解的观点来求解子空间，提出了一种低复杂度的鲁棒子空间估计方法，具体为高效的迭代重加权算法，该方法不需要任何的步长搜索，每次迭代都是一个子问题的最优闭式解。该方法能够适用于多种不同类型的脉冲噪声，为子空间估计提供了靠性保证。

本发明的技术方案是：

一种低复杂度的鲁棒子空间估计方法，该方法用于对脉冲噪声兼具鲁棒性的子空间估计，其特征在于，包括以下步骤：

S1、获取接收信号如下公式1：

Y＝Y_t+E (公式1)

公式1中，为无污染的测量信号，E为脉冲噪声；

S2、根据鲁棒子空间估计的目的，即找到和R＜＜min{M,N}，且使得Y_t尽可能的逼近UV^H，如下公式2：

Y_t≈UV^H (公式2)

公式2中，(·)^H表示共轭转置运算；

将目标U和V分解为2R个块变量，且每个块变量为U和V的列矢量，采用如下迭代方法获取U和V：

S21、设初始迭代数k＝0，并随机初始化初始化权重矩阵W为空矩阵；

S22、将第k次迭代估计值U(k)和V(k)赋给中间变量和

S23、通过公式r＝(k mod R)+1获取每次要更新的列r；其中，mod表示求余运算；

S24、根据获得的列r，通过如下公式3获得矩阵H_r：

公式3中，和分别表示矩阵和的第i个列向量；

S25、对权重矩阵W进行更新：

首先，将权重矩阵W更新为然后，如果|W_ij|≤ε，那么将W_ij更新为否则，W_ij更新为其中，W_ij为W的第i行、j列元素，ε>0，p∈(0,2]；

S26、对列向量u_r(k+1)进行更新：

采用如下公式4更新列向量u_r(k+1)中的M个元素：

其中，[u_r(k+1)]_i表示U(k+1)的第r个列向量的第i个元素，1≤i≤M，H_r[i,:]表示矩阵H_r的第i行，W[i,:]表示矩阵W的第i行，为Hadmard积；

S27、采用如步骤S25的方法再次更新权重矩阵W：

首先，将权重矩阵W更新为然后，如果|W_ij|≤ε，那么将W_ij更新为否则，W_ij更新为

S28、对列向量v_r(k+1)进行更新：

采用如下公式5更新列向量v_r(k+1)中的N个元素：

公式5中，[v_r(k+1)]_j表示V(k+1)的第r个列向量的第j个元素，H_r[:,j]表示矩阵H_r的第j列，W[:,j]表示矩阵W的第j列；

S29、采用如下公式6更新U(k+1)，V(k+1)：

S210、更新迭代次数k:k＝k+1；

S211、判断迭代次数k是否达到预设值，若是，则进入步骤S3，若否，则返回到步骤S22继续迭代更新；

S3、将S2中的输出目标值U(k)和V(k)分别赋给U和V；Orth(U)为所需的列子空间，其中Orth表示正交化。

为了处理常规范数不具有二阶可微分性，本发明利用文献《Sparsegeneralized eigenvalue problem via smooth optimization》中的改进的l_p范数

其中，|·|表示绝对值或求模运算(对复数变量)，ε>0为平滑参数(同时提供了鲁棒性)。

基于通过求解如下最优化问题找到所需的子空间：

其中，γ>0为正则化参数，||·||_F为Frobenius范数，且

其中，[·]_ij表示对应矩阵的第(i,j)元。

直接求解优化问题是非常困难的，特别是对于较高的数据维数，即当M或N很大时，每一步迭代中使用常规的梯度下降法都需要搜索对应的步长。为了减小计算复杂度，一种减少此计算复杂度的有效策略是将待优化变量分解为比较小的块变量，然后对相应较容易求解的字问题进行优化。为此，本发明提出了一种有效的块迭代重新加权算法来求解问题，并且，每个子问题都有唯一的闭式解。因此，相对于常规的块坐标下降(BCD，BlockCoordinate Descent)方法而言，此方法还不需要任何步长的搜索。

本发明的法涉及两个关键：1)将未知变量分解为一系列块变量并选择块(索引)更新规则；2)找到适当的对应于原始目标函数的上界函数。首先，我们将未知变量(U，V)分解成2R个块成分，且每个块变量为U和V的列矢量。采用这种块分解策略的主要原因是其可以导致一种高效的算法，且该算法在每次迭代中只需要基本的向量四则运算操作，因此具有较低的计算复杂度。具体来说，假设在当前迭代的估计为(U(k),V(k))，则在下一次迭代中，我们只需要更新两个列矢量变量，其分别为u_r(k+1)和v_r(k+1)，其中mod为求余运算。这种块更新的方式称之为round-robin(circle)fashion更新规则。在数学上，下一次迭代估计为求解如下两个优化问题：

其中，

值得注意的是，虽然上式仍然没有一个闭式解，但是本发明提出的方法可以找到这两个优化问题的逼近解，从而避免步长因子的搜索问题。

本发明的有益效果是：既具有对各种异常点噪声的鲁棒性，又不需要步长的搜索且具有较低的复杂度，能够适用于多种不同类型的脉冲噪声，为子空间估计提供了靠性保证。

附图说明

图1为数据降维的方向图实例。

图2为信源定位的空间谱图实例。

具体实施方式

实例1下面以鲁棒(实数)数据降维为例子，详细阐述本发明中所提出的鲁棒子空间方法的具体应用。首先，我们产生150个二维真实数据样本(如图1中的五角星)，且服从如下高斯分布：

其中，表示正太分布，0₂(2维列矢量)为均值，表示其方差。同时，产生5个异常数据点，o₁,…,o₅，如图1中的左上角圆圈点。同时，5个随机异常数据点服从如下高斯分布：

在此实例中，将算法1中的Y定义为同时，设ε＝1，γ＝0.5，p＝1。U(0)和V(0)的初始化分别为常规PCA(奇异值分解)得到的结果，设最大的迭代次数为200，R＝1。最终结果如图1所示，据图可知，本发明提出的方法与真实数据的主方向更接近，这说明了本发明提出的方案具有对异常点数据的鲁棒性。ACO方法为文献《L_p-MUSIC:Robustdirection-of arrival estimator for impulsive noise environments》中的方案；常规PCA方法具有最差的性能，这是因为l₂对异常点不具有鲁棒性而造成的。

实例2下面阐述本发明方案在传感器阵列信源定位中的应用。设K个波达方向角为θ＝[θ₁,...,θ_K]的窄带信源位于M(M>K)个均匀线性阵列结构(uniform linear array,ULA)的远场，则N个快拍阵列接受信号模型为

x(t)＝A(θ)s(t)+n(t),t＝1,…,N

式中，s(t)为零均值信号矢量，为阵列流型矢量，a(θ_k)为M维列矢量且其元素为m＝1,...,M，θ_k∈(-π/2,π/2)，λ与d分别表示信号波长与阵元间距。本案例中，设天线数M＝8，信源数K＝4，快拍数N＝100，信源位置为θ＝[-20°,0°,6°,35°]^T，且d＝λ/2；叠加噪声为高斯混合脉冲噪声，信噪比为-5dB。为了使本发明方案适用于信源定位，先将数学符号与发明内容中符号进行对应，设接收信号样本为R＝K＝4。本案例中，设p＝γ＝ε＝0.5。在得到最终的子空间Orth(U)之后，再利用类似于文献《Multiple emitter location andsignal parameter estimation》中标准MUSIC算法实现信源定位。实验结果如图2所示，据图可知，本发明方案具有更好的谱分辨率特性，分别在-20°,0°,6°,35°这4个位置都出现了较锐利的谱峰。特别地，其它两种方案在0°,6°这两个较近位置只出现了一个谱峰，即不具有分辨能力。

Claims (1)

1.一种低复杂度的鲁棒子空间估计方法，该方法用于对脉冲噪声兼具鲁棒性的子空间估计，其特征在于，包括以下步骤：

S1、获取接收信号如下公式1：

Y＝Y_t+E (公式1)

公式1中，为无污染的测量信号，E为脉冲噪声；

S2、根据鲁棒子空间估计的目的，即找到和R＜＜min{M,N}，且使得Y_t尽可能的逼近UV^H，如下公式2：

Y_t≈UV^H (公式2)

公式2中，(·)^H表示共轭转置运算；

将目标U和V分解为2R个块变量，且每个块变量为U和V的列矢量，采用如下迭代方法获取U和V：

S21、设初始迭代数k＝0，并随机初始化初始化权重矩阵W为空矩阵；

S22、将第k次迭代估计值U(k)和V(k)赋给中间变量和

S23、通过公式r＝mod(k,R)+1获取每次要更新的列r；其中，mod表示求余运算；

S24、根据获得的列r，通过如下公式3获得矩阵H_r：

公式3中，和分别表示矩阵和的第i个列向量；

S25、对权重矩阵W进行更新：

首先，将权重矩阵W更新为然后，如果|W_ij|≤ε，那么将W_ij更新为否则，W_ij更新为其中，W_ij为W的第i行、j列元素，ε>0，p∈(0,2]；

S26、对列向量u_r(k+1)进行更新：

采用如下公式4更新列向量u_r(k+1)中的M个元素：

其中，[u_r(k+1)]_i表示U(k+1)的第r个列向量的第i个元素，1≤i≤M，H_r[i,:]表示矩阵H_r的第i行，W[i,:]表示矩阵W的第i行，⊙为Hadmard积；

S27、采用如步骤S25的方法再次更新权重矩阵W：

首先，将权重矩阵W更新为然后，如果|W_ij|≤ε，那么将W_ij更新为否则，W_ij更新为

S28、对列向量v_r(k+1)进行更新：

采用如下公式5更新列向量v_r(k+1)中的N个元素：

公式5中，[v_r(k+1)]_j表示V(k+1)的第r个列向量的第j个元素，H_r[:,j]表示矩阵H_r的第j列，W[:,j]表示矩阵W的第j列；

S29、采用如下公式6更新U(k+1)，V(k+1)：

S210、更新迭代次数k:k＝k+1；

S211、判断迭代次数k是否达到预设值，若是，则进入步骤S3，若否，则返回到步骤S22继续迭代更新；

S3、将S2中的输出目标值U(k)和V(k)分别赋给U和V；Orth(U)为所需的列子空间，其中Orth表示正交化。