CN107450045B - 基于focuss二次加权算法的doa估计方法 - Google Patents

Abstract

本发明公开了一种基于FOCUSS二次加权算法的DOA估计方法，针对信源DOA估计问题，将传统信号模型扩展为稀疏信号模型，基于稀疏信号模型，在FOCUSS框架基础上采用二次加权算法进行迭代计算，在每一次迭代过程中进行两次数据重构；其中，第二次数据重构利用初步重构数据进行加权矩阵二次赋值，获得二次加权矩阵，并对输入信号再次重构，得到更加精确的信号重构结果；通过不断逼近目标信号值直至迭代终止，获得最终收敛解。在最终收敛解中可以找到波达方向。本发明比DOA估计的传统方法、最小l₁范数方法、FOCUSS方法更为稳健和精确。

Description

基于FOCUSS二次加权算法的DOA估计方法

技术领域

本发明涉及信号处理领域，尤其涉及一种基于FOCUSS框架的二次加权算法的波达方向(DOA)估计方法。

背景技术

传统的DOA估计方法(如MUSIC方法、ESPRIT方法等)，在信源不相干的前提下利用空间分解后的正交性构造出“针状”空间谱，若信源相干则会严重影响分辨和测向精度。空间中信源个数有限且一般小于阵列传感器数量，符合空间“稀疏”定义，因此稀疏重构算法可应用于DOA估计。由于稀疏空间中只有“0”或“1”的区别，而无需区分相干源或非相干源，也不会用到噪声空间和信号空间的正交特性，因此，在处理DOA估计问题时具有极强的稳健性；同时，稀疏重构方法只是在它认为信号可能存在的精确位置不断地聚集能量同时抑制无信号区域的能量，从而达到高分辨效果。

目前，可用于DOA估计的稀疏重构算法有很多，包括利用每步迭代中残差与采样矩阵的匹配运算后剪枝逼近最优解的贪婪算法(如OMP、CoSaMP等)；用SVD方法将观测矩阵分解，在保持信号子空间维度的前提下对噪声子空间降维以降低运算量的

方法；以及在每步迭代中令加权

范数最小的FOCUSS方法。贪婪算法虽然运算量低、实现简单，但是对信号模型传感矩阵的RIP(Restricted Isometry Property)性要求高，而且在稀疏项数量较少情况下才能保证较高的恢复成功率；采用最小

范数原理的

等方法，基于凸优化方法进行运算和重构，算法复杂、运算速度慢、运算硬件成本高等特性使之不适用于大规模计算；FOCUSS及其衍生算法类中，原始FOCUSS方法重构精度不高，TLS-FOCUSS方法不能用于MMV模型的DOA估计，SD-FOCUSS方法不能对DOA模型中模型噪声进行有效建模而降低了估计性能。

本发明继承了FOCUSS方法收敛速度快、计算量低的特点，对“稀疏”的限制和传感矩阵RIP的限制更为宽松，通过重新加权的

范数最小化过程有效提高了DOA估计的精度。

发明内容

有鉴于此，本发明针对信源DOA估计问题，将传统信号模型扩展为稀疏信号模型，利用FOCUSS方法的处理原理，在每次迭代中将上一次迭代的估计值作为本次迭代的加权值，不断逼近目标信号值直至迭代终止，比DOA估计的传统方法、最小

范数方法、FOCUSS方法更为稳健和精确。

为了解决上述技术问题，本发明是这样实现的：

步骤(1)构建稀疏信号模型

将传统的DOA估计信号模型，通过对无目标区域的填“0”扩展，构建成相应的DOA估计稀疏信号模型。其步骤为：

步骤(1.1)构建传统模型：考虑一个M-阵元的信号接收阵列，假设U个信号源分别以θ₁,…,θ_U的入射角投射到阵列接收机上，则离散后的信号观测模型表达为：

