CN110109050A - 嵌套阵列下基于稀疏贝叶斯的未知互耦的doa估计方法 - Google Patents

Abstract

本发明公开了一种嵌套阵列下基于稀疏贝叶斯的未知互耦的DOA估计方法，属于自适应阵列信号处理领域的DOA技术。本发明基于阵列非均匀结构的特殊性，在未知阵元之间互耦矩阵的情况下，通过将物理阵列存在互耦的情况下的接收数据转换到虚拟阵列上，从而增加了信号的自由度；然后通过互耦矩阵与导向矩阵的关系，设置新的导向矩阵；接着使用基于稀疏贝叶斯的方法进行估计信号功率后验分布的均值和方差；最后通过求解信号功率的二范数的均值来估计信号的DOA。本发明降低了阵列的字典矩阵的复杂度，使得估计的DOA更加准确，并且还充分利用了嵌套阵列的高自由度的优势，使得估计信号的个数超出了物理阵列的孔径；且本发明不需要对信号的存在互耦矩阵进行消除。

Description

嵌套阵列下基于稀疏贝叶斯的未知互耦的DOA估计方法

技术领域

本发明涉及自适应阵列信号处理领域的波达方向角(DOA，direction-of-arrival)技术，具体是涉及在嵌套阵列条件下阵列接收信号存在未知互耦的DOA估计方法。

背景技术

阵列信号处理在雷达，声纳，通信领域中有很重要的应用。近几年来，在有限的阵列数的情况下，测向自由度高，测向精度高受到了广大研究人员的关注。然而，对于均匀阵列，其阵列自由度很低，测向信号的个数受限于物理阵列的个数。2010年，Piya Pal和P.P.Vaidyanathan等人提出的一种稀疏阵列结构，其方法是通过接收数据的相关矩阵向量化处理后得到的虚拟均匀阵列结构，从而极大的增加了阵列的自由度。随后，一系列的稀疏阵列被提出，包括嵌套阵列，互质阵列，k-次扩展互质阵列等。但是，原有的基于稀疏阵列均是在理想的假设条件下进行的。在实际应用中，存在各种实际误差，例如，阵元互耦效应等。传感器之间的互耦对参数估计有不利的影响，一些经典的算法，例如l1-svd，MUSIC，ESPRIT等算法不在适用。

在嵌套阵列存在互耦的情况下，考虑二级嵌套阵列结构，其中子阵一的阵元数为M₁,子阵二的阵元数为M₂，且有M₁+M₂＝M。假设有K个来自不同方向的窄带远场信号，其入射角度分别为θ_k,k＝1,2,...,K。在信号存在互耦的情况下n时刻阵列的接收数据矩阵为

x(n)＝CA(θ)s(n)+e(n),n＝1,2,…,N

其中，A(θ)＝[a(θ₁)a(θ₂)…a(θ_K)],s(n)＝[s₁(n),s₂(n),…,s_K(n)]^T，(·)^T表示转置，导向阵列矢量其中e表示自然底数，λ表示波长，d_i(i＝1,2,…,M)表示阵元间距，信号s(n)与噪声e(n)统计独立。

对于嵌套阵列，互耦矩阵随着阵元的间距不同而有所不同。通常情况下，互耦效应与相邻单元之间的单元间距成反比。假定这里的互耦矩阵的自由度为m＝4，即阵元间距大于1.5λ时阵元互耦可以忽略为0，根据嵌套阵列的特点，可以得到在嵌套阵列下的互耦矩阵为

故此，上述等式的相关矩阵为：

其中，ρ_k(k＝1,2,…,K)表示第k个信号的功率，为了表示方便，用A代替A(θ)，表示噪声功率，I_M表示M×M的单位矩阵，(·)^H表示共轭转置。

将上述嵌套阵列转化到虚拟均匀阵列上，即将上述相关矩阵向量化后，可以得到：

其中， e_m表示一个包含M个元素的列向量，其第m个元素的值等于1，其余的元素为零。(·)^*表示共轭，表示Kronecker积,表示Khatri-Rao积。

