CN109783960A - 一种基于网格部分细化的波达方向估计方法 - Google Patents

Abstract

本发明公开一种基于网格部分细化的波达方向估计方法，通过学习和裂变过程，裂变过程通过产生新格点对网格进行细化，学习过程不断逼近波达方向，在真实的波达方向附近网格被细划分，而在远离波达方向的网格区域被粗划分，这就实现了网格的部分细化，不但保证了估计精度，而且与之前的离格DOA估计算法相比，网格点数大大减少，计算量随之减小。本发明不需要将信源数作为先验，以尽量稀疏的网格划分，减少了网格数，从而减小计算复杂度，算法耗时少；在非常稀疏的初始格点划分的条件下，通过最小间隔阈值自定义，保证了算法的估计精度。

Description

一种基于网格部分细化的波达方向估计方法

技术领域

本发明涉及阵列信号处理技术领域，具体涉及一种基于网格部分细化的波达方向估计方法。

背景技术

波达方向(Direction-of-arrival,DOA)估计技术在雷达、声纳、气象等众多领域都有着广泛应用。传统的波达方向估计技术中最为众所周知的就是MUSIC和ESPRIT等算法，但这些传统算法的性能受到诸多因素的制约，比如，高信噪比和大快拍数。压缩感知理论的提出使得DOA估计问题有了重大突破，尽管无论是基于网格模型还是离格模型的DOA估计算法都比传统算法有一定优势，但这些稀疏信号重建模型都是将空间角度网格等间距划分，为追求估计精度，将网格密集划分，从而导致算法计算量巨大。

发明内容

本发明现有技术在DOA估计中计算复杂度和估计精度方面的不足，提供一种基于网格部分细化的波达方向估计方法。

为解决上述问题，本发明是通过以下技术方案实现的：

一种基于网格部分细化的波达方向估计方法，包括以下步骤：

步骤1、基于初始划分的网格点，构建当前观测数据的初始离格波达方向估计模型；

步骤2、利用基于根稀疏贝叶斯学习的离格波达方向估计算法，对当前离格波达方向估计模型进行一次学习，在学习过程中，通过更新网格的划分，以更新离格波达方向估计模型；

步骤3、判断当前离格波达方向估计模型是否满足学习停止条件，即或当前次数是否达到预设的次数阈值；其中α^j表示当次学习后的信号的先验参数，α^j-1表示当次学习前的信号的先验参数，τ表示预设的迭代停止阈值；如果满足，则当前离格波达方向估计模型中的后验概率的均值即为所求的结果；否则，转入步骤4；

步骤4、计算当前离格波达方向估计模型中每个网格点位置信号的平均功率，并从中选出平均功率较大的M-1个网格点来作为待裂变的网格点；其中M是均匀天线阵列中阵元个数；

步骤5、当所选出的待裂变的网格点中，存在至少一个待裂变的网格点的左右间隔同时大于预设的间隔阈值，即存在至少一个需要裂变的网格点时，则转入步骤6；否则，返回步骤2；

步骤6、对当前离格波达方向估计模型进行一次裂变，在裂变过程中，在每一个需要裂变的网格点的前一个网格点和后一个网格点的中点位置各插入一个新的网格点，以更新离格波达方向估计模型；

步骤7、返回步骤2。

上述方案中，所述预设的迭代停止阈值τ的取值范围在10^-4～10^-3之间。

上述方案中，所述预设的次数阈值的取值范围在500～1000之间。

上述方案中，所述预设的间隔阈值的取值范围在1°～3°之间。

与现有技术相比，本发明具有如下特点：

(1)通过学习和裂变过程，裂变过程通过产生新格点对网格进行细化，学习过程不断逼近波达方向，在真实的波达方向附近网格被细划分，而在远离波达方向的网格区域被粗划分，这就实现了网格的部分细化，不但保证了估计精度，而且与之前的离格DOA估计算法相比，网格点数大大减少，计算量随之减小。

(2)本发明直接对所有待裂变网格点中，左右间隔同时大于预设的网格间隔阈值的网格点进行裂变，利用裂变过程，空间角度网格被选择性地进行细化，减小了网格点的总数，从而降低计算复杂度；