y(l)＝Φs(l)+n(l)＝[a(θ₁),…,a(θ_U)][s(1,l),…,s(U,l)]^T+n(l),

其中，l∈[1,…,L]代表第l个快拍单元，L为总快拍数，y(l)∈C^M是阵列输出的第l个快拍，s(u,l)为第u个信源的信号幅度，a(θ_u)为第u个信源的信号导引矢量，n(l)为信道中的噪声向量；y(l)、a(θ_u)、n(l)均为M×1的列向量；

步骤(1.2)等价填“0”扩展：对除s(1,l),…,s(U,l)之外的空间中无目标区域进行填“0”处理，同时扩展对应观测矩阵Φ，得到对应稀疏信号模型y(l)＝Ax(l)+n(l)，其中

为空间方向全角度或区域角度的网格划分，N为划分的格点数量；x(l)＝[x(1,l),…,x(N,l)]^T；当

等于某一信源真实的入射方向θ_u(u∈[1,…U])时，x(n,l)的值不为0，否则其值为0，即x(n,l)的值为

其构造原理如图2所示；

步骤(1.3)构建多快拍稀疏信号模型：Y＝AX+N，其中Y＝{y(l)|_l＝1,…,L}，X＝{x(l)|_l＝1,…,L}，及N＝{n(l)|_l＝1,…,L}。

步骤(2)参数初始化

步骤(2.1)求取信号估计初始值。为降低运算量，提高计算效率，初始值的位置应尽量接近真值，其选取公式为：X₀＝A^H(AA^H)^-1Y。

步骤(2.2)设置收敛因子p。p为(0,1]之间任意的常数，p的值越小，估计结果的稀疏性就越好，运算量随之增大；p的值越大，估计结果的稀疏性就越差，运算量随之降低(仿真实验中p＝0.5)。

步骤(2.3)求取噪声因子λ。其求取公式为：

其中σ为已知的观测噪声方差；β为系统的贝叶斯广义方差，

Γ(·)为标准gamma分布函数。

下面采用步骤(3)～(6)进行FOCUSS算法的迭代计算，在每一次迭代过程中进行两次数据重构。

步骤(3)求取迭代中初次加权矩阵

在第k次迭代中，上次迭代求出的估计结果为X_k-1，则第k次迭代中的初次加权矩阵为

其中

下角标k-1均表示上次迭代求出的估计结果；

步骤(4)迭代中初步估计

在第k次迭代中，X初步估计

的求取方法为：

其中，

步骤(5)求取迭代中二次加权矩阵

在第k次迭代中，求取能改善估计性能的二次加权矩阵W_k，其求取方法为：首先令W_k为N×N的对角矩阵，对角元素ω_n的取值为：若n属于坐标集U，则ω_n的幅度为1；否则，根据

幅度的大小设置ω_n，

幅度越大，ω_n的幅度值越大，ω_n幅度值取值范围为[0,1)。

在一个优选实施例中，令对角元素的值满足以下要求

其中坐标集

为序列

中前2U个幅度最大值对应的基坐标，U是步骤(1.1)中的信源个数。

是通过对

的行向量分别求取

范数得到的。

步骤(6)迭代中二次估计

在第k次迭代中，X二次估计X_k的求取方法为：

其中，A_k＝AW_k。

步骤(7)验证迭代终止条件

如果

则迭代终止，转入步骤(8)；否则，迭代次数加1，即k＝k+1，转入步骤(3)。其中ε为预设的终止门限(仿真实验中ε＝0.01)。

步骤(8)查找X中的非零行

其方法如下：对X_k的行向量分别求取

范数，即

找出其中前U个最大值对应的基坐标T＝[t₁,…,t_U]，就是X中非零行的行坐标。

步骤(9)找出目标源的入射角

在步骤(1.2)中构造的传感矩阵A的列坐标集中，查找T所在位置对应的角度值

即为本方法所估计的目标源入射角。

有益效果：