由于此时，阵列的接收数据存在未知互耦，未进行互耦校正的测向方法是无法估计出DOA值。目前国内外一些研究学者对于此问题展开了深度研究。Junpeng Shi等人提出了一种广义互质阵列，即通过一种具有两个柔性共素因子的GNA结构，去扩大两个级联均匀线性子阵的阵元间距，从而降低了阵列的互耦效应；Chun-Lin Liu等人提出了一种超嵌套阵列，其方法是通过将嵌套阵列第一级的稠密子阵进行重新排列，去降低第一级子阵的稠密性，从而降低了阵列的互耦效应，随后，Chun-Lin Liu又给出了高阶超嵌套阵列的具体表示形式，并指出二阶超嵌套阵列是其特殊的形式；陈璐等人通过对嵌套阵列结构进行优化，提出两种平移嵌套阵列结构，该结构使原嵌套阵列的一级阵的阵元稀疏度大大提高，降低了阵元间互耦效应。然而，这些方法都是通过布阵的形式进行降低阵列的互耦。为了在实际的稀疏阵列模型下进行估计未知互耦的DOA，ElieBouDaher等人提出了一种开路环的方法进行解耦，通过解耦合矩阵的逆来消除阵列之间的互耦效应，然而其方法并没有考虑非网格的问题；J.Dai等人提出了一种迭代的方法进行估计DOA，然而，此方法对待非均匀阵列作为均匀线性阵列的一个子集，并没有发挥出非均匀阵列的高自由度的优势。

发明内容

本发明的发明目的在于：针对现有的嵌套阵列DOA方法无法解决信号数超过阵列孔径且阵元间存在互耦的情况下，本发明公开了一种接收数据存在未知互耦时的DOA估计方法。

本发明的嵌套阵列下基于稀疏贝叶斯的未知互耦的DOA估计方法，包括下列步骤：

步骤1：基于嵌套阵列得到的接收数据x(n)，计算数据的相关矩阵R_X，并通过向量化操作，得到参量z＝vec(R_X)；

由于实际应用中不可避免的存在各种误差通过引入这种误差，可以得到虚拟阵列的接收数据即本发明中，设置虚拟阵列的接收数据为y。

其中，相关矩阵N表示快拍数，符号(·)^H表示共轭转置；

符号(·)^T表示转置；

误差服从渐进正态分布，且的均值为0，方差为单位矩阵；

步骤2：基于稀疏贝叶斯估计K个来自不同方向的信号的波达方向角DOA：

步骤201：将感兴趣的信号的角度区域划分为D份，得到D个网格角度：

用表示每个网格角度的网格失配量，设置网格失配矩阵Δ_D为：

其中，I_m×m表示m×m的单位矩阵，m表示互耦矩阵的自由度；

设置第一导向矢量为：

并设置第二导向矢量为：其中

参量其中i＝l₁,l₂,…,l_D；

且

其中，符号[·]_j,k表示矩阵第j行第k列的元素，[·]_j表示向量第j个元素，a(θ_i)表示入射角为θ_i的阵列导向矢量；

步骤202：初始化迭代参量，包括噪声矢量α₀，向量参数α，网格失配矩阵Δ_D；以及设置迭代收敛阈值τ和迭代次数上限；

其中，α₀为误差的渐进正态分布的精度，其概率密度函数为：p(α₀；c,d)＝Γ₁(α₀|c,d)，其中参数c和d为接近于零的固定值；

向量参数α中的每个元素α_i服从分布：α_i～Γ₂(1,ρ),i＝1,2,…,m²D，其中参数ρ为预设值，其取值大于0；

步骤203：求解信号的后验分布的均值μ和方差Σ：

其中，

且表示噪声功率，e_m表示一个包含M个元素的列向量，其第m个元素的值等于1，其余的元素为零；其中M表示嵌套阵列的阵元数；

步骤204：更新向量参数α和噪声矢量α₀；