(3)本发明通过将裂变过程中格点间隔阈值设置为自定义的值，保证了算法的估计精度；

(4)本发明采用稀疏贝叶斯学习模型，不需要事先知道信源数，更加适用于应用到实际环境中。

附图说明

图1为一种基于网格部分细化的波达方向估计方法的流程图。

图2是信号源的空间谱示意图。

图3是本发明裂变过程的示意图。

图4是本发明所提算法和现有算法的空间谱示意图。

图5是本发明所提算法与现有算法均方根误差(RMSE)对比曲线图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实例，并参照附图，对本发明进一步详细说明。

一种基于网格部分细化的波达方向估计方法，如图1所示，其具体包括如下步骤：

步骤1：参数初始化：

基于初始划分的网格点，构建当前观测数据的离格波达方向(DOA)估计模型。

为了将DOA估计转化为稀疏重建问题，一般令为在角度空间[-90°,90°]上的等间隔划分后的网格点，其中N表示网格点的数目。

对于网格失配的情况，对导向矢量采用泰勒展开式方法计算，即：

其中，表示离θ_k最近的网格点，令Φ(β)＝A+Bdiag(β),其中，β是除了第nk个元素为其余元素都为0的N维向量。

则观测数据模型可以重新写成：

y(t)＝Φ(β)x(t)+e(t),t＝1,2,…,T

将上式写成矩阵的形式，得到离格DOA估计模型：

Y＝Φ(β)X+E

其中，Y＝[y(1),y(2),…,y(T)],X＝[x(1),x(2),…,x(T)],E＝[e(1),e(2),…,e(T)]。

步骤2：学习过程：

通过已知的OGRSBL(基于根稀疏贝叶斯学习的离格DOA估计)算法对当前观测数据的离格波达方向估计模型进行一次学习。

上述OGRSBL算法如下：

(1)稀疏贝叶斯模型：

首先假定噪声矢量服从高斯分布，则

其中，α₀＝σ^-2，σ²表示噪声方差。

得到阵列接收数据矩阵的后验概率分布为：

假设噪声方差未知，由于Gamma先验是高斯分布的共轭先验，故假定α₀服从Gamma先验分布，即：

p(α₀；c,d)＝Γ(α₀|c,d)

其中，Gamma分布满足

传统的稀疏贝叶斯模型假定信号矢量X的第i行服从方差为α_i的高斯分布，令α＝[α₁,α₂,…,α_N]^T,Δ＝diag(α)，则有

为了实现分层先验，需要定义超参数α，同样α的分布选择Gamma分布：

其中，ρ一般为非常小的正的常数。

假定β服从均匀分布，即

联合概率分布为：

p(X,Y,α,α₀,β)＝p(Y|X,α₀,β)p(X|α)p(α)p(α₀)p(β)

(2)贝叶斯推理：

参数的后验分布p(X,α,α₀,β|Y)不能被直接求解出来，EM算法可以用来实现贝叶斯推理。首先，很容易得到X的后验分布：

均值μ(t)和方差Σ分别为：

μ(t)＝α₀ΣΦ^Hy(t),t＝1,2,…,T

Σ＝(α₀Φ^HΦ+Δ^-1)^-1

计算μ(t)和Σ需要知道α₀、α和β的值。根据其他参考文献可以得出参数α和α₀的更新为：

其中，Ξ_t＝μ(t)(μ(t))^H+Σ，tr(·)表示矩阵的迹。

针对的更新，采用多项式求根的方法，对联合概率分布取对数，然后求关于的导数，并令导数为零得到：

其中，z_i为多项式所求根中最接近单位圆的点。

假设X第i行能量比较小，对应的就可以忽略，所以在每次迭代中不必更新每一个网格点，而由M个阵元组成的阵列可以分辨的最大目标数为M-1，因此选择M-1个能量最大的格点作为可调节的网格点，X第i行的平均功率为：

其中，μ_ti可以由μ(t)＝α₀ΣΦ^Hy(t)求得。

