CN109307855A - 基于网格误差模型的无网格稀疏近似最小方差doa估计方法 - Google Patents

Abstract

本发明涉及一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，建立基于网格误差模型的阵列接收信号模型，并基于该模型依次对信号、噪声功率和网格误差进行迭代求解，使稀疏近似最小方差DOA估计方法的定位精度不再受空间网格划分精度的限制，当目标信号的DOA与网格点不匹配时仍能实现较高的DOA估计精度。

Description

基于网格误差模型的无网格稀疏近似最小方差DOA估计方法

技术领域

本发明属于信号处理等领域，涉及一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，通过对网格与信号方位之间误差的估计，避免网格划分与信号方位不匹配而引起的DOA估计误差。

背景技术

阵列信号处理在雷达、声纳等领域上具有广泛的应用，目标方位(Direction ofarrival，DOA)估计是阵列信号处理一大主要任务。稀疏信号处理类的DOA估计算法是近十年来发展起来的DOA估计算法，与常规的DOA估计算法相比，该类DOA估计算法能够适用于小快拍和低信噪比的情况下，性能远远优于常规的DOA估计算法。

稀疏信号处理算法主要可以分为正则参数类算法和非正则参数类算法。正则参数类算法如稀疏谱估计算法和l₁-SVD算法使用l₁范数作为惩罚项，通过正则参数将惩罚项与模型拟合误差相结合，实现稀疏DOA估计。正则参数控制着结果稀疏性和模型拟合误差之间的平衡，随着正则参数的增加，DOA估计结果的稀疏性逐渐增大，模型拟合误差也逐渐增加。在实际数据处理时正则参数的选取通常较为困难，而不合适的正则参数将导致该类算法的性能大大下降，这使得该类算法在使用时受到了较大的限制。稀疏近似最小方差(Sparseasymptotic minimum variance，SAMV)算法(H.Abeida，Q.Zhang，J.Li，et al.Iterativesparse asymptotic minimum variance based approaches for array processing[J].Transactions on Signal Processing，2013，61(4)：933-944)则是一种常见的非正则参数类稀疏DOA估计算法，该算法利用近似最小方差准则给出信号和噪声功率的迭代关系式，通过迭代的方式使协方差矩阵模型接近采样协方差矩阵。算法的整个求解过程中仅需提供迭代停止门限，而无需任何类似正则参数的超参数，因此相较于正则参数类的稀疏信号处理算法来说，该算法在实际信号处理中更易被实现。

SAMV算法在进行DOA估计前，需将观测空间划分为离散的网格，在该网格上进行DOA估计。该算法存在一大缺陷，即当信号的真实方位与所划分的网格点不一致时，信号方位的估计值为距离该信号最近的网格点所对应的方位角，这导致该信号方位角的估计值与真实值之间存在一个恒定的误差。虽然增加空间网格的划分精度在一定程度上能够减轻目标方位与网格点不匹配的问题，但却大大增加了运算量。因此，需选择合适的方式对网格点与目标信号方位角之间的误差进行估计，使得SAMV算法的定位精度不受网格划分的限制。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，使SAMV算法的DOA估计精度不受网格划分的限制，通过对网格与目标方位之间误差的估计解决SAMV算法由于信号方位与网格点不匹配造成的DOA估计不精准的问题。

技术方案

一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，其特征在于估计步骤如下：

步骤1：采用阵元间距为半波长的M元均匀线列阵接收窄带信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n)，1≤n≤N，i＝1，...，M；

将观测空间[-90°，90°]划分为Q个网格，所述90°为端射方向，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁，θ₂，...，θ_Q]，在该网格上，阵列的接收信号模型表示为x(n)＝Φ(β)s(n)+e(n)，n＝1，...，N

其中：和分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；Φ(β)＝A(Θ)+B(Θ)diag(β)，所述的为阵列流形矩阵，a(θ)＝[1 e^{-j2πdsin(θ)f/c} ... e^{-j2πd(M-1)sin(θ)f/c}]^T，f为窄带信号的中心频率，c为声速；所述的B(Θ)＝[a′(q₁)，...，a′(q_Q)]，a′(q)表示a(q)在θ处的导数；所述的β＝[β₁，...，β_Q]^T为网格误差，将与空间信号距离最近的网格点组成的集合记为对于属于的网格点，对应的网格误差为真实的信号方位与网格点所代表方位之差，其余网格点对应的网格误差为0；diag(·)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵；