其中，更新后的向量参数α的每个元素为：

再基于更新后的向量参数α更新噪声矢量α₀，得到更新后的噪声矢量

即，首先计算每个元素的更新值保存向量参数α更新前的值并更新向量参数α的每个元素从而得到更新后的向量参数α，再基于当前向量参数α计算并令

步骤205：更新网格失配矩阵量Δ_D；

估计网格失配矩阵Δ_D中D个网格角度的更新值其中n＝l₁,l₂,…,l_D；

并根据得到更新后的网格失配矩阵量Δ_D；

更新值

参量

参量

其中，F(n,：)表示第n行元素，()_-n表示去除第n个位置的元素，符号(·)^*表示共轭，diag()表示取对角元素，Re{·}表示取实部；

步骤206：判断是否满足迭代收敛条件，若是，则基于当前得到的均值μ和方差Σ，以及网格失配矩阵Δ_D执行步骤207；否则，返回步骤203；

所述迭代收敛条件为：||α^t+1-α^t||²/||α^t||²＜τ或迭代次数达到预设的迭代次数上限；

其中，α^t+1、α^t分别表示更新前、后的向量参数α；

步骤207：计算经过去耦合预处理后的信号的功率的数学期望其中Tr表示矩阵的迹；

通过使用分块的思想进行预处理，即对m²×D个构成的向量进行分块处理，得到D个大小相同的块；

再基于D个块进行谱峰搜索，将搜索到的前K个最大的峰值对应的网格角度分别加上各网格角度对应的网格失配量，作为K个来自不同方向的信号的波达方向角的估计结果。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

在嵌套阵列阵元之间存在互耦的情况下,既能解决DOA估计的问题，也考虑了信号的离格、快拍数有限所导致的误差等问题，从而降低了阵列的字典矩阵的复杂度，使得估计的DOA更加准确。是一种新的DOA估计方法。本发明充分利用了嵌套阵列的高自由度的优势，使得估计信号的个数超出了物理阵列的孔径，并且，本发明不需要对互耦矩阵进行预估计，更不需要对阵元位置进行重新排列。

附图说明

图1为贝叶斯网络示意图；

图2本发明与其他算法的对比图；

图3本发明与未经预白化处理的LASSO算法的对比图；

图4 RMSE随SNR的性能对比图；

图5 RMSE随snapshot的性能对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

本发明针对现有的嵌套阵列DOA方法无法解决信号数超过阵列孔径且阵元之间存在互耦的情况下，公开了一种信号基于稀疏贝叶斯的未知互耦的DOA估计方法，该方法可测信号数超过了实际物理阵列个数。

本发明基于嵌套阵列结构的特殊性，在阵元之间存在互耦的情况下，通过三部分处理实现DOA估计：(1)重构出向量化后的接收数据导向矩阵；(2)在考虑了各种实际误差的情况下，重构出新的接收信号向量的形式；(3)利用基于稀疏贝叶斯的方法来求得最终的DOA值。本发明不用提前预估互耦矩阵的互耦系数，而是利用互耦矩阵与导向矢量的关系，重新构造出新的导向矩阵，从而进行对信号的波达方向角进行估计。在信噪比(SNR)较低时，也依然能够对信号进行准确的估计。其各部分具体为：

(1)重构向量化后的接收数据的导向矩阵。

根据带状对称的toeplitz矩阵的特性，经过推算，可以得到

其中，m为互耦矩阵的自由度。对于嵌套阵列，参量T(θ_i)可通过下面的方式得到：

其中，

上述公式中的符号[·]_j,k表示矩阵第j行第k列的元素，[·]_j表示向量第j个元素，a(θ_i)表示入射角为θ_i的阵列导向矢量。

从而，可以得到

故此，上述嵌套阵列转化到虚拟均匀阵列的接收数据可以转化为：

其中，

DOA估计的问题能够转化为下面的优化问题：

