CN110850359B - 一种基于原子范数的互质阵欠定测向方法 - Google Patents

Abstract

本发明公开了一种基于原子范数的互质阵欠定测向方法，本发明先接收端天线按照互质阵列进行架构；根据阵列接收信号进行建模；计算虚拟阵列输出信号；初始化内插输出信号；设计基于原子范数的优化方法恢复内插输出信号；重建内插后虚拟均匀线性阵列信号协方差矩阵；根据重建协方差矩阵进行波达方向估计。本发明通过内插将互质阵的不连续虚拟线阵插值到虚拟均匀性阵，充分了利用了互质阵的自由度和阵元信息；在向量化协方差矩阵得到虚拟阵列输出信号时，对同一位置的信号进行了平均，从而比直接选取同一位置的第一个信号更能充分利用信号的信息；可以避免传统稀疏表示方法的基不匹配问题，并且在优化问题时，引入了去噪操作，提高了测向的准确性。

Description

一种基于原子范数的互质阵欠定测向方法

技术领域

本发明属于阵列信号处理技术领域，尤其涉及对于入射信号的波达方向的测定，具体是一种基于原子范数的互质阵欠定测向方法。

背景技术

测向，即波达方向(Direction-of-Arrival，DOA)估计，是指利用阵列接收天线接收空域信号，并通过一系列信号处理技术及优化方法估计信号的来波方向，在雷达、语音、无线通信领域等军用和民用领域有着重要的应用。

测向的自由度是指阵列可有效利用的阵元数。传统DOA估计方法阵列自由度受物理阵元数限制。当信号源数大于阵元数时(欠定)，传统方法失效。为了解决欠定测向问题，稀疏阵列结构被提出。在稀疏阵列结构中，互质阵由于具有可解析的阵列结构和较小的阵元间互耦，受到了人们的广泛关注。

现有的基于互质阵列的测向方法主要通过利用质数的性质将互质阵列推广到虚拟域，并形成虚拟阵列接收信号，在此基础上进行测向。由于虚拟阵列中包含的虚拟阵元的个数大于实际阵元的个数，因此自由度得到了有效的提升。但是由互质阵列推导来的虚拟阵列属于非均匀阵列，现有的基于均匀线性阵列的信号处理方法仅能使用虚拟线阵的一段均匀子阵进行DOA估计，从而造成阵列自由度和阵元信息的浪费，如SS-MUSIC方法(见文献：Pal P,Vaidyanathan P P.Coprime sampling and the music algorithm[C].2011Digital Signal Processing and Signal Processing Education Meeting(DSP/SPE).IEEE,2011)和LRD方法(见文献：Pal P,Vaidyanathan P P.A Grid-Less Approachto Underdetermined Direction of Arrival Estimation Via Low Rank MatrixDenoising[J].IEEE Signal Processing Letters.vol.21,no.6,pp.737-741,2014)。为了克服该问题，人们根据压缩感知理论又提出了基于稀疏表示类的方法，如LASSO算法(见文献：Y.D.Zhang,M.G.Amin，B.Himed.Sparsity-based DOA estimation using co-primearrays[C].IEEE Int.Conf.Acoust.,Speech,Signal Process..IEEE,2013.)。虽然此方法提高了阵列自由度的利用率，但由于压缩感知中的稀疏字典是离散的，存在会网格失配问题。为了解决该问题，研究人员提出了基于无网格稀疏的方法，如NMM算法(见文献：C.Liu,Pal P,Vaidyanathan P P..Coprime coarray interpolation for DOA estimation vianuclear norm minimization[C].IEEE Int.Symp.Circuits Syst.(ISCAS).IEEE,2016.)和CA-VAI算法(见文献：C.Zhou,Y.Gu,X.Fan,Z.Shi,G.Mao,Y.D.Zhang.Direction-of-Arrival Estimation for Coprime Array via Virtual Array Interpolation[J].IEEETrans.Signal Process.,vol.66,no.22,pp.5956-5971,2018.)。虽然这两种算法均可以实现无网格的欠定测向，但是在NMM算法中，对缺失的虚拟阵元位置处的虚拟信号进行插值时，并未对噪声进行处理。而CA-VAI算法在得到虚拟阵列输出信号时丢弃了协方差矩阵中的一些元素，未充分利用信号所提供的全部信息，而且此方法是对协方差矩阵进行内插，形成了一个多测量矢量的优化方法，导致了计算复杂度的增加。因此，对于如何使用无网格化的方法提高欠定测向的性能成为一个值得关注的问题。