假设环境噪声为均匀白噪声情况，则阵列接收信号的采样协方差矩阵计算为R＝E[x(n)x^H(n)]＝Φ(β)R_sΦ^H(β)+s²I_M

其中：E[·]为期望算子；R_s为信号协方差矩阵；上标“H”为共轭转置符号；σ²为信号功率；I_M代表M维的单位矩阵；一般地，协方差矩阵由采样协方差矩阵所代替；

对协方差矩阵两端进行向量化得

其中：vec(·)为矩阵向量化算子；矩阵的第i列为其中Φ_i为Φ(β)的第i列，上标“*”表示求共轭，表示Kronecker积；p＝diag(R_s)为R_s对角线元素组成的向量；

步骤2：采用迭代的方式计算信号和噪声功率：

其中：上标(i)表示第i次迭代；表示第i次迭代得到的第q个网格点对应的方位角； 所述Θ⁽ⁱ⁾为第i次迭代所求网格点集合，所述Tr(·)为矩阵求迹符号；

迭代的初始值为：

其中：||·||_F表示Frobenius范数；

步骤3：在第i次迭代中，完成信号和噪声计算后，计算新的协方差矩阵估计量为

对中最大的(M-1)个峰值所对应的网格进行网格误差求解。采用多项式求根的方式对网格误差进行求解，所述求根方程为

其中： Re[·]为复数求实部；

由于求解网格误差多项式为3次多项式，因此存在3个根，故按以下准则对所求的根进行筛选：

1)网格误差必为实数；

2)所求的网格误差需满足：

3)若根据前两个准则筛选完毕后，仍存在两个及以上的根，则选择绝对值最小的根作为该次迭代的网格误差量；

当计算完所有网格误差后，新的网格点计算为：

依次对信号功率、噪声功率以及网格误差进行迭代计算，当所计算的网格误差的l₂范数小于所选门限值ρ时，停止对网格误差的迭代计算，仅对信号功率和噪声进行求解；当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于所选的迭代停止门限值ε时，迭代终止，所求的信号功率的峰值所对应的方位角即为所估计的信号的DOA。

有益效果

本发明提出的一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，建立基于网格误差模型的阵列接收信号模型，并基于该模型依次对信号、噪声功率和网格误差进行迭代求解，使稀疏近似最小方差DOA估计方法的定位精度不再受空间网格划分精度的限制，当目标信号的DOA与网格点不匹配时仍能实现较高的DOA估计精度。

本发明的有益效果是：建立基于网格误差模型的阵列接收信号模型，并将其应用于近似最小方差准则中，依次对信号、噪声功率和网格误差进行求导并置零，得到信号、噪声功率的迭代关系式和网格误差的求解方程，使SAMV算法的定位精度不再受空间网格划分精度的限制，当目标信号方位与网格点不匹配时仍可实现较高的定位精度。

附图说明

图1：为采用OGSAMV算法进行信号DOA估计的总体流程图，图中上标(new)表示最新迭代得到的估计值，(old)表示上一次迭代的估计值；

图2：(a)OGSAMV算法和(b)SAMV算法的DOA估计结果，其中红色圆点为所估计的信号方位，蓝色虚线为真实的信号方位。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案包括以下步骤：

1)建立基于网格误差模型的阵列接收信号模型

采用阵元间距为半波长的M元均匀线列阵作为接收阵列，接收单频信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n)，0≤n≤N，i＝1，...，M。将观测空间[-90°，90°](其中90°为端射方向)划分为Q个网格，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁，θ₂，...，θ_Q]。当信号分布在该离散网格上时，阵列的接收信号可表示为x(n)＝A(Θ)s(n)+e(n)，n＝1，...，N，其中和分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；为阵列流形矩阵，对于均匀线列阵来说，a(θ)＝[1 e^{-j2πdsin(θ)f/c} ... e^{-j2πd(M-1)sin(θ)f/c}]^T，f为信号的中心频率，d为阵元间距，c为声速。

假设空间中实际存在的信号个数为K个，它们的方位角组成的向量表示为当第k个信号的方位角不在所划分的网格点上时，利用上述阵列接收信号模型进行DOA估计的结果为距离该信号最近的网格点所对应的方位角导致了该信号方位的估计值与真实值之间存在一定的误差。为避免该误差的产生，将阵列流形在上进行泰勒展开，并保留至一阶项，则信号真实方位角所对应的阵列流形由下式估计：