(1)本发明通过构建DOA估计的稀疏信号模型，并对FOCUSS方法通过重新加权进行改进，提高了DOA估计的稳定性和精度。

(2)本发明设计了二次加权矩阵的具体求取方式，以

中各元素的值与

最大幅值的比值作为权值，简单、准确的实现了权值的表达，便于二次加权矩阵的求解。

附图说明

图1是基于FOCUSS重新加权算法的DOA估计方法流程图。

图2是稀疏信号模型构建原理图，图2(a)是不含稀疏信号的观测模型，图2(b)是与图2(a)等价的含稀疏信号的观测模型。

图3是本发明方法与传统DOA估计(MUSIC)方法对比，图3(a)是2个非相干信源的估计结果，图3(b)是2个相干信源的估计结果。

图4是本发明方法与其他稀疏重构算法DOA估计的RMSE(均方根误差)曲线图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明将传统的DOA估计模型构建为具有稀疏性的信号处理模型，对压缩传感/感知技术应用于DOA估计的方法进行拓展和改进，在信号处理的每次迭代运算中，利用FOCUSS方法求解出最小

范数加权矩阵，并对输入信号进行初步重构；再利用初步重构数据进行改善性能的加权矩阵二次赋值并对输入信号再次重构，得到更加精确的信号重构结果，最后利用重构信号的行稀疏性与传感矩阵列向量之间的映射关系进行DOA估计。基于FOCUSS重新加权算法的DOA估计方法流程如图1所示。具体实现步骤如下：

步骤(1)：构建稀疏信号模型

考虑一个M-阵元的信号接收阵列，如果有U个信源分别以θ₁,…,θ_U的入射角投射到阵列接收机上，则离散后的信号观测模型表达为

y(l)＝As(l)+n(l)＝[a(θ₁),…,a(θ_U)][s(1,l),…,s(U,l)]^T+n(l),

其中，l∈[1,…,L]代表第l个快拍单元，L为总快拍数，y(l)是阵列输出的第l个快拍，s(u,l)为第u个信源的信号幅度，

为第k个信源的信号导引矢量，d为阵元间距，λ₀为载波波长，n(l)为信道中的噪声向量；y(l)、a(θ_u)、n(l)均为M×1的列向量；

上述模型中并未体现出稀疏特性，然而空间中信源个数是有限的，可认为没有信源的空间中其信号强度为0。通过对除s(1,l),…,s(U,l)之外的空间中无目标区域进行填“0”处理，同时扩展对应观测矩阵A，得到对应稀疏信号模型

y(l)＝Ax(l)+n(l)

其中

为空间方向全角度(或区域角度)的网格划分，N为划分的格点数量，

表示第n个网格对应的空间角度；x(l)＝[x(1,l),…,x(N,l)]^T，且x(n,l)|_n＝1,…,N为分别对应每个网格的信源信号幅度，若网格对应的空间角度

有真实的信源信号s(u,l)射入，则该网格的值x(n,l)为相应入射角所对应信源的信号值s(u,l)，对于没有信源信号的网格，其值为0，具体为：

重建后的多快拍稀疏信号模型为

Y＝AX+N

其中Y＝{y(l)|_l＝1,…,L}，X＝{x(l)|_l＝1,…,L}，N＝{n(l)|_l＝1,…,L}。若某个算法能对X进行重构恢复，或找到X中非零行的位置T＝[t₁,…,t_U]，就能够在观测矩阵A中找到对应的信号导引矢量

即可对该信源空间中的波达方向做出估计。

步骤(2)：参数初始化

在贝叶斯理论观点中，X可认为是独立于噪声N的随机、具有行稀疏性的未知矩阵。在此假设下，可得到X的最大后验估计

式中，P(X|Y)表示在Y发生的条件下X发生的概率，P(Y|X)表示在X发生的条件下Y发生的概率，P(X)表示X发生的概率。假设噪声N中各元素是服从方差为σ²的高斯白噪声，则可知