其中，ε是一个很小的预设定的参数，是一个维度为m²D×1的向量，该DOA估计问题可以通过LASSO(Least Absolute Shrinkage and Selection 0perator)块匹配的算法进行求解。然而，LASSO算法难以解决字典矩阵的非网格问题。故本发明通过采用一种稀疏贝叶斯压缩感知的方法进行估计DOA。

(2)考虑各种实际误差的情况下新的接收信号向量的形式。

由于实际生活中，接收信号通过有限的快拍数进行测量的。测试的结果数据会存在误差，并且这个误差值满足的是渐进正态分布。故此，得到虚拟阵列的接收数据为

其中，表示在有限的快拍数下得到的信号的协方差矩阵，Δz表示偏差量或误差，其服从渐进正态分布，可以简写为Δz～AsN(0,W)，其中0表示Δz的均值，W表示Δz的方差。

由于实际应用中是通过有限快拍数来得到阵列的协方差矩阵的，故此，使用代替上述的W，一种经过预白化处理后得到的虚拟阵列的单快拍接收数据为

其中，

假设整个DOA空间均匀的划分为D份，标记为当信号的角度不在划分的网格上时，它不可避免的产生网格失配的情况，从而会带来一定的误差。如果真实的角度是θ_k，它能够表示为下面的表达式

其中，表示最邻近θ_k的网格角度，表示每个网格角度对应的网格失配量，l_k∈{l₁,l₂,…,l_D}，通常假设其在区间上服从均匀分布，其中，即l表示网格角度间隔。通过使用一阶泰勒展开式，并且忽略高于一阶的泰勒展开式，则有接收信号的新的导向矢量为

其中，k∈{1,2,…,K},表示的一阶导数，记为则网格上所有的新的导向矢量可以表示为：

其中，

其中，I_m×m表示m×m的单位矩阵，即Δ_D表示网格失配量(网格失配矩阵)，也可称为网格误差或者网格偏移量。

因此，可以得到整个网络的失配量的导向矩阵的表达式为

从而可以得到虚拟阵列的接收数据可重新表示为

由于是需要进行估计的，故此，整个网络的稀疏贝叶斯的字典矩阵可以表示为

其中， D表示在l_d处的网格失配量。

在压缩感知理论中，上述的信号模型可以表示为

其中，虚拟阵列的接收数据y也可以称为观测向量样本，可以表示为是一个维度为m²×1的向量。即DOA估计就是基于块匹配的思想。当是需要估计的角度θ_k(k∈{1,2,…,K})的时候，向量中的第l_k块的维度为m²×1的值近似等于否则，此块的值为零向量。通过前述的推导，可以得到即误差服从渐进正态分布。则可以得到其概率密度函数为

其中，α₀表示误差满足的渐进正态分布的精度，并且α₀服从Gamma分布，其概率密度函数为

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

其中，c,d是一个固定的参数，通常要求c,d接近于零。从而可以得到y的条件概率密度函数为

其中，

为了将一个压缩感知的稀疏信号问题变为一个求最大后验分布的贝叶斯学习问题，需要通过观测向量样本y去求解未知稀疏信号矢量和α₀的最大后验分布。

要想求得未知互耦条件下的稀疏信号矢量参数与未知误差的似然函数，还需要知道未知参数的先验信息。而在基于贝叶斯学习的压缩感知理论中，比较常用的参数的先验信息为拉普拉斯密度分布。因此，当给定了观测向量样本y，并且在已知过完备字典矩阵的情况下，假设观测值样本y关于未知参数与α₀的似然函数满足渐进正态分布模型，而稀疏信号矢量参数的先验信息满足拉普拉斯密度分布，便可以通过贝叶斯学习方法求得未知稀疏信号矢量参数的最大后验分布。不过由于服从拉普拉斯密度分布的未知参数先验信息的求解过程十分复杂。因此，本发明首先假设未知的稀疏信号矢量参数中的每一个元素均服从渐进正态分布，并且该渐进正态分布的均值为0，方差为即