其中表示在处的导数，为信号方位和与之距离最近的网格点之间的偏差。对于任意在已划分的空间网格中总可以找到与之距离最近的网格点对应的网格误差记为将这些网格点组成的向量记为对于网格点所对应的网格点偏差则为β_q＝0。将所有网格点对应的网格偏差记为β＝[β₁，...，β_Q]^T，对离散网格修正后的阵列流形矩阵可表示为

Φ(β)＝A(Θ)+B(Θ)diag(β) (2)

其中B(Θ)＝[a′(θ₁)，...，a′(θ_Q)]，diag(g)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵。建立基于式(2)所示的网格误差模型的阵列接收信号模型：

x(n)＝Φ(β)s(n)+e(n)，n＝1，...，N (3)

根据式(3)，计算阵列接收信号的协方差矩阵：

R＝E[x(n)x^H(n)]＝Φ(β)R_sΦ^H(β)+R_e (4)

其中R_s和R_e分别为信号协方差矩阵和噪声协方差矩阵，上标“H”为共轭转置；E[·]为求期望算子。考虑噪声为均匀高斯白噪声时，式(4)的协方差矩阵重新表示为：

R＝Φ(β)R_sΦ^H(β)+σ²I_M (5)

其中σ²为噪声功率。对式(5)两端进行向量化，可得

vec(·)为矩阵向量化算子；矩阵的第i列为其中Φ_i为Φ(β)的第i列，上标“*”表示求共轭，表示Kronecker积；p_s＝diag(R_s)，p_e＝σ²vec(I_M)。

定义一个新矩阵和新向量p＝[p_s ^T σ²]^T，则式(6)可重新表示为：

r(p)＝Sp (7)

2)信号功率和噪声功率的求解

根据近似最小方差准则，最优的p的估计量可由下式得到：

其中为采样协方差矩阵；

令r(p)＝r′_q+p_qs_q，其中r′_q为r(p)除去第q个信号分量后的向量，s_q为矩阵S的第q列，将其带入至式(8)的f(p)中可得：

对p_q求导并置零，可得信号和噪声估计值：

Tr(·)为矩阵求迹。由于信号和噪声功率的求解须已知p_q、σ²和R，故采用迭代的方式进行求解。迭代式如下：

其中上标(i)表示第i次迭代；是第i次迭代后的协方差矩阵估计量，

为防止迭代中出现信号功率和噪声功率出现小于0的情况，令和σ²＝Tr(R^-1)/Tr(R^-2)，带入至式(11)中可得

由于Φ_q中存在未知变量β_q，因此在第(i+1)次迭代中计算信号和噪声功率时可固定即考虑到为第i次迭代第q个网格点对应的阵列流形，同时还表示为 表示第i次迭代得到的第q个网格点对应的方位角，将其带入至式(12)中得：

由于第i次迭代得到的阵列流形矩阵可表示为因此式(13)中也可由所计算。信号和噪声功率的迭代初始值可由和给出，其中||·||_F为Frobenius范数。

3)网格误差计算

在第(i+1)次迭代中，按式(13)得到信号与噪声功率后，此时的协方差矩阵估计量为将信号和噪声功率带入至式(9)中得在中仅存在β_q一个未知变量。为求得合适的β_q以使最小，令得到方程：

其中 Re[·]为复数求实部。

对于方程(14)来说，最多存在3个不相同的根，故需要一个合适的选取准则以选取一个最优的根作为该次迭代所估计的网格误差量。

首先从网格误差的计算式可以看出，所求的网格误差必为实数，因此需排除复数根。其次，为防止所求的网格误差b_q的绝对值过大，使得第q个网格点大于第(q+1)个网格点或者小于第(q-1)个网格点，所得到的b_q必须满足如下关系式：

根据上述两条准则舍去不满足条件的根，若仍存在2个及以上的根，则选择绝对值最小的根作为该次迭代的网格误差量，即选择距离该网格点最近的网格点作为新的网格点。为加快OGSAMV算法的计算速度，并不需要对所有网格点进行误差求解。考虑到M元的均匀线列阵最多可分辨的(M-1)个信号，因此只需对中最大的(M-1)个峰值对应的网格进行误差求解即可。