发明内容

本发明的目的在于针对上述现有技术的不足，提出了一种基于原子范数的互质阵欠定测向方法，充分利用了虚拟阵列所提供的全部信息和自由度并引入去噪操作，从而提高测向性能。

本发明是通过以下技术方案来实现的：一种基于原子范数的互质阵欠定测向方法，包含以下步骤：

(1)接收天线使用N+2M-1个阵元，并按照扩展互质阵列进行架构，其中N和M为互质整数；

(2)假设入射信号为K个波达方向为θ＝[θ₁,θ₂,...,θ_K]^T的远场窄带非相干的平面波信号，入射到互质阵列S上，则t时刻阵列的接收信号可建模为：

其中，

为阵列流型矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为t时刻的接收信号，n(t)＝[n₁(t),n₂(t),...,n_K(t)]^T为与入射信号独立的加性高斯白噪声分量，(·)^T表示转置操作。a(θ_k)表示第k个信源的导引向量，可表示为：

其中，

为虚部单位，u_k表示实际阵元的位置，且u₁＝0。当采样快拍数为T时，阵列接收信号的协方差矩阵可表示为：

其中，(·)^H表示共轭转置操作；

(3)计算虚拟阵列输出信号：将阵列接收信号的协方差矩阵

向量化可以得到：

其中，p_k为第k个入射信号的功率，

为克罗内克积(Kronecker product)，(·)^*为取共轭操作，vec(·)表示向量化算子，且i＝vec(I)。虚拟阵元的位置可以表示为

其中：

且

|·|表示集合的势。由于互质阵列是一种部分可扩展阵列，因此其对应的虚拟阵列

是不连续的，即存在一些“孔洞”。通过将信号y中对应于

的相同位置的信号做平均后即可得到虚拟阵列输出信号

其中，

为虚拟阵列流型矩阵，

表示中心元素为1，其余元素全为0的向量；

(4)初始化内插输出信号：由于虚拟阵列具有一些“孔洞”，且实际上“孔洞”位置处的输出信号是不存在的。可以通过假设在“孔洞”位置处有虚构阵元，而其对应的输出信号可以通过内插得到。因此，将虚构阵元与虚拟阵元结合起来就可以得到一个阵元位置为

的虚拟均匀线性阵列，其中

首先将虚构阵元出的输出信号初始化为0，则虚拟均匀线性阵列的初始化内插输出信号可以表示为：

其中，<z〉_i表示对应于阵元位置为id处的信号；

(5)设计基于原子范数的优化算法恢复内插输出信号：首先构建表示信号z的原子集

则内插输出信号的原子范数

可以定义为在原子集

中能够表示信号z的最小原子数，即：

其中，inf表示下确界。因此，最优值z^*可以通过下面的原子范数最小化问题得到：

其中，η²表示噪声上限，Ω表示由物理阵元推导的非连续虚拟阵元位置的集合。将上述原子范数

转换为以下半定规划(SDP)形式，就可以使用SDP求解器高效地求解：

其中，Tr(·)表示矩阵的迹，

表示

作为第一列的Toeplitz矩阵。因此，该原子范数

最小化问题可以通过下式计算：

其中，

t为调节参数。进而基于该原子范数的优化问题可以表示为：

其中，

将上述优化问题转化为对偶问题即可得到内插输出信号z^*；

(6)重建内插后虚拟均匀线性阵列信号协方差矩阵：由于内插虚拟信号是秩为1的信号。因此，该问题等价于相干信号或单快拍信号的DOA估计问题。令

重建的信号协方差矩阵

可通过下式得到：

(7)根据重建协方差矩阵

进行波达方向估计。