学习过程就是根据上述过程，每次学习，更新角度空间集合从而更新网格的划分，同时更新后验概率的均值μ(t)，后验概率的方差Σ，信号的先验参数α和噪声的先验参数α₀，使新的网格点更加逼近真实的波达方向。

步骤3、学习停止判断：

判断当次学习后的离格波达方向估计模型是否满足学习停止条件，即：或者当前迭代次数达到预设的最大值；其中α^j表示当次学习后的信号的先验参数，α^j-1表示当次学习前的信号的先验参数，τ表示预设的迭代停止阈值：

如果满足，则当前离格波达方向估计模型中的后验概率的均值μ(t)即为所求的结果；

若不满足，转入步骤4。

步骤4：确定待裂变的网格点：

根据公式计算平均功率来选择待裂变的网格点。平均功率P(i)越大，说明第i个位置附近存在真实波达方向的可能性越大。而由M个阵元组成的阵列可以分辨的最大源目标数为M-1，故选择前M-1个最大的P(i)对应的角度θ作为待裂变的网格点。但当某一信号源能量较大时，使其旁瓣能量可能大于其他信源的主瓣能量，如图2所示，i₁位置的信源产生的旁瓣能量大于i₂位置信源的主瓣能量，这样在选择裂变网格点时，可能会丢失第i₂个网格点，所以在前M-1个最大的平均功率P中选择局部最大值对应的角度位置作为待裂变的网格点。

步骤5、裂变进入判断：

当所选出的待裂变的网格点中，存在至少一个需要裂变的网格点(待裂变的网格点的左右间隔同时大于预设的网格间隔阈值)时，也即r(i-1)>θ₀且r(i)>θ₀时，需要进行裂变过程，继续步骤6，进行裂变，其中θ₀为用户自定义的值，与波达方向的估计精度有关，网格点之间的间隔定义为：相反，如果所有待裂变的网格点均满足r(i-1)≤θ₀或r(i)≤θ₀，跳过裂变过程，则返回步骤2，继续学习。

步骤6、裂变过程：

通过裂变过程加入新的网格点，选择性细化角度空间[-90°,90°]。

首先，不同于现有技术，在确定新增网格点位置时需要复杂的判断过程，且每个网格点每次裂变新增一个网格点的方式，本发明直接对所有待裂变网格点中，左右间隔同时大于预设的网格间隔阈值，也即r(i-1)>θ₀且r(i)>θ₀的网格点，即需要裂变的网格点进行裂变。对于需要裂变的每个网格点，在每一个网格点与其前一个和后一个网格点的中点位置各插入一个新的网格点。比如选定裂变点为第i个网格点，那么裂变后新的网格点的位置为第i个格点左右间隔(即r(i-1)和r(i))的中点。如图3所示，用其中一个裂变格点作为例子说明，第一次裂变，A点分裂成A₁和A₂，第二次裂变时，A₁点的平均功率比较大，故其产生分裂，分裂成A₁₁和A₁₂，其中A₁、A₂、A₁₁和A₁₂均为其相邻两个格点的中点。当然分裂点A和A₁同时通过学习过程不断向θ₁逼近。通过不断地学习和迭代，网格点会逐渐向真实波达方向逼近。

然后，更新X和根据上述裂变新增网格点的过程，增加新的更新角度空间集合根据公式X的先验分布是α的函数，所以更新X就是更新α。而由图2可知，一个网格点分裂后总共新增了两个网格点，比如B分裂后产生B₁,B₂,C三个点，那么根据分裂前后的分布分别为：

为了保证分裂前后X的分布保持不变，那么令即在这里取

步骤7、完成网格裂变之后，直接返回步骤2，对当前观测数据的离格波达方向估计模型进行继续学习。

通过上述的学习和裂变过程，裂变过程通过产生新格点对网格进行细化，学习过程不断逼近波达方向，结果正如图2所示，在真实的波达方向附近网格被细划分，而在远离波达方向的网格区域被粗划分，这就实现了网格的部分细化，不但保证了估计精度，而且与之前的离格DOA估计算法相比，网格点数大大减少，计算量随之减小。