当计算完所有的网格误差后，新的网格点可由得到。当前后两次迭代所得的网格点满足||Θ⁽ⁱ⁺¹⁾-Θ⁽ⁱ⁾||₂≤ρ时，其中r为网格误差迭代停止门限，即可停止对网格误差的求解，而仅对信号功率进行迭代求解；当前后两次迭代所计算的信号和噪声功率满足其中e为功率计算的迭代停止门限，则完成了对目标信号的DOA估计。

具体实施如下：

1)采用阵元间距为半波长的M元均匀线列阵作为接收阵列，接收单频信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n)，0≤n≤N，i＝1，...，M。将空间[-90°，90°](其中90°为端射方向)划分为Q个网格，各网格点所代表的方向角组成的向量记为Θ＝[θ₁，θ₂，...，θ_Q]。假设空间中实际存在的信号个数为K，它们的方位角组成的向量表示为当第k个信号的方位角不在所划分的网格点上时，寻找距离该信号最近的网格点所对应的方位角将阵列流形在上进行泰勒展开，并保留至一阶项，则信号真实方位角所对应的阵列流形可由计算，其中表示在处的导数，为网格点与信号的偏差。对于任意在已划分的空间网格中总可以找到与之距离最近的网格点所对应的网格误差为将这些网格点组成的向量记为对于网格点所对应的网格点偏差为β_q＝0。将所有网格点对应的网格偏差组成的向量记为β＝[β₁，...，β_Q]^T，则对离散网格修正后的阵列流形矩阵可表示为Φ(β)＝A(Θ)+B(Θ)diag(β)，其中B(Θ)＝[a′(θ₁)，...，a′(θ_Q)]，diag(·)表示由矩阵主对角线元素组成的向量或以向量元素为主对角线的对角矩阵。

基于网格误差模型的阵列接收信号模型可表示为x(n)＝Φ(β)s(n)+e(n)，n＝1，...，N，而所对应的阵列接收信号的协方差矩阵模型为：R＝E[x(n)x^H(n)]＝Φ(β)R_sΦ^H(β)+R_e，其中R_s和R_e分别为信号协方差矩阵和噪声协方差矩阵，上标“H”为共轭转置。考虑噪声为均匀高斯白噪声时，协方差矩阵模型可重新写为：R＝Φ(β)R_sΦ^H(β)+σ²I_M，其中s²为噪声功率。对等式两端进行向量化，可得其中vec(·)为矩阵向量化算子；矩阵的第i列为其中Φ_i为Φ(β)的第i列，上标“*”表示求共轭，表示Kronecker积；p_s＝diag(R_s)，p_e＝σ²vec(I_M)。定义一个新矩阵和新向量p＝[p_s ^T σ²]^T，则协方差矩阵模型可重新表示为r(p)＝Sp。

2)根据近似最小方差准则，最优的p的估计量可由下式得到：其中 为采样协方差矩阵。

令r(p)＝r′_q+p_qs_q，r′_q为r(p)除去第q个信号分量后的向量，s_q为矩阵S的第q列，f(p)可写为：

对p_q求导并置零，可得信号和噪声的迭代关系式：

其中上标(i)表示第i次迭代；是第i次迭代后的协方差矩阵估计量，由于Φ_q中存在未知变量b_q，因此在第(i+1)次迭代中可令进行求解，即考虑到为第i次迭代第q个网格点对应的阵列流形，故同时可表示为 表示第i次迭代得到的第q个网格点对应的方位角，则信号和噪声的迭代式可重新写为：

由于第i次迭代得到的阵列流形矩阵可表示为因此协方差矩阵也可由式所计算。信号和噪声功率的迭代初始值可由和给出。

3)在第(i+1)次迭代中，当得到信号和噪声功率后，协方差矩阵估计量可计算为为求得合适的β_q以使f(p_q)最小，令并使得到关于β_q的方程，求解方程即可完成对网格误差的计算。为加快OGSAMV算法的计算速度，并不需要对所有网格点进行误差求解。考虑到M元的均匀线列阵最多可分辨(M-1)个信号，因此只需对中最大的(M-1)个峰值对应的网格进行误差求解即可。

从网格误差的求解公式中可以看出，所求得的β_q必须满足为实数。同时为防止所求的网格误差β_q过大，使得第q个网格点大于第(q+1)个网格点或者小于第(q-1)个网格点，所得到的β_q必须满足如下关系式：