进一步地，步骤(1)所述的扩展互质阵列结构可具体描述为：首先选取一对互质整数N和M，然后构造一对稀疏的均匀线性子阵，其中第一个子阵列包含2M个间距为Nd的天线阵列，其位置为0，Nd，…，(2M-1)Nd，第二个子阵列包含N个间距为Md的天线阵列，其位置为0，Md，…，(N-1)Md，单位间距取接收信号的半波长，即d＝λ/2。然后将两个子阵列按照首个阵元重叠的方式进行子阵列组合，即可获得实际包含有N+2M-1个天线阵列的扩展互质阵列架构。

进一步地，步骤(5)中的优化问题可以将其转化为对偶问题，并利用SDP的标准求解器SDPT3更高效地求解。根据拉格朗日分析法，原问题的对偶问题可以表示为如下形式：

其中，η表示正则化参数，

表示取实部操作，||·||₂表示l₂范数，T^*(·)表示T(·)的伴随算子。根据对偶性，对应于对偶问题中的第一个约束条件的拉格朗日乘子就是原问题中的第一个约束矩阵，其中包含了内插输出信号z。因此在解决对偶问题后，即可获得内插输出信号z^*。

进一步地，步骤(7)还可采用多重信号分类算法(MUSIC)及其改进算法、基于旋转不变技术估计信号参数(ESPRIT)算法、加权子空间拟合(WSF)算法等多种方法进行波达方向估计。

进一步地，步骤(7)中，通过多重信号分类方法进行波达方向估计，具体为计算虚拟域空间谱：

其中，a^H(θ)为虚拟接收信号的导引矢量，对应于位置为由0到L的一段虚拟均匀阵列，U_N是L×(L-K)维的噪声子空间；通过谱峰搜索寻找空间谱P_MUSIC的峰值即为波达方向估计结果。

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

首先，本发明采用互质阵列进行波达方向角度估计，克服了采用均匀线阵造成的自由度受限于物理阵元数的缺点，并且通过内插将不连续虚拟阵列插值到虚拟均匀线性阵列，充分了利用了互质阵的自由度和阵元信息；其次，本发明在向量化协方差矩阵得到虚拟阵列输出信号时，对同一位置的信号进行了平均，从而比直接选取同一位置的第一个信号更能充分利用信号的信息；最后，本发明的优化方法属于无网格方法，可以避免传统稀疏表示算法的基不匹配问题，并且在优化问题时，引入了去噪操作，提高了测向的准确性。

附图说明

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

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

图3是本发明中扩展互质阵列的结构示意图。

图4是用于体现本发明所提方法自由度性能的空间功率谱示意图。

图5是本发明所提方法与现有各方法均方根误差与信噪比之间的对比结果，采用快拍数为500。

图6是本发明所提方法与现有各方法均方根误差与快拍数之间的对比结果，信噪比为20dB。

具体实施方式

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

现有的互质阵列波达方向估计方法通常通过计算和处理虚拟阵列输出信号，以突破物理阵元数对于自由度的限制，实现自由度性能的提升。但是由于虚拟阵列存在“孔洞”，使得所有的阵元不能被充分利用，一种普遍的做法是采用虚拟阵列中的一段均匀阵列进行处理，从而导致了丢弃阵元上信息的缺失。为了充分利用互质阵列接收信号所包含的全部信息，本发明提供了一种基于原子范数的互质阵欠定测向方法，参照图1，本发明的实现步骤如下：

步骤一：在接收端使用N+2M-1个天线阵元架构扩展互质阵列。首先选取一组互质整数M和N；然后构造一对稀疏均匀线性子阵列，其中第一个子阵列包含2M个阵元间距为Nd的天线阵元，其位置为0，Nd，…，(2M-1)Nd；第二个子阵列包含N个阵元间距为Md的天线阵元，其位置为0，Md，…，(N-1)Md；单位间距取接收信号的半波长，即d＝λ/2。然后将两个子阵列的首个天线阵元视作参考阵元，参照图2和图3，将两个子阵的参考阵元重叠以实现子阵列组合，即可获得实际包含N+2M-1个天线阵元的非均匀扩展互质阵列架构。

