基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法与流程

本发明属于信号处理技术领域，尤其涉及对雷达信号、声学信号及电磁信号的波达方向估计，具体是基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法，可用于无源定位和目标探测。

背景技术：

波达方向(direction-of-arrival,doa)估计是阵列信号处理领域的一个重要分支，它是指利用阵列天线接收空域信号，并通过现代信号处理技术和各类优化方法实现对接收信号统计量进行有效处理，以实现信号的doa估计，在雷达、声呐、语音、无线通信等领域有着重要的应用价值。

doa估计方法的自由度是指其能够同时分辨的入射信号源的个数，作为实际系统应用中的一个重要衡量指标，决定着系统的总体复杂度。现有的doa估计方法通常采用均匀线性阵列进行信号的接收与建模，但是基于均匀线性阵列方法的自由度是受实际天线阵元个数限制的。具体而言，对于一个包含l个天线阵元的均匀线性阵列，其自由度为l-1，即最多只能分辨l-1个入射信号。因此，当某个空域范围内入射信号源的个数大于或等于阵列中天线阵元的个数时，现有采用均匀线性阵列的方法将无法进行有效的doa估计。

为了增加自由度，传统方法需要通过增加物理天线阵元及相应的射频模块来实现，这造成了系统计算复杂度和硬件复杂度的增加。因此，现有采用均匀阵列的doa估计方法在自由度性能与计算复杂度之间存在着一定的利弊权衡问题。如何在物理天线阵元个数一定的情况下提升波达方向估计方法的自由度，对于提升方法在实际系统应用中的经济性和实用性有着重要的意义。

技术实现要素：

本发明的目的在于针对上述现有技术存在的不足，提出一种基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法，通过利用互质阵列的特性推导虚拟域等价接收信号，并通过虚拟阵列协方差矩阵重建的方法重建空间功率谱，以提升估计方法的自由度，实现在信号源个数大于或等于物理天线阵元个数情况下的有效doa估计。

本发明的目的是通过以下技术方案来实现的：基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法，该方法包含以下步骤：

(1)在接收端使用q个物理天线阵元架构互质阵列，并通过互质阵列接收入射信号；

(2)假设有k个来自θ1,θ2,…,θk方向的远场窄带非相干信号源，则q×1维互质阵列接收信号y(t)可建模为：

其中，sk(t)为信号波形，n(t)为噪声分量，与各信号源相互独立，a(θk)为θk方向的导引向量，表示为

其中，uq,q＝1,2,…,q表示互质阵列中第q个物理天线阵元的实际位置，且u1＝0，λ表示信号波长，[·]^t表示转置操作；共采集t个采样快拍，得到采样协方差矩阵

这里(·)^h表示共轭转置；

(3)计算互质阵列接收信号所对应的等价虚拟信号：向量化互质阵列接收信号的采样协方差矩阵获得虚拟阵列等价接收信号z：

其中，为q²×k维虚拟阵列导引矩阵，p＝[p1,p2,…,pk]^t包含k个入射信号源的功率，为噪声功率，i＝vec(iq)。这里，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)^*表示共轭操作，表示克罗内克积，iq表示q×q维单位矩阵。向量z对应的虚拟阵列中各虚拟阵元的位置为s：

s(i,j)＝{ui-uj|i,j＝1,2,…,q}。

去除集合s中各位置上重复的虚拟阵元，得到一个非均匀的虚拟阵列sn，其对应的等价虚拟信号可通过选取向量z中相对应的元素获得；

(4)构造虚拟阵列协方差矩阵：首先选取非均匀虚拟阵列sn中以0为中心连续均匀排列的一段虚拟阵元，形成一个包含l个虚拟阵元的均匀虚拟阵列su，其相应的虚拟阵元位置为-lvd到lvd之间的连续位置，其中，d为单位间隔，取入射窄带信号波长的一半，即d＝λ/2，

相应地，该均匀虚拟阵列的等价信号可通过截取中与该l个虚拟阵元所对应的位置上的元素获得，维度为l×1。接着，虚拟阵列协方差矩阵rv可根据构建如下toeplitz结构的(lv+1)×(lv+1)维矩阵得到：

其中，表示位置为id的虚拟阵元所对应的等价接收信号；