舍弃不满足上述条件的根，若仍存在2个及以上的根，则选择绝对值最小的根作为该网格点在该次迭代中的网格误差量。

4)当计算完网格点误差后，新的网格点可由得到。当前后两次迭代所得的网格点之差满足||Θ⁽ⁱ⁺¹⁾-Θ⁽ⁱ⁾||₂≤ρ时，其中ρ为网格误差计算的停止门限，即可停止对网格误差的求解，而仅对信号功率进行迭代求解。当前后两次迭代所计算的信号和噪声功率满足其中e为功率迭代停止门限，则完成了对目标信号的DOA估计。

从仿真结果中可以看出，SAMV算法所估计出的两个信号的方位分别为19°和24°，相应的估计误差为而OGSAMV算法所估计出的两个信号的方位分别为19.44°和24.31°，相应的估计误差为由此可以看出，OGSAMV算法的定位精度远远高于SAMV算法。

Claims (1)

1.一种基于网格误差模型的无网格稀疏近似最小方差DOA估计方法，其特征在于估计步骤如下：

步骤1：采用阵元间距为半波长的M元均匀线列阵接收窄带信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),1≤n≤N,i＝1,...,M；

将观测空间[-90°,90°]划分为Q个网格，所述90°为端射方向，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]，在该网格上，阵列的接收信号模型表示为x(n)＝Φ(β)s(n)+e(n),n＝1,...,N

其中： 和分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；Φ(β)＝A(Θ)+B(Θ)diag(β)，所述的为阵列流形矩阵，a(θ)＝[1 e^{-j2πdsin(θ)f/c} ... e^{-j2πd(M-1)sin(θ)f/c}]^T，f为窄带信号的中心频率，c为声速；所述的B(Θ)＝[a'(q₁),...,a'(q_Q)]，a'(q)表示a(q)在θ处的导数；所述的β＝[β₁,...,β_Q]^T为网格误差，将与空间信号距离最近的网格点组成的集合记为对于属于的网格点，对应的网格误差为真实的信号方位与网格点所代表方位之差，其余网格点对应的网格误差为0；diag(·)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵；

假设环境噪声为均匀白噪声情况，则阵列接收信号的采样协方差矩阵计算为R＝E[x(n)x^H(n)]＝Φ(β)R_sΦ^H(β)+s²I_M

其中：E[·]为期望算子；R_s为信号协方差矩阵；上标“H”为共轭转置符号；σ²为信号功率；I_M代表M维的单位矩阵；一般地，协方差矩阵由采样协方差矩阵所代替；

对协方差矩阵两端进行向量化得

其中：vec(·)为矩阵向量化算子；矩阵的第i列为其中Φ_i为Φ(β)的第i列，上标“*”表示求共轭，表示Kronecker积；p＝diag(R_s)为R_s对角线元素组成的向量；

步骤2：采用迭代的方式计算信号和噪声功率：

其中：上标(i)表示第i次迭代；表示第i次迭代得到的第q个网格点对应的方位角；所述Θ⁽ⁱ⁾为第i次迭代所求网格点集合，所述Tr(·)为矩阵求迹符号；

迭代的初始值为：和

其中：||·||_F表示Frobenius范数；

步骤3：在第i次迭代中，完成信号和噪声计算后，计算新的协方差矩阵估计量为

对中最大的(M-1)个峰值所对应的网格进行网格误差求解。采用多项式求根的方式对网格误差进行求解，所述求根方程为

其中：c₁＝γ²，c₂＝3βγ, Re[·]为复数求实部；

由于求解网格误差多项式为3次多项式，因此存在3个根，故按以下准则对所求的根进行筛选：

1)网格误差必为实数；

2)所求的网格误差需满足：

3)若根据前两个准则筛选完毕后，仍存在两个及以上的根，则选择绝对值最小的根作为该次迭代的网格误差量；

当计算完所有网格误差后，新的网格点计算为：

依次对信号功率、噪声功率以及网格误差进行迭代计算，当所计算的网格误差的l₂范数小于所选门限值ρ时，停止对网格误差的迭代计算，仅对信号功率和噪声进行求解；当前后两次迭代得到的信号与噪声功率组成的向量差值的l₂范数小于所选的迭代停止门限值ε时，迭代终止，所求的信号功率的峰值所对应的方位角即为所估计的信号的DOA。