在未知的稀疏信号矢量参数中，其各个参数之间均是相互独立的。所以，未知稀疏信号矢量参数的先验分布可以表示为

其中，假设未知向量参数α中的每一个元素[α]_i服从一个Γ分布，即

α_i～Γ(1,ρ),i＝1,2,…,m²D

其中，ρ是一个未知的参数，通常要求ρ＞0。同样，在未知向量参数α中，每个元素之间是相互独立的，所以，可以得到未知向量参数α的先验分布为

故此，可以得到未知稀疏信号矢量参数关于参数ρ的先验分布为

由于未知稀疏信号矢量服从于均值为0的渐进正态先验分布，所以可以保证该参数中的大部分元素取值为0。这样，不仅可以保证了压缩感知理论中要求原始信号是稀疏的这一前提条件，同时也降低了接下来求解未知稀疏信号矢量参数的最大后验分布过程的计算量。其可以表示为

因为观测样本y与未知参数向量α无关，而未知稀疏信号矢量与误差的方差参数α₀无关，故上式可写为

对于上式，通过前面的推导，知道上述公式是已知的，也就是说，根据观测向量样本y的高斯似然函数与未知稀疏信号矢量参数的先验分布可以得到未知稀疏信号矢量参数的后验分布。

通过上述的推导，可以得出这样一个结论：利用稀疏贝叶斯DOA估计算法是可行的。故此，可以通过后验分布来求解相应的未知参数。为了表述方便，写出上述的联合概率分布为

其贝叶斯网络如图1所示。

(3)基于稀疏贝叶斯DOA估计。

基于上述c部分推导的后验分布由于已知的为

其中，Λ＝diag(α)。如何利用上述已知的先验信息去推导想要的DOA估计结果的具体处理过程为：

首先，先进行推导上述后验分布的均值和方差，即

从而，可以得到

其中，对于最后一个式子与倒数第二个式子进行比较得到

即

即

也就是说，均值为方差为写成概率密度函数即为

通过矩阵求逆引理(即Woodbury matrix identity)的应用，有

因此，可以得到，上述后验分布的均值和方差为

接下来，再利用上述的均值和方差来获取未知参数α₀，α和网格失配量Δ_D的表达式，以便于对上述的后验分布进行更新迭代，从而得到一个更加准确的DOA信息。利用最大似然估计对上式求解，即可以得到未知参数α₀和α的估计值，为

由于故此，上式可以转化为

两边同时取ln函数，有

对其关于α₀求导，让其等于零，可以得到

进而，可以得到

同理，对其关于α求导，让其等于零，可以得到

上式可以转化为

从而可以得到

虽然现在已经得到了未知参数α₀和α的表达式，但是，无法将其与前式中的均值与方差进行迭代。为了进行更新上述的后验分布的均值与方差，本发明使用期望值最大化算法(Expectation Maximization,EM)来求解最大似然函数，即，将上述推导的关于α₀和α的表达式的两边求数学期望，可以得到

而对于网格失配量Δ_D，它的估计是通过最大化联合概率密度来得到的，即

而不包含网格失配量Δ_D，故此，只需最大化即可。即

而E{lnp(Δ_D)}为常数矩阵，有负号的存在，上式实际上可以转化为最小化，即

其中，

δ＝diag(Δ_D)

对上式关于δ求导，则有2Fδ-2ν，令其为零，即Fδ＝ν，由于F是一个奇异矩阵，因此，根据矩阵展开式，有

其中，F(n,：)表示F的第n行。()_-n表示去除第n的位置的元素。由于将每个预估的角度的网格失配量规定在内，如果超出，则通过下式进行更新

至此，通过未知稀疏信号矢量的后验信息均值μ和方差Σ的函数去更新 然后再将其带入未知稀疏信号矢量的后验信息均值μ和方差Σ的函数，而均值μ和方差Σ又是的函数，故这两种过程进行反复的迭代，直至收敛。收敛的判决条件为||αⁱ⁺¹-αⁱ||²/||αⁱ||²＜τ，其中τ为预设阈值，其一般为1e-3。