步骤二：利用扩展互质阵列接收入射信号并建模。假设入射信号为K个波达方向为θ＝[θ₁,θ₂,…,θ_K]^T的远场窄带非相干的平面波信号，入射到互质阵列S上，则t时刻阵列的接收信号可建模为：

其中，

为阵列流型矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为t时刻的接收信号，n(t)＝[n₁(t),n₂(t),...,n_K(t)]^T为与入射信号独立的加性高斯白噪声分量，(·)^T为转置操作。a(θ_k)表示第k个信源的互质阵列导引向量，可表示为：

其中，

为虚部单位，u_k表示实际阵元的位置，且u₁＝0。当采样快拍数为T时，互质阵列接收信号的协方差矩阵可表示为：

其中，(·)^H表示共轭转置操作。

步骤三：计算虚拟阵列输出信号。将阵列接收信号的协方差矩阵

向量化可以得到：

其中，p_k为第k个入射信号的功率，

为克罗内克积(Kronecker product)，(·)^*为取共轭操作，vec(·)表示向量化算子，且i＝vec(I)。虚拟阵元的位置可以表示为

其中：

且

|·|表示集合的势。由于互质阵列是一种部分可扩展阵列，因此其对应的虚拟阵列

是不连续的，即存在一些“孔洞”。通过将信号y中对应于

的相同位置的信号做平均后即可得到虚拟阵列输出信号

其中，

为虚拟阵列流型矩阵，

表示中心元素为1，其余元素全为0的向量。

步骤四：初始化内插输出信号。由于虚拟阵列具有一些“孔洞”，且实际上“孔洞”位置处的输出信号是不存在的。可以通过假设在“孔洞”位置处有虚构阵元，而其对应的输出信号可以通过内插得到。因此，将虚构阵元与虚拟阵元结合起来就可以得到一个阵元位置为

的虚拟均匀线性阵列，其中

首先将虚构阵元处的输出信号初始化为0，则虚拟均匀线性阵列的初始化内插输出信号可以表示为：

其中，<z>_i表示对应于阵元位置为id处的信号。

步骤五：设计基于原子范数的优化算法恢复内插输出信号。首先构建表示信号z的原子集

则内插输出信号的原子范数

可以定义为在原子集

中能够表示信号z的最小原子数，即：

其中，inf表示下确界。因此，最优值z^*可以通过下面的原子范数最小化问题得到：

其中，η²表示噪声上限，Ω表示由物理阵元推导的非连续虚拟阵元位置的集合。将上述原子范数

转换为以下半定规划(SDP)形式，就可以使用SDP求解器高效地求解：

其中，Tr(·)表示矩阵的迹，

表示

作为第一列的Toeplitz矩阵。因此，该原子范数

最小化问题可以通过下式计算：

其中，

t为调节参数。进而基于该原子范数的优化问题可以表示为：

其中，

通过将其转化为对偶问题，可以利用SDP的标准求解器SDPT3更加高效地求解。根据拉格朗日分析法，上式的对偶问题可以表示为如下形式：

其中，η表示正则化参数，

表示取实部操作，||·||₂表示l₂范数，T^*(·)表示T(·)的伴随算子。根据对偶性，对应于对偶问题中的第一个约束条件的拉格朗日乘子就是原问题中的第一个约束矩阵，其中包含了内插输出信号z。因此在解决对偶问题后，即可获得内插输出信号z^*。

步骤六：重建内插后虚拟均匀线性阵列信号协方差矩阵。由于内插虚拟信号是秩为1的信号。因此，该问题等价于相干信号或单快拍信号的DOA估计问题。尽管空间平滑可以去除信号的相干性，但是还有一种更加有效的方法可以重构秩为K的信号协方差矩阵。令

重建的信号协方差矩阵

可通过下式得到：

步骤七：根据重建的内插输出信号协方差矩阵

进行波达方向估计。得到重建的内插输出信号协方差矩阵后，可以通过经典的波达方向估计方法，如多重信号分类(MUSIC)算法及其改进算法、基于旋转不变技术估计信号参数(ESPRIT)算法、加权子空间拟合(WSF)算法等多种方法均求得波达方向估计结果。以多重信号分类方法为例，其空间功率谱P_MUSIC可表示为：