(5)构建虚拟阵列协方差矩阵稀疏重建优化问题并求解：首先，将波达方向角的角度域范围等间距地划分为个网格点即然后，根据虚拟阵列协方差矩阵rv构建如下以向量和噪声功率为变量的优化问题：

其中，为维虚拟阵列导引矩阵，其对应于虚拟阵元位置为0到lvd之间连续的一段虚拟均匀阵列；向量包含个潜在来波方向上的信号功率，其对应的对角矩阵为∈为阈值常数，用于约束协方差矩阵的重建误差；保证了各方向上的信号功率值大于等于零；‖·‖0和‖·‖f分别表示0范数和f范数，为(lv+1)×(lv+1)维单位向量。将上述非凸优化问题转化为凸优化问题，并求得最优解

(6)通过谱峰搜索获得波达方向估计结果：以x轴为个均匀分布的空间网格点来波方向，y轴为优化值中所包含的元素，画出空间功率谱。寻找空间功率谱上的峰值，并将这些峰值所对应的响应值从大到小排列，取前k个峰值所对应的x轴角度方向，即为波达方向估计结果。

进一步地，步骤(1)所述的互质阵列结构可描述为：首先选取一对互质整数m、n；然后，构造一对稀疏均匀线性子阵列，其中第一个子阵列包含m个间距为nd的天线阵元，其位置为0,nd,…,(m-1)nd，第二个子阵列包含n个间距为md的天线阵元，其位置为0,md,…,(n-1)md；接着，将两个子阵列按照首个阵元重叠的方式进行子阵列组合，获得实际包含q＝m+n-1个天线阵元的非均匀互质阵列架构。

进一步地，步骤(1)所述的互质阵列结构可描述为：首先选取一对互质整数m、n，且m<n；然后，构造一对稀疏均匀线性子阵列，其中第一个子阵列包含2m个间距为nd的天线阵元，其位置为0,nd,…,(2m-1)nd，第二个子阵列包含n个间距为md的天线阵元，其位置为0,md,…,(n-1)md；接着，将两个子阵列按照首个阵元重叠的方式进行子阵列组合，获得实际包含q＝2m+n-1个天线阵元的非均匀互质阵列架构。

进一步地，步骤(4)所述的虚拟阵列协方差矩阵rv可通过以下方式等价获得：

进一步地，步骤(4)所述的虚拟阵列协方差矩阵rv可通过以下方式等价获得：首先将向量分割为lv+1个维度为(lv+1)×1的子向量，每个子向量包含向量中的第i个至第i+lv个元素，即：

则

由此，rv可通过求上式中四阶统计量的矩阵开平方获得。

进一步地，步骤(5)中所构建的非凸优化问题可通过凸松弛技术，将优化问题中的0范数替换为1范数，得到以下以向量和噪声功率为变量的凸优化问题：

其中‖·‖1表示1范数。

进一步地，步骤(5)中所构建的非凸优化问题可转化为如下以向量和噪声功率为变量的基追踪去噪优化问题：

其中ξ为正则化参数，用于权衡虚拟阵列协方差矩阵重建误差和向量的稀疏性。

本发明与现有技术相比具有以下优点：

(1)本发明充分利用了互质阵列的互质特性，通过虚拟阵列等价接收信号的推导，实现了在虚拟阵列协方差矩阵的构建。由于虚拟阵列所包含的虚拟阵元个数大于实际物理天线阵元的个数，基于虚拟域的信号处理为自由度性能的提升奠定了基础；

(2)本发明考虑了信号在空间域范围内呈现稀疏性这一特点，在虚拟域上构建协方差矩阵稀疏重建优化问题，以实现在信号源个数大于或等于物理天线阵元个数情况下的有效doa估计；

(3)与现有采用均匀阵列的方法相比，本发明所提方法在自由度、阵列孔径等方面的优势能够有效减少实际系统中物理天线阵元及射频模块单元的个数，在实际应用中体现出较为理想的经济性和实用性。

附图说明

图1是本发明的方法总体流程框图。

图2是本发明中组成第一类互质阵列的一对稀疏均匀子阵列结构示意图。

图3是本发明中第一类互质阵列的结构示意图。

图4是本发明中组成第二类互质阵列的一对稀疏均匀子阵列结构示意图。

图5是本发明中第二类互质阵列的结构示意图。

图6是本发明中采用第一类互质阵列重建的空间功率谱示意图。

图7是本发明中采用第二类互质阵列重建的空间功率谱示意图。

