基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法与流程

本发明属于信号处理技术领域，尤其涉及对雷达信号、声学信号及电磁信号的波达方向估计，具体是一种基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法，可用于无源定位和目标探测。

背景技术

波达方向(direction-of-arrival,doa)估计是阵列信号处理领域的基本问题之一，指的是利用阵列天线接收空域信号，并通过一系列统计信号处理技术及优化方法实现对接收信号的有效处理，从而得到信号的doa信息，在雷达、声呐、语音、无线通信等领域均有着重要而广泛的应用。

自由度是doa估计方法中的重要指标之一，指的是其能够估计的入射信号的个数。目前的doa估计方法通常基于均匀线性阵列进行信号的接收和建模，并围绕该模型形成了成熟的理论和广泛的应用。然而，基于均匀线阵的doa估计方法自由度受到实际天线阵元个数的限制，当空间中入射信号源的个数大于或等于阵列中物理天线阵元个数时，现有基于均匀线阵的doa估计方法将无法进行有效的估计。

互质阵列是一种具有系统化结构的非均匀稀疏阵列，能够在天线阵元个数一定的情况下增加doa估计的自由度，突破了传统均匀线阵自由度受限的瓶颈，实现了doa估计方法自由度性能的提升，受到了学术界的广泛关注。现有的基于互质阵列的doa估计方法主要利用质数的性质，对互质阵列进行推导以形成一个二阶等价虚拟均匀线性阵列接收信号，进而实现doa估计。由于虚拟阵列中的虚拟阵元数大于实际的天线阵元数，因此doa估计的自由度能够得到有效的提升。然而，从互质阵列直接推导而来的虚拟阵列是非均匀的，为了避免非均匀阵列引起的信号模型失配问题，现有的互质阵列doa估计方法大多简单地利用虚拟阵列中连续阵元部分形成虚拟均匀线阵进行doa估计，这引起了部分原始信息的丢失以及相关估计性能的降低。

此外，目前众多doa估计方法均包含网格化的设计步骤，即在空间谱峰搜索过程中设置一定的搜索步长或者在优化问题中引入预先定义的波达方向过完备空间网格点。此类设计造成了doa估计精度和计算复杂度之间的矛盾，即随着doa估计精度要求的提高，网格化doa估计方法需要设定更加密集的网格，导致计算复杂度的急剧增加。不仅如此，在实际中，即使网格点较为密集，在很多情况下信号的波达方向无法完全落在预先设定的网格点上，这种基不匹配现象也导致了固有的波达方向估计误差。

技术实现要素：

本发明的目的在于针对上述现有技术存在的不足，提出一种基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法。一方面，通过虚拟阵列内插技术充分利用了非均匀虚拟阵列所包含的全部信号信息，从而进一步提高了doa估计的自由度；另一方面，利用内插虚拟阵列协方差矩阵子空间的旋转不变性，给出了波达方向估计的闭式解，从而实现了无网格化的高精度doa估计。

本发明的目的是通过以下技术方案来实现的：一种基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法，包含以下步骤：

(1)接收端使用m+n-1个天线，并按照互质阵列结构进行架构；其中m与n为互质整数；

(2)互质阵列接收信号建模并计算其等价虚拟信号：假设有k个来自θ1,θ2,…，θk方向的远场窄带非相关信号源，采用步骤(1)架构的非均匀互质阵列接收入射信号，得到(m+n-1)×1维互质阵列接收信号x(t)，可建模为：

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

其中，zi,i＝1,2,…,m+n-1,表示互质阵列中第i个物理天线阵元的实际位置，且z1＝0；j为虚数单位，[·]^t表示转置操作，为互质阵列各物理阵元位置的集合，可表示为：

其中，λ为入射信号的波长。共采集t个采样快拍，得到互质阵列接收信号的采样协方差矩阵

其中，(·)^h表示共轭转置操作。向量化互质阵列接收信号的采样协方差矩阵获得虚拟阵列等价接收信号