其中，a^H(θ)为虚拟接收信号的导引矢量，对应于位置为由0到L的一段虚拟均匀阵列，U_N是L×(L-K)维的噪声子空间；通过谱峰搜索寻找空间谱P_MUSIC的峰值即为波达方向估计结果。

本发明采用扩展互质阵列进行波达方向角度估计，克服了采用均匀线性阵列造成的自由度受限于物理阵元数的缺点，而且充分利用了虚拟阵列提供的全部信息和自由度，避免了由于虚拟阵列非均匀性造成的信息损失，有效提高了DOA估计的性能。

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

仿真实例1：选取M＝3，N＝5，即采用含有N+2M-1＝10个物理阵元的扩展互质阵列接收入射信号，阵列孔径为25d；假设有21个不相关的信号源均匀分布在[25.84°，154.15°]的角度范围内。信噪比设置为20dB，采样快拍数为500，且正则化参数设置为1。本发明提出的基于原子范数的互质阵欠定测向方法的空间功率谱如图4所示，图中垂直虚线的位置为入射信源的实际方向。可以看出，本发明所提方法能够有效分辨这21个信源。传统的采用均匀线阵波达方向估计方法，采用10个物理阵元只能分辨9个入射信号，而只采用扩展互质阵列的一段虚拟均匀线性子阵列的自由度为MN+M-1，即MN+M-1＝17＜21。以上结果表明了本发明所提方法实现了互质阵列自由度的提高。

仿真实例2：采用扩展互质阵列接收入射信号，其参数同样选取M＝3，N＝5，即架构的扩展互质阵列含有N+2M-1＝10个物理天线阵元；假设12个远场窄带平面波信号的入射方向均匀分布于25°到155°之间，且在每个入射方向上加入随机分布于-0.1°-0.1°偏移量，以便更好展现算法的性能。本发明所提的算法与LASSO算法、SS-MUSIC算法，低秩矩阵降噪(LRD)算法、核范数最小化(NNM)算法和基于虚拟阵元内插(CA-VAI)算法等一些现存的使用互质阵列的算法进行了对比。LASSO算法、LRD算法、CA-VAI算法和本发明所提算法中的正则化参数分别为1.5、10、0.25和1，其中LASSO算法的网格密度设置为Δθ＝0.1°。各方法的均方根误差随信噪比和采样快拍数之间的变化如图5、图6，其中图5采用的快拍数为500，图6采用的信噪比为20dB。对于多重信号分类方法的空间谱搜索范围为[0°,180°]，搜索间距为0.01°。与此同时，克拉美罗界也在图5和图6中给出，用于指示估计性能的下界。针对每一组仿真参数的设置，均方根误差均通过1000次蒙特卡洛试验求平均获得。

由图5和图6的均方根误差性能对比示意图可以看出，本发明多提算法性能优于其他算法性能，体现了所提算法估计的准确度。

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

Claims (5)

1.一种基于原子范数的互质阵欠定测向方法，其特征在于，该方法包括以下步骤：

(1)接收天线使用N+2M-1个阵元，并按照扩展互质阵列进行架构，其中N和M为互质整数；

(2)假设入射信号为K个波达方向为θ＝[θ₁,θ₂,...,θ_K]^T的远场窄带非相干的平面波信号，入射到互质阵列

上，则t时刻阵列的接收信号建模为：

其中，

为阵列流型矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为t时刻的接收信号，n(t)＝[n₁(t),n₂(t),...,n_K(t)]^T为与入射信号独立的加性高斯白噪声分量，(·)^T表示转置操作；a(θ_k)表示第k个信源的导引向量，可表示为：

其中，

为虚部单位，u_k表示实际阵元的位置，且u₁＝0；当采样快拍数为T时，阵列接收信号的协方差矩阵可表示为：

其中，(·)^H表示共轭转置操作；

(3)计算虚拟阵列输出信号：将阵列接收信号的协方差矩阵

向量化得到：

其中，p_k为第k个入射信号的功率，