式中C表示常数。

文献1(Rao B D,Engan K,Cotter S F,Palmer J and Kreutz-Delgado K.Subsetselection in noise based on diversity measure minimization.IEEE Trans.SignalProcess.,2003,51(3):760–770)中讨论了稀疏向量的分布假设问题，认为稀疏向量x服从独立同分布的广义高斯分布，给出了x的概率分布函数：

其中，p为预设的收敛因子，

Γ(·)代表标准gamma分布函数。当p→∞时，该概率分布函数接近于均匀分布；当p→0时，该概率分布函数接近于一个冲击响应函数。根据DOA估计的实际情况，X相同行间的元素高度相关，本发明可将X的概率分布函数估计为：

其中，

综上可得X的最大后验估计为：

其中，λ为噪声因子且

为降低运算量，提高计算效率，初始值的位置应尽量接近真值，其选取方法为：X₀＝A^H(AA^H)^-1Y。

步骤(3)：求取迭代中初次加权矩阵

步骤(2)中最优解需满足的必要条件是J(X)对X的偏导为0，即

其中,

可得到

((AW(X*))^HAW(X_*)+λI)W^-1(X_*)X_*＝(AW(X_*))^HY

因此，最优解应满足

X_*＝W(X_*)((AW(X_*))^HAW(X_*)+λI)^-1(AW(X_*))^HY

利用迭代松弛算法，可得到

可得到初次加权矩阵

及

步骤(4)：迭代中初步估计

在步骤(3)中利用矩阵求逆公式，可得到迭代中X的初步估计

步骤(5)：求取迭代中二次加权矩阵

在步骤(4)中，令

则Q_k可通过以下求解过程获得

因此，

可等价地看作是以下优化问题的最优解：

X_k＝argminJ_k(X),其中

在上式中，

控制噪声的残差最小，

约束重构结果的稀疏性。稀疏性的约束过程本质上是通过X_k与X_k-1对应元素间幅值比确定的，最终使得非零项的幅值越来越大，零项的幅值越来越小。可见，FOCUSS算法实际上等价于在每步迭代中寻找X的最优近似解，可看做将前一个加权最小

范数得到近似解X_k-1，作为后一个迭代过程的初始化输入再通过加权最小

范数过程得到X_k。在文献2(Candes E J,Wakin M B,Boyd SP.Enhancing Sparsity by Reweighted

Minimization.J Fourier Anal Appl,2008,14:877-905.)中，Candes论述了在输入数据中可能为零的位置乘以较大的权值系数，在可能非零的位置乘以较小的权值系数，会令结果的稀疏性更加明显、更加接近原始稀疏数据，但文献中并未给出加权系数的显性表达。

根据上述加权原理，本发明提出了一种具有显性表达式的加权方法，即在第k次迭代中，通过

的值求取能改善估计性能的二次加权矩阵W_k。将

按行取

范数得到的

对FOCUSS方法的加权矩阵赋值方式进行改进，按照

各元素的幅值比重对加权矩阵W_k的对角元素进行赋值，即：

首先令W_k为n×n对角矩阵，然后令对角元素的值满足

其中坐标集

为序列

中前2U个幅度最大值对应的基坐标，U是信源个数。本发明给出了一种权值的显式表达，在其他环境中也能鲁棒性的适应这个表达式，改善系统估计性能，而且经验证其效果很好。

步骤(6)：迭代中二次估计

得到二次加权矩阵后，在第k次迭代中对X进行二次估计，将步骤(4)中的加权矩阵

替换为新加权矩阵W_k，其估计公式为：

其中，A_k＝AW_k。

步骤(7)：验证迭代终止条件

在迭代运算中，若达到收敛条件则迭代终止，其收敛条件为

其中，ε为预设的终止门限。若未满足收敛条件，则转入步骤(3)继续执行迭代算法。

步骤(8)：查找X中的非零行