通过上述的方法，其DOA估计的方法是没有直接去估计信号的功率，而是估计信号的功率的均值与方差。故此，上述方法需要将其转化到去估计与信号的功率相关的问题上来才可以进行DOA估计。使用求信号的功率的二范数的数学期望来进行估计信号的DOA，即有

其中，表示去耦合处理后的信号的功率的数学期望。

由于信号存在互耦的情况，在经过去耦合预处理后，每一个信号入射角所对应的方向的维度从变为了这样不利于进行谱峰搜索。因此，本发明需要通过分块的思想去估计信号的DOA。即如果有一个方向上存在信号的入射，那么此块便存在非零值，否则，此块的值为零。通过上述方法进行估计后，还需要进行预处理才能得到最终的入射信号DOA估计的结果。也就是说，为了得到最终的DOA估计结果，需要对DOA估计分块进行预处理，即利用l₂范数来得到每一块的值，从而将每一块的维度从降为即将每一块从一个向量变为一个具体的值，从而便于进行谱峰搜索。同理，上述算法进行估计的网格失配量Δ_D也需要通过降维的思路进行分块处理。通过分块预处理后所得到的每一个角度的和网格失配量，便可以得到最终的DOA估计的信息。

即本发明的嵌套阵列条件下阵列接收信号存在互耦的DOA估计方法，具体包括下列步骤：

步骤1：通过嵌套阵列得到接收数据x(n)，进而得到数据的相关矩阵R_X；通过向量化操作，得到z＝vec(R_X)，当考虑快拍数有限所导致的误差时，即本发明中，考虑到实际工程应用，将得到虚拟阵列的接收数据表示为

考虑到实际工程应用，本发明中相关矩阵求解方法为N表示快拍数，即估计的次数，符号(·)^H表示共轭转置。

表示快拍数有限所导致的阵列误差，其满足的概率密度函数为并且，

步骤2：采用基于稀疏贝叶斯的在线学习的思想，对基于稀疏贝叶斯估计算法的先验信息α₀，α，Δ_D进行初始化，设定ρ,c,d,τ的初始值，迭代的总次数iters以及第一次迭代的次数iter＝1。

本具体实施方式中，取ρ＝0.01；c＝d＝10^-4；τ＝10^-3。

完成上述初始化后，再将感兴趣的信号的角度区域划分为D份，即即目标假设在划分的角度区域内，从而，得到导向矩阵为

当考虑快拍数有限所导致的数据误差时，有

由于在实际应用中，估计的角度或许不在所划分的网格上，故其存在网格失配(网格误差)的情况。对于网格失配的情况下，接收信号的新的导向矢量为

其中，

步骤3：求解后验分布的均值μ和方差Σ：

步骤4：求解α和噪声矢量α₀：

为了公式的简洁表述，定义表示噪声功率，

步骤5：更新网格失配量Δ_D，即有

其中，F(n,：)表示第n行元素，()_-n表示去除第n个位置的元素，符号(·)*表示共轭，diag()表示取对角元素，Re{·}表示取实部，则有

步骤6：判断是否继续更新迭代，即判断||αⁱ⁺¹-αⁱ||²/||αⁱ||²＜τ，如果满足条件，则返回步骤3，继续执行；否则，结束，得到网格失配量Δ_D。

步骤7：通过将上述的估计值做进一步的预处理，便可以得到真实的信号的方向角度。

首先计算经过去耦合处理后的信号的功率的数学期望其中j＝1,2,…,m²·D，Tr表示矩阵的迹；

通过使用分块的思想进行预处理后，再基于D个块进行谱峰搜索，即从D个块中查找前K个最大的峰值，再基于前K个最大的峰值所对应的网格角度加上各网格角度对应的网格失配量得到信号的方向角度。

为使本发明的目的、技术方案和技术效果更加清楚，通过仿真实验对本发明作进一步地详细描述。