其中，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，为向量对应的非均匀虚拟阵列中各虚拟阵元位置的集合，可表示为：

(3)构造内插虚拟阵列并初始化其虚拟域信号统计量：对于非均匀的虚拟阵列在保留其原有虚拟阵元位置不变的前提下，向其中存在孔洞的位置插入若干虚拟阵元，从而将非均匀虚拟阵列转化为间距为d、阵列孔径与相同、且虚拟阵元数目增加的均匀虚拟阵列

其中，表示集合中的最小值，表示集合中的最大值；内插虚拟阵列共包含个虚拟阵元，其中|·|表示集合的势。内插虚拟阵列对应的等价虚拟信号可通过在向量中与孔洞相对应的位置填充0获得，可表示为：

其中，表示向量中对应于虚拟阵元位置为的元素。若中存在多个不同元素与同一虚拟阵元位置对应的情况，则计算这些元素的平均值作为对应的内插虚拟阵列的初始化采样协方差矩阵可以构建为：

其中，为集合中非负元素组成的子集，为的子向量，其中的元素为与中虚拟阵元位置对应的中的元素，表示以为第一列的厄米特对称toeplitz矩阵；

(4)构造投影矩阵并定义投影运算：定义一个与维度相同的投影矩阵ω，如果中某一位置上的元素是0，则投影矩阵ω中相同位置的元素值也为0；反之则投影矩阵ω中相应位置的元素值为1。定义为投影运算，其中括号内变量为与ω维度相同的矩阵，投影运算通过变量矩阵的每一个元素与投影矩阵ω中相应位置上的元素一一相乘实现，得到一个与矩阵ω维度相同的矩阵；

(5)设计基于核范数最小化的优化问题并求解，以重建内插虚拟阵列协方差矩阵：构建如下以向量u为变量的优化问题：

其中，表示的核范数，表示以向量u为第一列的厄米特对称toeplitz矩阵；为阈值常数，用于约束协方差矩阵的重建误差；保证了重建的协方差矩阵满足半正定的条件；‖·‖f表示frobenius范数。求解上述凸优化问题可得到最优化解相应地，toeplitz矩阵为重建的内插虚拟阵列协方差矩阵；

(6)形成重建的内插虚拟阵列协方差矩阵的信号子空间：对进行特征值分解：

其中，λs为k×k维对角矩阵，其对角线上包含将的特征值从大到小排列后的前k个特征值，为包含该k个特征值所对应特征向量的维矩阵；相应地，λn为维对角矩阵，其对角线上包含剩余的个的特征值，为包含剩余的个特征值所对应特征向量的维矩阵；为的信号子空间；

(7)构造具有旋转不变性的两个虚拟子阵信号子空间：将虚拟均匀线性阵列分为两个虚拟子阵和其中，为包含中位置为0到的一段虚拟均匀线性子阵列，为包含中位置为d到的一段虚拟均匀线性子阵列；将矩阵分别去除最后一行和第一行后得到两个维矩阵和以形成这两个虚拟子阵和对应的信号子空间，该两个信号子空间具有旋转不变性；

(8)计算两个虚拟子阵信号子空间的转换矩阵γ：两个虚拟子阵和的信号子空间和之间的k×k维转换矩阵γ可由下式得到：

其中，表示求伪逆运算；

(9)形成波达方向估计的闭式解：信号的波达方向估计值可计算为：