对得到的最终收敛解X_*，需找出其行稀疏特性。对X_*的行向量分别求取

范数，即

找出

中前U个幅度最大值对应的基坐标T＝[t₁,…,t_U]，就是X_*中非零行的行坐标。

步骤(9)找出目标源的入射角

利用重构信号的行稀疏性与传感矩阵列向量之间的映射关系(即二者相乘时的行列对应关系)，在构造的传感矩阵A的列坐标集中，查找T所在位置对应的角度值

查找原理如图2所示，(a)是稀疏前的，(b)是稀疏后的。传感矩阵A中的第a列，对应X中的第a行。

图3给出的是本发明方法与传统的DOA估计方法(以MUSIC方法为例)单次试验的分辨力对比。其中，空间谱最大值经过了归一化处理，信噪比为10dB，阵元数M＝10，快拍数L＝40。图3(a)是对2个相近非相干信源的DOA估计结果，可发现本发明方法和MUSIC方法都能分辨出两个信源，且本发明方法分辨力更高；图3(b)是对2个相近相干信源的DOA估计结果，可发现MUSIC方法已经失效，而本发明方法仍然能清晰辨别出两个目标。

图4给出的是本发明方法(RewFOCUSS)、RM-FOCUSS算法、

方法、SD-FOCUSS方法应用于DOA估计时，得到的估计角与真实角之间的统计均方根误差(RMSE)曲线图，下表1给出了本发明改进FOCUSS方法(RewFOCUSS)与其他稀疏重构算法DOA估计的平均运行时间对比(单位：秒)

表1

其中，信源为2个相干信源，入射角分别为6.05°和-6.31°，M＝40，L＝100，观测矩阵A中空间角度的网格划分方式为-90°:0.1°:90°，独立蒙特卡洛实验次数为500，仿真软件为Matlab R2012a，计算机处理器为Intel Core i3-4170 3.70GHz，内存为4GB。通过图4可以看出，本发明方法和

方法所得结果的RMSE最小，即估计精度最高。然而由于

方法的运算过程是基于凸优化算法完成的，运算量较大，从表1的各算法完成平均时间对比中可发现，

方法所用时间约是本发明方法的20倍。

基于FOCUSS重新加权算法的DOA估计方法，能够有效克服传统DOA估计不能分辨相干信源的缺陷，通过对稀疏信号精确重构能够高分辨估计出信源的入射角度；与其他稀疏重构方法相比，该方法的估计精度较高，同时能降低硬件运算成本，其综合性能是最高的。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于FOCUSS二次加权算法的DOA估计方法，其特征在于，包括：

步骤一、将传统的DOA估计信号模型，通过对无目标区的填零扩展，构建为DOA估计稀疏信号模型Y＝AX+N；其中，Y为观测量矩阵，A为观测矩阵，X为具有行稀疏性的信号输入矩阵，N为噪声矩阵；

步骤二、基于DOA估计稀疏信号模型，在FOCUSS框架基础上采用二次加权算法进行迭代计算，在每一次迭代过程中进行两次数据重构：

先通过求解最小

范数得到加权矩阵

并对信号输入矩阵X进行初步重构，获得初步重构数据

再利用初步重构数据

进行加权矩阵二次赋值，获得二次加权矩阵W_k，并对输入信号再次重构，得到更加精确的信号重构结果X_k；

步骤三、利用迭代获得的最终收敛解，找到X中的非零行的行坐标，记为T＝[t₁，…，t_U]；从观测矩阵A中找到列号为T＝[t₁，…，t_U]的列对应的角度，即为DOA估计获得的波达方向。

2.如权利要求1所述的方法，其特征在于，所述步骤一为：

步骤1.1：构建传统模型：对于一个具有M个阵元的信号接收阵列，假设U个信源分别以θ₁,…,θ_U的入射角投射到阵列接收机上，则离散后的信号观测模型表达为：