仿真条件：假设均匀线阵由M＝10个阵元组成，阵元间隔为d＝λ/2，快拍数T＝30，设置ρ＝0.01，c＝d＝1×10^-4。令α₀＝1，α中所有元素均为1，学习过程的迭代停止条件中的阈值τ＝1×10^-4，最大迭代次数为500。本发明算法的初始网格点间隔为r₀＝180°/(M-1)＝20°，裂变终止条件中的最小网格点间隔为2°。

图4展示了本发明算法与现有算法的空间谱图。两个窄带远场信号源的来波方向分别是-17.4°和13.7°，信噪比SNR＝10dB。从图4可以看出，相比其他算法，本发明提出的算法空间谱更窄，而且更加接近来波方向，说明该算法的角度分辨率更高。

图5展现了本发明算法与现有算法的均方根误差随信噪比SNR变化的情形。考虑K＝2个信号源θ₁和θ₂分别来自间隔[-10°,0°]和[20°,30°]，对于每一个SNR都进行R＝500次蒙特卡洛试验，均方根误差定义为：其中为第k个信号源的估计值。均方根误差的数值对应算法的优劣，从图5可以看出，所提算法在大部分SNR下都比其他算法的均方根误差小，说明本发明所提算法比现有算法的估计精度高。

表1不同格点间隔下的算法耗时

表1呈现了各种算法在信噪比SNR＝10dB不同格点间隔下的CPU平均耗时。现有l₁-SVD算法和现有OGRSBL算法的初始格点间隔依次取[1°,2°,4°,6°,8°]，现有GEDOA算法的初始格点间隔为20°，与现有l₁-SVD算法和现有OGRSBL算法相比，由于本发明PRGDOA算法的网格划分比较粗糙，所以CPU耗时较短。虽然本发明PRGDOA算法比现有GEDOA算法的CPU耗时稍微多一点，因为本发明PRGDOA算法比现有GEDOA算法的格点数多一些，但同样条件下本发明PRGDOA算法的估计精度更高。

本发明不需要将信源数作为先验，以尽量稀疏的网格划分，减少了网格数，从而减小计算复杂度，算法耗时少；在非常稀疏的初始格点划分的条件下，通过最小间隔阈值自定义，保证了算法的估计精度。

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。

Claims (4)

1.一种基于网格部分细化的波达方向估计方法，其特征是，包括以下步骤：

步骤1、基于初始划分的网格点，构建当前观测数据的初始离格波达方向估计模型；

步骤2、利用基于根稀疏贝叶斯学习的离格波达方向估计算法，对当前离格波达方向估计模型进行一次学习，在学习过程中，通过更新网格的划分，以更新离格波达方向估计模型；

步骤3、判断当前离格波达方向估计模型是否满足学习停止条件，即或当前次数是否达到预设的次数阈值；其中α^j表示当次学习后的信号的先验参数，α^j-1表示当次学习前的信号的先验参数，τ表示预设的迭代停止阈值；如果满足，则当前离格波达方向估计模型中的后验概率的均值即为所求的结果；否则，转入步骤4；

步骤4、计算当前离格波达方向估计模型中每个网格点位置信号的平均功率，并从中选出平均功率较大的M-1个网格点来作为待裂变的网格点；其中M是均匀天线阵列中阵元个数；

步骤5、当所选出的待裂变的网格点中，存在至少一个待裂变的网格点的左右间隔同时大于预设的间隔阈值，即存在至少一个需要裂变的网格点时，则转入步骤6；否则，返回步骤2；

步骤6、对当前离格波达方向估计模型进行一次裂变，在裂变过程中，在每一个需要裂变的网格点的前一个网格点和后一个网格点的中点位置各插入一个新的网格点，以更新离格波达方向估计模型；

步骤7、返回步骤2。

2.根据权利要求1所述的一种基于网格部分细化的波达方向估计方法，其特征是，所述预设的迭代停止阈值τ的取值范围在10^-4～10^-3之间。

3.根据一种基于网格部分细化的波达方向估计方法，其特征是，所述预设的次数阈值的取值范围在500～1000之间。

4.根据权利要求1所述的一种基于网格部分细化的波达方向估计方法，其特征是，所述预设的间隔阈值的取值范围在1°～3°之间。