具体实施方式

以下参照附图，对本发明的技术方案和效果作进一步的详细说明。

对于doa估计方法在实际系统中的应用，我们往往希望能够采用较少的天线设备估计更多的入射信号源，但是受限于天线阵列结构、阵元个数等因素，现有方法在这两个方面无法同时实现优化，往往存在着一定的利弊权衡关系。为了提升doa估计方法的自由度，本发明提供了基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法，参照图1，本发明的实现步骤如下：

步骤一：在接收端使用q个物理天线阵元架构互质阵列。互质阵列结构主要包含以下两类，均适用于本发明所提doa估计方法。

第一类互质阵列结构如下：首先选取一对互质整数m、n；然后，参照图2，构造一对稀疏均匀线性子阵列，其中第一个子阵列包含m个间距为nd的天线阵元，其位置为0,nd,…,(m-1)nd，第二个子阵列包含n个间距为md的天线阵元，其位置为0,md,…,(n-1)md；单位间隔d取为入射窄带信号波长的一半，即d＝λ/2；接着，将两个子阵列按照首个阵元重叠的方式进行子阵列组合，参照图3，获得实际包含m+n-1个天线阵元的非均匀互质阵列架构。此时，物理天线阵元个数q＝m+n-1。

第二类互质阵列结构如下：首先选取一对互质整数m、n，且m<n；然后，参照图4，构造一对稀疏均匀线性子阵列，其中第一个子阵列包含2m个间距为nd的天线阵元，其位置为0,nd,…,(2m-1)nd，第二个子阵列包含n个间距为md的天线阵元，其位置为0,md,…,(n-1)md；接着，将两个子阵列按照首个阵元重叠的方式进行子阵列组合，参照图5，获得实际包含2m+n-1个天线阵元的非均匀互质阵列架构。此时，物理天线阵元个数q＝2m+n-1。

步骤二：采用互质阵列接收信号并建模。假设有k个来自θ1,θ2,…,θk方向的远场窄带非相干信号源，采用步骤一架构的非均匀互质阵列接收入射信号，得到q×1维互质阵列接收信号y(t)，可建模为：

其中，sk(t)为信号波形，n(t)为与各信号源相互独立的噪声分量，a(θk)为θk方向的导引向量，表示为

其中，uq,q＝1,2,…,q表示互质阵列中第q个物理天线阵元的实际位置，且u1＝0，[·]^t表示转置操作。共采集t个采样快拍，得到采样协方差矩阵

这里(·)^h表示共轭转置。

步骤三：计算互质阵列接收信号所对应的等价虚拟信号。向量化互质阵列接收信号的采样协方差矩阵获得虚拟阵列等价接收信号z：

其中，为q²×k维虚拟阵列导引矩阵，p＝[p1,p2,…,pk]^t包含k个入射信号源的功率，为噪声功率，i＝vec(iq)。这里，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)^*表示共轭操作，表示克罗内克积，iq表示q×q维单位矩阵。向量z对应的虚拟阵列中各虚拟阵元的位置为s：

s(i,j)＝{ui-uj|i,j＝1,2,…,q}。

去除集合s中各位置上重复的虚拟阵元，得到一个非均匀的虚拟阵列sn，其对应的等价虚拟信号可通过选取向量z中相对应的元素获得。

步骤四：构造虚拟阵列协方差矩阵。首先，选取非均匀虚拟阵列sn中以0为中心连续均匀排列的一段虚拟阵元，形成一个包含l个虚拟阵元的均匀虚拟阵列su(由于su中的虚拟阵元以零位对称分布，l始终为奇数)，其相应的虚拟阵元位置为-lvd到lvd之间的连续位置，其中

相应地，该均匀虚拟阵列的等价信号可通过截取中与该l个虚拟阵元所对应的位置上的元素获得，维度为l×1。接着，虚拟阵列协方差矩阵rv可根据构建如下toeplitz结构的(lv+1)×(lv+1)维矩阵得到：

其中，表示位置为id的虚拟阵元所对应的等价接收信号。由于均匀虚拟阵列su中的虚拟阵元关于原点对称，其对称阵元所对应的二阶等价接收信号统计量互为共轭关系，因此rv也可以等价表示为：

此外，rv亦可通过空间平滑技术获得，具体为：将向量分割为lv+1个维度为(lv+1)×1的子向量，每个子向量包含向量中的第i个至第i+lv个元素，即：