其中，imag(·表示复数的虚部，γk是矩阵γ的第k个特征值。

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

进一步地，步骤(5)中的凸优化问题可转化为如下以向量u为变量的优化问题：

其中μ为正则化参数，用于在最小化过程中权衡矩阵的重建误差和的核范数。

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

(1)本发明通过对互质阵列推导得到的非均匀虚拟阵列进行内插，不仅保留了原始非均匀虚拟阵列所包含的全部接收信号信息，而且此方法构建的虚拟均匀线性阵列等价虚拟信号满足奈奎斯特采样定律；

(2)本发明基于内插虚拟阵列协方差矩阵核范数最小化的思想设计矩阵重建优化问题并求解，以获得内插虚拟阵列协方差矩阵，并保证了优化求解结果为厄米特对称的半正定toeplitz矩阵。由于均匀线性阵列非相关接收信号的理论协方差矩阵满足toeplitz结构，因此利用它的toeplitz特性作为先验条件进行协方差矩阵的重建，可以使重建结果与真实值差异更小，从而提高doa估计的性能；

(3)本发明利用重建所得内插虚拟阵列协方差矩阵的子空间旋转不变性，实现了无网格化的波达方向估计，避免了现有网格化波达方向估计方法中的谱峰搜索以及预定义空间网格点等设计步骤引起的估计精度与计算复杂度之间的矛盾；通过计算波达方向估计的闭式解，避免了网格化波达方向估计方法中基不匹配现象带来的固有估计误差，同时保证了波达方向估计的准确度和计算效率。

附图说明

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

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

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

图4是本发明中内插虚拟阵列的结构示意图。

图5是用于体现本发明所提方法自由度性能的空间谱示意图，假设入射信号源个数为9

图6是用于体现本发明所提方法自由度性能的空间谱示意图，假设入射信号源个数为11

图7是用于体现本发明所提方法doa估计准确度的均方根误差性能对比示意图

具体实施方式

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

对于doa估计在实际系统中的应用，互质阵列由于可以通过等价虚拟阵列信号的计算和统计信号处理，突破物理阵元数量对自由度的限制而备受关注。但是受限于虚拟阵列的非均匀性，目前很多方法都会选择利用其中连续部分的虚拟阵元进行doa估计，从而造成了信息损失。同时，很多方法在进行doa估计之前会预先设置假定的信号波达方向空间网格点，这造成了固有的估计误差以及计算复杂度与估计准确度之间的矛盾。为了充分利用非均匀虚拟阵列中所包含的所有信号信息，并避免由预定义空间网格点所造成的估计准确度受限问题，本发明提供了一种基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法，参照图1，本发明的实现步骤如下：

步骤一：在接收端使用m+n-1个天线阵元架构互质阵列。首先，选取一组互质整数m、n；然后，参照图2，构造一对稀疏均匀线性子阵列，其中第一个子阵列包含m个间距为nd的天线阵元，其位置为0,nd,…,(m-1)nd；第二个子阵列包含n个间距为md的天线阵元，其位置为0,md,…,(n-1)md；单位间距d取为入射窄带信号波长λ的一半，即d＝λ/2；接着，将两个子阵列的首个天线阵元视为参考阵元，参照图3，将两个子阵的参考阵元重叠以实现子阵列组合，获得实际包含m+n-1个天线阵元的非均匀互质阵列架构，各阵元位置的集合定义为可表示为

步骤二：互质阵列接收信号建模并计算其等价虚拟信号。假设有k个来自θ1,θ2,…,θk方向的远场窄带非相关信号源，采用步骤一架构的非均匀互质阵列接收入射信号，得到(m+n-1)×1维互质阵列接收信号x(t)，可建模为：

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

其中，zi,i＝1,2,…,m+n-1,表示互质阵列中第i个物理天线阵元的实际位置，且z1＝0；j为虚数单位，[·]^t表示转置操作。共采集t个采样快拍，得到互质阵列接收信号的采样协方差矩阵

其中，(·)^h表示共轭转置操作。向量化互质阵列接收信号的采样协方差矩阵获得虚拟阵列等价接收信号

其中，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，为向量对应的非均匀虚拟阵列中各虚拟阵元位置的集合，可表示为：

步骤三：构造内插虚拟阵列并初始化其虚拟域信号统计量。参照图4，对于非均匀的虚拟阵列在保留其原有虚拟阵元位置不变的前提下，向其中存在孔洞的位置插入若干虚拟阵元(如图4中的空心圆所示)，从而将非均匀虚拟阵列转化为间距为d、阵列孔径与相同、且虚拟阵元数目增加的虚拟均匀线性阵列

其中，表示集合中的最小值，表示集合中的最大值。内插虚拟阵列共包含个虚拟阵元，其中|·|表示集合的势。内插虚拟阵列对应的等价虚拟信号可通过在向量中孔洞的相应位置填充0获得，可表示为：

其中，表示向量中对应于虚拟阵元位置为的元素。若中存在多个不同元素与同一虚拟阵元位置对应，则计算这些元素的平均值作为对应的根据和toeplitz矩阵结构的关系，内插虚拟阵列的初始化采样协方差矩阵可构建为：

其中，为集合中非负元素组成的子集，为的子向量，其中的元素为与中各虚拟阵元位置对应的中的元素，表示以为第一列的厄米特对称toeplitz矩阵。

步骤四：构造投影矩阵并定义投影运算。根据步骤三，由于中包含初始化而来的对应于内插虚拟阵元的零元素，所以中相应位置对角线上的元素全部为零。根据这样的结构定义一个与维度相同的投影矩阵ω，如果中某一位置上的元素是0，则投影矩阵ω中相同位置的元素值也为0；反之则投影矩阵ω中相应位置的元素值为1。定义为投影运算，其中括号内变量为与ω维度相同的矩阵，投影运算通过变量矩阵的每一个元素与投影矩阵ω中相应位置上的元素一一相乘实现，得到一个与矩阵ω维度相同的矩阵。

步骤五：设计基于核范数最小化的优化问题并求解，以重建内插虚拟阵列协方差矩阵。利用步骤三得到的初始化内插虚拟阵列采样协方差矩阵作为参考值，寻找一个核范数最小的toeplitz矩阵作为重建的内插虚拟阵列协方差矩阵，且要求其与的拟合误差小于某一阈值，可构建如下以向量u为变量的优化问题：

其中，表示的核范数，表示以向量u为第一列的厄米特对称toeplitz矩阵；为阈值常数，用于约束协方差矩阵的重建误差；保证了重建的协方差矩阵满足半正定的条件；‖·‖f表示frobenius范数。上述凸优化问题可转化为以下以向量u为变量的优化问题：

其中μ为正则化参数，用于在最小化过程中权衡矩阵重建误差和的核范数。求解上述凸优化问题可得到最优化解相应地，toeplitz矩阵为重建的内插虚拟阵列协方差矩阵。

步骤六：形成重建的内插虚拟阵列协方差矩阵的信号子空间。对进行特征值分解：

其中，λs为k×k维对角矩阵，其对角线上包含将的特征值从大到小排列后的前k个特征值，为包含该k个特征值所对应特征向量的维矩阵；相应地，λn为维对角矩阵，其对角线上包含剩余的个的特征值，为包含该个特征值所对应特征向量的维矩阵。其中，的信号子空间与噪声子空间正交，故存在一个唯一的、非奇异的k×k维满秩矩阵t满足其中为虚拟阵列的导引矩阵。

步骤七：构造具有旋转不变性的两个虚拟子阵的信号子空间。将虚拟均匀线性阵列分为两个虚拟子阵，定义为和其中，为包含中位置为0到的一段虚拟均匀线性子阵列，为包含中位置为d到的一段虚拟均匀线性子阵列。于是，这两个虚拟子阵的导引矩阵理论值分别为维矩阵和其中，和可由步骤六中虚拟均匀线性阵列导引矩阵分别移除最后一行和第一行得到。由于和具有完全相同的均匀线性阵列结构，仅由子阵间的单位位移d产生的虚拟子阵移不变性，则两个虚拟子阵的导引矩阵和之间的关系可由一个旋转因子φ表示：

其中，diag{·}表示对角矩阵；相应地，将步骤六中得到的矩阵分别去除最后一行和第一行后得到的两个维矩阵和形成虚拟子阵和的信号子空间，则这两个虚拟子阵的信号子空间具有旋转不变性。

步骤八：计算两个虚拟子阵信号子空间的转换矩阵γ。两个虚拟子阵和的信号子空间和之间的k×k维转换矩阵γ可由下式得到：

其中，表示伪逆运算。

步骤九：计算波达方向估计结果。根据步骤六中和步骤七中的两个虚拟子阵信号子空间的旋转不变性，和分别满足以及故转换矩阵γ满足φ＝tγt^-1，其中(·)^-1表示矩阵的逆。根据步骤七中φ的定义，结合转换矩阵γ的特征值和旋转因子φ中对角元素上所包含的波达方向信息，我们可以得到各入射信号的波达方向估计：

其中，imag(·)表示复数的虚部，γk是γ的第k个特征值。

本发明基于虚拟阵列内插的思想，在推导的非均匀虚拟阵列基础上内插入虚拟阵元从而得到一个虚拟均匀线性阵列，在有效利用原始非均匀虚拟阵列所包含的全部信号信息的同时，避免了因虚拟阵列的非均匀性所引起的信号模型失配问题。引入基于虚拟阵列协方差矩阵核范数最小化的思想设计优化问题，实现了内插虚拟阵列协方差矩阵的重建，并利用其子空间旋转不变性实现了虚拟域中的无网格化波达方向估计。

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

仿真实例1：采用互质阵列接收入射信号，其参数选取为m＝3，n＝5，即架构的互质阵列共包含m+n-1＝7个物理阵元。首先假定入射窄带信号个数为9，且入射方向均匀分布于-50°至50°这一空间角度域范围内；之后假定入射窄带信号个数为11，且入射方向同样均匀分布于-50°至50°这一空间角度域范围内；两次仿真信噪比均设置为0db，采样快拍数均为t＝500；正则化参数μ设置均为0.25。

本发明所提出的基于内插虚拟阵列协方差矩阵子空间旋转不变性的互质阵列波达方向估计方法空间谱如图5和图6所示，其中横轴表示各入射信号源的序号，纵轴表示波达方向角度。从图5可以看出，本发明所提方法能够有效分辨这9个信号源。而对于传统使用均匀线性阵列的方法，利用7个物理阵元最多只能分辨6个入射信号，该结果体现了本发明所提方法实现了波达方向估计自由度的增加。此外，9个信号源也超过了非均匀虚拟阵列中连续部分的虚拟阵元个数，说明了非均匀虚拟阵列中非连续部分的虚拟阵元也得到了有效的利用。从图6可以看出，本发明所提方法能够有效分辨11个入射信号，且该仿真条件下入射信号数量已经超过了非均匀虚拟阵列中的全部虚拟阵元数量，说明了本发明所提方法中引入的虚拟阵列内插技术进一步增加了自由度。

仿真实例2：采用互质阵列接收入射信号，其参数同样选取为m＝3，n＝5，即架构的互质阵列共包含m+n-1＝7个物理天线阵元；假定入射窄带信号个数为1，入射信号波达方向随机生成且满足高斯分布在不同信噪比情况下的doa估计均方根误差对比如图7所示，采样快拍数为500，仿真所示结果中每一个数值点均通过1000次蒙特卡洛试验求平均值获得。本发明所提方法将与基于矩阵填充思想的核范数最小化方法进行对比，与此同时，克拉美罗界也将在图7中同时给出，作为参考。可以看出，随着信噪比增大，本发明所提方法逼近克拉美罗界，充分体现了无网格化doa估计高准确度的优势。

综上所述，本发明所提方法通过虚拟阵列内插技术充分利用了非均匀虚拟阵列中所包含的全部信号信息，能够在信号源个数大于等于物理天线个数甚至大于等于非均匀虚拟阵列阵元个数的情况下实现doa的有效估计，增加了doa估计的自由度。利用内插虚拟阵列协方差矩阵的子空间旋转不变性得到doa估计的闭式解，实现了无网格化的doa估计，保证了doa估计的准确度。此外，与传统采用均匀线性阵列的方法相比，本发明所提方法在实际应用中所需的物理天线阵元及射频模块也能够相应减少，体现了经济性和高效性。