y(l)＝Φs(l)+n(l)＝[a(θ₁),…,a(θ_U)][s(1,l),…,s(U,l)]^T+n(l)

其中，l∈[1,…,L]代表第l个快拍单元，L为总快拍数，y(l)是阵列输出的第l个快拍；s(l)＝s(1,l),…,s(K,l)，s(u,l)为第u个信源的信号幅度；a(θ_u)为第u个信源的信号导引矢量，n(l)为信道中的噪声向量；Φ为观测矩阵，且Φ＝[a(θ₁),…,a(θ_U)]；y(l)、a(θ_u)、n(l)均为M×1的列向量，U<M；

步骤1.2等价填零扩展：对除s(1,l),…,s(U,l)之外的空间中无目标区域进行填零处理，同时扩展对应观测矩阵Φ，得到对应稀疏信号模型y(l)＝Ax(l)+n(l)，其中

为空间方向全角度或区域角度的网格划分，N为划分的格点数量；x(l)＝[x(1,l),…,x(N,l)]^T，x(n,l)|_n＝1,…,N为分别对应每个网格的信源信号幅度；当

等于某一信源真实的入射方向θ_u,u∈[1,…U]时，x(n,l)的值不为0，否则其值为0，即x(n,l)的值为：

步骤1.3构建多快拍稀疏信号模型为：

Y＝AX+N，

其中Y＝{y(l)|_l＝1,…,L}，X＝{x(l)|_l＝1,…,L}，及N＝{n(l)|_l＝1,…,L}，Y为M×L的矩阵，A为M×N的矩阵，X为N×L的矩阵，N为M×L的矩阵。

3.如权利要求2所述的方法，其特征在于，步骤二利用松弛迭代方法得到所述初步加权矩阵

其中，

x_k-1(n,l)表示第k-1次迭代求出的x(n,l)的估计结果，而后得到初步重构数据

其中

λ为噪声因子，p为收敛因子，I为单位对角阵；

利用初步估计结果

再进行加权矩阵二次赋值，获得二次加权矩阵W_k，其方法为：

对

的行向量分别求取

范数，获得

取序列

中前2U个幅度最大值对应的基坐标，组成坐标集U；

令W_k为N×N的对角矩阵，对角元素ω_n的取值为：若n属于坐标集U，则ω_n的幅度为1；否则，根据

幅度的大小设置ω_n，

幅度越大，ω_n的幅度值越大，ω_n幅度值取值范围为[0,1)；

然后可得到迭代中二次重构数据X_k

其中，A_k＝AW_k。

4.如权利要求3所述的方法，其特征在于，对角元素ω_n的取值为：

5.如权利要求1所述的方法，其特征在于，FOCUSS迭代运算中，X的初始值的确定方式为：X₀＝A^H(AA^H)^-1Y。

6.如权利要求1所述的方法，其特征在于，FOCUSS迭代运算所用到的噪声因子λ的求取方式为：

其中σ为已知的观测噪声方差；β为贝叶斯广义方差，

Γ(·)为标准gamma分布函数，p为收敛因子。

7.如权利要求3或6所述的方法，其特征在于，FOCUSS迭代运算所用到的收敛因子p的求取方式为：根据估计结果的性能要求和允许的运算量设置p，p为(0,1]之间任意的常数，p的值越小，估计结果的稀疏性就越好，运算量随之增大；p的值越大，估计结果的稀疏性就越差，运算量随之降低。

8.如权利要求1或6所述的方法，其特征在于，所述步骤三为：

步骤3.1查找X中的非零行：对迭代获得的最终收敛解的行向量分别求取

范数，获得序列

找出其中前U个最大值对应的基坐标T＝[t₁,…,t_U]，就是X中非零行的行坐标；

步骤3.2找出目标源的入射角：在传感矩阵A的列向量中，分别查找第t₁,…,t_U列所在位置对应的角度值

即为所估计的目标源入射角。

Images