则

由此，rv可通过求上式中四阶统计量的矩阵开平方获得。

步骤五：设计虚拟阵列协方差矩阵稀疏重建优化问题并求解。首先，根据信号在空间域范围内的稀疏分布特性，将波达方向角的角度域范围等间距地划分为个网格点即然后，根据步骤四计算得出的虚拟阵列协方差矩阵rv构建如下以向量和噪声功率为变量的优化问题：

其中，为维虚拟阵列导引矩阵，其对应于虚拟阵元位置为0到lvd之间连续的一段虚拟均匀阵列；向量包含个潜在来波方向上的信号功率，其对应的对角矩阵为∈为阈值常数，用于约束协方差矩阵的重建误差；保证了各方向上的信号功率值大于等于零；‖·‖0和‖·‖f分别表示0范数和f范数，为(lv+1)×(lv+1)维单位向量。由于上述优化问题包含0范数这一非凸项，这将导致求解困难；为了获得优化解，我们考虑引入凸松弛技术，将上述优化问题中的0范数替换为1范数，得到以下以向量和噪声功率为变量的凸优化问题：

其中‖·‖1表示1范数。上述凸优化问题可等价写为以下以向量和噪声功率为变量的基追踪去噪优化问题：

其中ξ为正则化参数，用于权衡虚拟阵列协方差矩阵重建误差和向量的稀疏性。求解上述凸优化问题可得到最优化值

步骤六：通过谱峰搜索获得波达方向估计结果。以x轴为个均匀分布的空间网格点来波方向，y轴为优化值中所包含的元素，画出空间功率谱。寻找空间功率谱上的峰值，并将这些峰值所对应的响应值从大到小排列，取前k个峰值所对应的x轴角度方向，即为波达方向估计结果。

本发明一方面充分利用了互质阵列能够增加doa估计方法自由度的优势，突破了均匀线性阵列自由度受限的瓶颈，通过虚拟阵列等价信号的计算实现了在天线阵元个数一定的条件下估计更多个数的入射信号源；另一方面引入了协方差矩阵稀疏重建的思想，并将其应用于虚拟域以实现空间功率谱的重建和doa估计。

下面结合仿真实例对本发明的效果做进一步的描述。

仿真实例1：采用第一类互质阵列接收入射信号，其参数选取为m＝3，n＝5，即架构的互质阵列共包含q＝m+n-1＝7个物理天线阵元。假定入射窄带信号个数为7，且入射方向均匀分布于-60°至60°这一空间角度范围内；波达方向角的角度域范围为[-90°,90°]，空间域网格点均匀采样间距设置为0.1°；正则化参数ξ设置为0.25；信噪比设置为0db，采样快拍数k＝500。

本发明所提出的基于虚拟阵列协方差矩阵稀疏重建的互质阵列波达方向估计方法空间功率谱如图6所示，其中垂直虚线代表入射信号源的实际方向。在本实例的参数设置下，虚拟阵列上连续的虚拟阵元个数为l＝15，相应地，lv＝7。可以看出，本发明所提方法能够有效分辨这7个入射信号源；此外，功率谱上各信号方向的响应值能够体现实际信号源的功率信息，可同时估计各信号的波达方向信息和功率信息。相比于传统采用均匀线性阵列的方法，利用7个物理天线阵元最多只能分辨6个入射信号，以上结果体现了本发明所提方法在自由度性能上的增加。

仿真实例2：采用第二类互质阵列接收入射信号，其参数选取为m＝3，n＝5，即架构的互质阵列共包含q＝2m+n-1＝10个物理天线阵元；假定入射窄带信号个数为15，其余参数设置与仿真实例1保持一致。此时，虚拟阵列上连续的虚拟阵元个数为l＝35，相应地，lv＝17。由图7所示的空间功率谱可以看出，本发明所提方法仅采用10个物理天线阵元就能够有效地分辨15个入射信号源的波达方向和角度信息，体现了自由度性能上的优势。

综上所述，本发明所提方法能够在信号源个数大于等于物理天线个数的情况下实现入射信号的有效估计，增加了自由度和计算效率。此外，与传统采用均匀线性阵列的方法相比，采用本发明所提方法在实际应用系统中所需的物理天线阵元及射频模块也能够相应减少，体现了经济性和高效性。