仿真实验条件一：采用二级嵌套阵列排列，阵元个数为6，子阵1的阵元数为M₁＝3，子阵2的阵元数为M₂＝3，阵元位置为[1，2，3，4，8，12]，互耦系数为：c1＝0.2121+0.2121i；c2＝-0.0882+0.1214i；c3＝-0.0588+0.0809i；互耦矩阵的自由度为m＝4；信号均为非相干的。网格搜索间距为1°；利用稀疏贝叶斯的方法进行的估计的参数为ρ＝10^-2，网格失配量为D为网格的个数，稀疏贝叶斯迭代次数为2000次，初始噪声功率为利用LASSO的压缩感知方法进行估计的参数为λ_t＝1.28，信号的角度为[-21.1,

-5.2,10.1,30.1],快拍数1000，信噪比SNR＝10dB，仿真结果如图2所示。

通过图2，可以明显的发现，传统的物理阵列、预白化处理后的LASSO的压缩感知方法均无法估计出信号的DOA，然而，本发明可以很好的估计出信号的角度，从而也验证了本发明的可行性。

仿真实验条件二：考虑信号的角度为[-45.3,-30.51,-14.98,0.5,15.5,30.1]，网格搜索间距为0.5°，其余的仿真条件与仿真实验条件一相同，本发明与未进行预白化处理后的LASSO的压缩感知方法进行了对比，仿真结果如图3所示。

通过图3可以明显的发现，虽然未进行预白化处理后的LASSO的压缩感知方法可以估计出信号的角度，然而，会出现很多的虚假的伪峰，此伪峰与真实的信号角度的功率相当，无法进行去除伪峰操作，其不利于进行估计信号，无法进行分辨出哪些是想要的信号的方向。此外，其也无法估计信号的角度不在网格上的DOA，估计的精度不准确。但是，本发明可以很好的估计出信号的角度，并且，通过离格估计，可以很好的降低网格的稀疏度，其估计信号的方向明显优于未进行预白化处理后的LASSO算法。

仿真实验条件三：研究RMSE随快拍数、信噪比(SNR)变化的情况下的性能对比情况。当仿真RMSE随快拍数变化时，SNR选择为SNR＝10dB；当考虑RMSE随SNR变化时，快拍数选择为1000。信号角度为-5.24°，16.15°,其余的条件同仿真实验条件一。本发明与未进行预白化处理后的LASSO算法、OMP算法进行性能对比。为消除实验的随机性对实验的影响，本次实验的蒙特卡罗次数为500次。从图4，5中可以明显的发现，未进行预白化处理后的LASSO算法无法估计出信号离格的问题，误差较大；OMP算法由于受到最小二乘原理及正交完备基的制约，角度的设置受限，当两个角度间隔很近的时候，从而导致估计精度很差，然而，本发明却可以很好的进行估计信号的角度，在快拍数及SNR很低的情况下，性能也很好。

本发明公开了在嵌套阵列下基于稀疏贝叶斯的未知互耦的DOA估计方法，基于阵列非均匀结构的特殊性，在未知阵元之间互耦矩阵的情况下，通过将物理阵列存在互耦的情况下的接收与导向矩阵的关系，设置新的导向矩阵；接着使用基于稀疏贝叶斯的方法进行估计信号功率后验分布的均值和方差；最后通过求解信号功率的二范数的均值来估计信号的DOA。本发明不仅可以估计出信号的DOA，而且还考虑了信号的离格问题，从而降低了阵列的字典矩阵的复杂度，使得估计的DOA更加准确，并且还充分利用了嵌套阵列的高自由度的优势，使得估计信号的个数超出了物理阵列的孔径，并且，本发明不需要对信号的存在互耦矩阵进行消除。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (6)

1.嵌套阵列下基于稀疏贝叶斯的未知互耦的DOA估计方法，其特征在于，包括下列步骤：

步骤1：基于嵌套阵列得到的接收数据x(n)，计算数据的相关矩阵R_X，并通过向量化操作，得到参量z＝vec(R_X)；