为克罗内克积，(·)^*为取共轭操作，vec(·)表示向量化算子，且i＝vec(I)；虚拟阵元的位置表示为

其中：

且

|·|表示集合的势；由于互质阵列是一种部分可扩展阵列，因此其对应的虚拟阵列

是不连续的，即存在一些“孔洞”；通过将信号y中对应于

的相同位置的信号做平均后即可得到虚拟阵列输出信号

其中，

为虚拟阵列流型矩阵，

表示中心元素为1，其余元素全为0的向量；

(4)初始化内插输出信号：由于虚拟阵列具有一些“孔洞”，且实际上“孔洞”位置处的输出信号是不存在的；通过假设在“孔洞”位置处有虚构阵元，而其对应的输出信号通过内插得到；因此，将虚构阵元与虚拟阵元结合起来得到一个阵元位置为

的虚拟均匀线性阵列，其中

首先将虚构阵元处的输出信号初始化为0，则虚拟均匀线性阵列的初始化内插输出信号表示为：

其中，<z>_i表示对应于阵元位置为id处的信号；

(5)设计基于原子范数的优化方法恢复内插输出信号：首先构建表示信号z的原子集

则内插输出信号的原子范数

定义为在原子集

中能够表示信号z的最小原子数，即：

其中，inf表示下确界；因此，最优值z^*通过下面的原子范数最小化问题得到：

其中，η²表示噪声上限，Ω表示由物理阵元推导的非连续虚拟阵元位置的集合；将上述原子范数

转换为以下半定规划形式，使用SDP求解器高效地求解：

其中，Tr(·)表示矩阵的迹，

表示

作为第一列的Toeplitz矩阵；因此，该原子范数

最小化问题通过下式计算：

其中，

t为调节参数；进而基于该原子范数的优化问题表示为：

其中，

将上述优化问题转化为对偶问题即可得到内插输出信号z^*；

(6)重建内插后虚拟均匀线性阵列信号协方差矩阵：由于内插虚拟信号是秩为1的信号；因此，该问题等价于相干信号或单快拍信号的DOA估计问题；令

重建的信号协方差矩阵

通过下式得到：

(7)根据重建协方差矩阵

进行波达方向估计。

2.根据权利要求1所述的一种基于原子范数的互质阵欠定测向方法，其中步骤(1)所述的扩展互质阵列结构可具体描述为：首先选取一对互质整数N和M，然后构造一对稀疏的均匀线性子阵，其中第一个子阵列包含2M个间距为Nd的天线阵列，其位置为0，Nd，…，(2M-1)Nd，第二个子阵列包含N个间距为Md的天线阵列，其位置为0，Md，…，(N-1)Md，单位间距取接收信号的半波长，即d＝λ/2；然后将两个子阵列按照首个阵元重叠的方式进行子阵列组合，即可获得实际包含有N+2M-1个天线阵列的扩展互质阵列架构。

3.根据权利要求1所述的一种基于原子范数的互质阵欠定测向方法，导出原子范数

的半正定规划形式，将内插输出信号原子范数最小化的优化问题转化为对偶问题，表示为如下形式：

其中，η表示正则化参数，

表示取实部操作，||·||₂表示l₂范数，T^*(·)表示T(·)的伴随算子；根据对偶性，对应于对偶问题中的第一个约束条件的拉格朗日乘子就是原问题中的第一个约束矩阵，其中包含了内插输出信号z；因此在解决对偶问题后，即可获得内插输出信号z^*。

4.根据权利要求1所述的一种基于原子范数的互质阵欠定测向方法，其中步骤(7)还采用多重信号分类算法及其改进算法、基于旋转不变技术估计信号参数算法、加权子空间拟合算法进行波达方向估计。

5.根据权利要求1所述的一种基于原子范数的互质阵欠定测向方法，其中步骤(7)中，通过多重信号分类方法进行波达方向估计，具体为计算虚拟域空间谱：

其中，a^H(θ)为虚拟接收信号的导引矢量，对应于位置为由0到L的一段虚拟均匀阵列，U_N是L×(L-K)维的噪声子空间；通过谱峰搜索寻找空间谱P_MUSIC的峰值即为波达方向估计结果。

Images