并设置虚拟阵列的接收数据为y，且

其中，相关矩阵N表示快拍数，符号(·)^H表示共轭转置；

符号(·)^T表示转置；

误差服从渐进正态分布，且的均值为0，方差为单位矩阵；

步骤2：基于稀疏贝叶斯估计K个来自不同方向的信号的波达方向角DOA：

步骤201：将感兴趣的信号的角度区域划分为D份，得到D个网格角度：

用表示每个网格角度的网格失配量，设置网格失配矩阵Δ_D为：

其中，I_m×m表示m×m的单位矩阵，m表示互耦矩阵的自由度；

设置第一导向矢量为：

并设置第二导向矢量为：其中

参量其中i＝l₁,l₂,…,l_D；

且

其中，符号[·]_j,k表示矩阵第j行第k列的元素，[·]_j表示向量第j个元素，a(θ_i)表示入射角为θ_i的阵列导向矢量；

步骤202：初始化迭代参量，包括噪声矢量α₀，向量参数α，网格失配矩阵Δ_D；以及设置迭代收敛阈值τ和迭代次数上限；

其中，α₀为误差的渐进正态分布的精度，其概率密度函数为：p(α₀；c,d)＝Γ₁(α₀|c,d)，其中参数c和d为接近于零的固定值；

向量参数α中的每个元素α_i服从分布：α_i～Γ₂(1,ρ),i＝1,2,…,m²D，其中参数ρ为预设值，其取值大于0；

步骤203：求解信号的后验分布的均值μ和方差Σ：

其中，

且 表示噪声功率，e_m表示一个包含M个元素的列向量，其第m个元素的值等于1，其余的元素为零；其中M表示嵌套阵列的阵元数；

步骤204：更新向量参数α和噪声矢量α₀；

其中，更新后的向量参数α的每个元素为：

再基于更新后的向量参数α更新噪声矢量α₀，得到更新后的噪声矢量

即，首先计算每个元素的更新值保存向量参数α更新前的值并更新向量参数α的每个元素从而得到更新后的向量参数α，再基于当前向量参数α计算并令

步骤205：更新网格失配矩阵量Δ_D；

估计网格失配矩阵Δ_D中D个网格角度的更新值其中n＝l₁,l₂,…,l_D；

并根据得到更新后的网格失配矩阵量Δ_D；

更新值

参量

参量

其中，F(n,：)表示第n行元素，()_-n表示去除第n个位置的元素，符号(·)^*表示共轭，diag()表示取对角元素，Re{·}表示取实部；

步骤206：判断是否满足迭代收敛条件，若是，则基于当前得到的均值μ和方差Σ，以及网格失配矩阵Δ_D执行步骤207；否则，返回步骤203；

所述迭代收敛条件为：||α^t+1-α^t||²/||α^t||²＜τ或迭代次数达到预设的迭代次数上限；

其中，α^t+1、α^t分别表示更新前、后的向量参数α；

步骤207：计算经过去耦合预处理后的信号的功率的数学期望其中Tr表示矩阵的迹；

由m²·D个构成向量并将向量均分为D个块；

再基于D个块进行谱峰搜索，将搜索到的前K个最大的峰值对应的网格角度分别加上各网格角度对应的网格失配量，作为K个来自不同方向的信号的波达方向角的估计结果。

2.如权利要求1所述的方法，其特征在于，步骤205中，更新值的计算方式可以替换为：

其中，l表示相邻网格角度的角度间隔值。

3.如权利要求1或2所述的方法，其特征在于，步骤202中，噪声矢量α₀的初始化方式为：初始化且参数c＝d＝10^-4。

4.如权利要求1或2所述的方法，其特征在于，步骤202中，向量参数α的初始化方式为：初始化且参数ρ＝0.01。

5.如权利要求1或2所述的方法，其特征在于，步骤202中，网格失配矩阵Δ_D的初始值设置为：

6.如权利要求1或2所述的方法，其特征在于，步骤202中，迭代收敛阈值τ＝10^-3。