基于信号子空间旋转不变性的互质阵列波达方向估计方法与流程

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

背景技术：

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

doa估计方法的自由度是指其能够分辨的入射信号源个数，作为实际系统应用中的一个重要性能指标，决定着系统的总体复杂度。现有的doa估计方法通常采用均匀线性阵列进行信号的接收与建模，但是基于均匀线性阵列方法的自由度受限于实际天线阵元个数。具体而言，对于一个包含l个天线阵元的均匀线性阵列，其自由度为l-1，即最多只能分辨l-1个入射信号。因此，当某个空域范围内入射信号源的个数大于或等于阵列中天线阵元的个数时，现有采用均匀线性阵列的方法将无法进行有效的doa估计。为了增加自由度，传统方法需要通过增加物理天线阵元及相应的射频模块来实现，这造成了系统计算复杂度和硬件复杂度的增加。因此，现有采用均匀线性阵列的doa估计方法在自由度性能与计算复杂度之间存在着一定的利弊权衡问题。

与均匀阵列相比，互质阵列能够采用相同个数的天线阵元取得更多的自由度，因而受到了广泛关注。作为互质采样技术在空间域上的一个典型表现形式，互质阵列提供了一个系统化的稀疏阵列架构方案，并能够突破传统均匀线性阵列自由度受限的瓶颈，实现doa估计方法自由度性能的提升。现有的基于互质阵列的doa估计方法可以通过利用质数的性质将互质阵列推导到虚拟域，并形成等价虚拟均匀线性阵列接收信号以实现doa估计。由于虚拟阵列中包含的虚拟阵元数大于实际的天线阵元数，自由度因此得到了有效的提升。但是，现有基于等价虚拟信号的方法通常需要通过设置预先定义的网格点来设计优化问题，从而进行doa估计。由于实际中信号的来波方向不会完全位于这些预先定义的网格点上，导致了doa估计的精度受限于网格密度。尽管可以通过减小网格间的采样间隔实现估计精度的提升，网格点的密集化将导致方法计算复杂度的指数型增长。此外，现有的doa估计方法的频谱图通常为伪谱，其doa估计方向的波峰响应无法用来表示信号功率；而事实上，信号功率也是描述信号源的一个重要参量，我们希望频谱图能够在估计波达方向的同时反映其功率信息。

技术实现要素：

本发明的目的在于针对上述现有技术存在的不足，提出一种基于信号子空间旋转不变性的互质阵列波达方向估计方法，通过利用互质阵列的特性推导虚拟域等价接收信号，以提升估计方法的自由度；并基于虚拟均匀线性阵列信号子空间的旋转不变性进行doa估计，在进行无网格化波达方向估计的同时实现了相应波达方向上的功率估计，从而在实际应用过程中降低了计算复杂度和硬件复杂度，提高了系统的整体效率和估计精确度。

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

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

(2)利用互质阵列接收k个方向为θ1,θ2,…,θk的远场窄带非相干信号源入射信号，则(2m+n-1)×1维互质阵列接收信号y(t)可建模为：

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

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

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

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

其中，为(2m+n-1)²×k维矩阵，包含k个入射信号源的功率，为噪声功率，i＝vec(i2m+n-1)。这里，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)^*表示共轭操作，表示克罗内克积，i2m+n-1表示(2m+n-1)×(2m+n-1)维单位矩阵。向量z对应的虚拟阵列中各虚拟阵元的位置为

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

(4)构造虚拟阵列协方差矩阵。选取非均匀虚拟阵列中以0为中心连续均匀排列的一段虚拟阵元，形成一个包含l个虚拟阵元的均匀虚拟阵列其相应的虚拟阵元位置为(-lv+1)d到(lv-1)d之间的连续位置，其中，d为入射窄带信号波长的一半，且

相应地，该均匀虚拟阵列的等价信号可通过截取中与该l个虚拟阵元所对应的位置上的元素获得，维度为l×1。虚拟阵列协方差矩阵rv可通过空间平滑技术获得：将向量分割为lv个相互重叠的子向量，每个子向量的维度为lv×1，包含向量中的第i个至第i+lv-1个元素，表示为

则rv可通过取四阶统计量的主平方根获得：

其中，rv的维度为lv×lv，对应于均匀虚拟阵列中位置为0到(lv-1)d的虚拟阵元，称为虚拟均匀线性阵列；且rv在理论上可表示为：

其中，为lv×k维的虚拟均匀线性阵列导引矩阵。∑为对角元素为各信号功率的对角矩阵，表示lv×lv维单位矩阵；

(5)形成虚拟阵列协方差矩阵rv的信号子空间。将获得的虚拟阵列协方差矩阵rv进行特征值分解：

rv＝ωsλsωs^h+ωnλnωn^h，

其中，λs为k×k维对角矩阵，包含将rv的特征值从大到小排列后的前k个特征值，ωs为包含该k个特征值所对应特征向量的lv×k维矩阵；相应地，λn为(lv-k)×(lv-k)维对角矩阵，包含余下的(lv-k)个rv的特征值，ωn为包含该(lv-k)个特征值所对应特征向量的lv×(lv-k)维矩阵；则ωsλsωs^h可被当作rv的信号子空间；

(6)构造具有旋转不变性的两个虚拟子阵信号子空间。可由虚拟均匀线性阵列得到两个包含(lv-1)个虚拟阵元的虚拟子阵x和y，其中，x包含位置为0到(lv-2)d的一段均匀虚拟阵元，y包含位置为d到(lv-1)d的一段均匀虚拟阵元；从理论上来看，x和y的(lv-1)×k维导引矩阵和可由分别移除最后一行和第一行得到，且导引矩阵和间的旋转不变性关系可由一个旋转因子φ表示：

其中，φ为包含这两个虚拟子阵相位差的k×k维对角矩阵；相应地，将实际获得的ωs分别去除最后一行和第一行得到的(lv-1）×k维矩阵ωx和ωy可作为x和y相应的虚拟子阵信号子空间；

(7)计算两个虚拟子阵信号子空间的转换矩阵ψ。两个虚拟子阵信号估计ωx和ωy之间的k×k维转换矩阵ψ可由下式得到：

其中，为ωx的伪逆；

(8)计算波达方向估计结果。信号的波达方向估计值为：

其中，imag(·)表示复数的虚部，ψk为转换矩阵ψ特征值分解后所得的特征值。同时，相应波达方向的信号功率估计值为：

其中，[·]k表示对角矩阵的第k个对角元素，p为包含特征值[ψ1,ψ2,…,ψk]对应特征向量的k×k维矩阵。

进一步地，步骤(1)所述的互质阵列由一对稀疏均匀线性子阵列组合而成，其中第一个子阵列包含2m个天线阵元，阵元间距为nd；第二个子阵列包含n个天线阵元，阵元间距为md。将两个子阵列以首个天线阵元重叠的方式进行组合，得到包含2m+n-1个物理天线阵元的互质阵列架构。

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

(1)本发明充分利用了互质阵列能够增加doa估计的自由度这一优势，将获得的信号模型推导至虚拟域，并通过虚拟阵列协方差矩阵实现doa估计，可分辨的入射信号源个数大于物理天线阵元个数，使得自由度得到提升，并降低了计算复杂度和硬件复杂度；

(2)本发明利用基于虚拟均匀线性阵列信号子空间旋转不变性的方法，无需预先设置网格点来进行波达方向估计，避免了传统方法预先定义网格点所导致的固有估计误差；

(3)本发明在有效估计入射信号doa的同时还能实现相应波达方向的功率估计，所得的空间谱能够同时反映入射信号的波达方向信息和功率信息。

附图说明

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

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

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

图4是本发明所提方法的空间功率谱示意图。

具体实施方式

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

对于doa估计方法在实际系统中的应用，估计精度、自由度和计算复杂度是重要的技术性能指标。现有方法在自由度性能上受限于物理天线阵元的个数，而计算复杂度和估计精度也与预先定义的网格点密度息息相关，且获得的空间谱响应无法有效反映入射信号的功率信息。为了能在增加自由度的同时进行无网格化的doa估计，本发明提供了一种基于信号子空间旋转不变性的互质阵列波达方向估计方法，参照图1，本发明的实现步骤如下：

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

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

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

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

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

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

其中，为(2m+n-1)²×k维虚拟阵列导引矩阵，包含k个入射信号源的功率，为噪声功率，i＝vec(i2m+n-1)。这里，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)^*表示共轭操作，表示克罗内克积，i2m+n-1表示(2m+n-1)×(2m+n-1)维单位矩阵。向量z对应的虚拟阵列中各虚拟阵元的位置为

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

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

相应地，该均匀虚拟阵列的等价信号可通过截取中与该l个虚拟阵元所对应的位置上的元素获得，维度为l×1。接着，虚拟阵列协方差矩阵rv可通过空间平滑技术获得：将向量分割为lv个维度为lv×1的子向量，每个子向量包含向量中的第i个至第i+lv-1个元素，即：

则rv可通过取四阶统计量的主平方根获得：

其中，rv的维度为lv×lv，对应于均匀虚拟阵列中位置为0到(lv-1)d的虚拟阵元，称为虚拟均匀线性阵列，且rv在理论上可以等价地表示为：

其中，为lv×k维的虚拟均匀线性阵列导引矩阵，对应于位置为0到(lv-1)d的虚拟均匀线性阵列。∑为包含各信号功率的对角矩阵，表示lv×lv维单位矩阵。

步骤五：形成虚拟阵列协方差矩阵rv的信号子空间。将获得的虚拟阵列协方差矩阵rv进行特征值分解，以区分信号子空间和噪声子空间：

rv＝ωsλsωs^h+ωnλnωn^h，

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

步骤六：构造具有旋转不变性的两个虚拟子阵信号子空间。将虚拟均匀线性阵列分为x和y两个虚拟子阵，其中，x包含位置为0到(lv-2)d的一段均匀虚拟阵元，y包含位置为d到(lv-1)d的一段均匀虚拟阵元。于是，这两个虚拟子阵的导引矩阵理论值分别为(lv-1)×k维矩阵和其中，和可由步骤四中虚拟均匀线性阵列导引矩阵分别移除最后一行和第一行得到。由于虚拟子阵x和y具有完全相同的均匀且线性的阵列结构，仅由子阵间的单位位移d产生的虚拟子阵移不变性，形成了两个虚拟子阵导引矩阵间的旋转不变性，故和之间的关系可由一个旋转因子φ表示：

其中，φ为包含两个虚拟子阵信号之间的相位差的k×k维对角矩阵；相应地，将步骤五中得到的矩阵ωs分别去除最后一行和第一行后得到的两个(lv-1)×k维矩阵ωx和ωy形成该两个虚拟子阵的信号子空间。

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

其中，为ωx的伪逆。

步骤八：计算波达方向估计结果。根据步骤五中和步骤六中的两个虚拟子阵信号子空间的旋转不变性，则ωx和ωy分别满足以及故转换矩阵ψ满足ψ＝t^-1φt，则将转换矩阵ψ进行特征值分解后得到的包含其特征值[ψ1,ψ2,…,ψk]的对角矩阵和包含对应特征向量的矩阵p分别对应于φ和t^-1。结合转换矩阵ψ的特征值和旋转因子φ对角元素包含的波达方向信息，我们可以得到信号的波达方向估计：

其中，imag(·)表示复数的虚部。同时，将步骤五中代入信号子空间并与步骤四中的对应就可以得到对应的信号功率估计：

其中，[·]k表示对角矩阵的第k个对角元素。

本发明一方面充分利用了互质阵列虚拟域信号处理能够增加doa估计方法自由度的优势，突破了均匀线性阵列自由度受限的瓶颈，实现了在天线阵元个数一定的条件下估计更多个数的入射信号源；另一方面利用了构造的虚拟子阵信号子空间之间的旋转不变性特点，使doa估计能在无网格化的条件下进行，大大降低了计算复杂度，并在进行波达方向估计的同时频谱图能够直观而准确地反映入射信号的功率信息。

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

仿真条件：互质阵列的参数选取为m＝3，n＝5，即架构的互质阵列共包含2m+n-1＝10个天线阵元。假定入射窄带信号个数为15，且入射方向均匀分布于-60^°至60^°这一区间范围内，信噪比为10db，采样快拍数为t＝500，波达方向角的角度域范围为[-90^°,90^°]。

仿真实例：本发明所提出的基于信号子空间旋转不变性的互质阵列波达方向估计方法的空间功率谱如图4所示，其中，虚线表示真实的入射信号波达方向。在本实例的参数设置下，对应于虚拟线性均匀阵列的虚拟阵元的位置为0到17d。可以看出，本发明所提方法能够在利用10个物理阵元的情况下有效地分辨这15个入射信号源，体现了本发明方法在自由度上的优势；此外，空间功率谱的响应值可反映相应波达方向上的信号功率，说明了本发明所提方法能够同时估计各信号的波达方向信息及其对应的功率信息。

综上所述，本发明主要解决了现有技术在doa估计自由度性能与计算复杂度方面存在的不足，一方面充分利用互质阵列的特性在虚拟域进行信号处理以实现自由度的增加；另一方面基于信号子空间旋转不变性的无网格化doa估计有效地降低了计算复杂度，且能够同时估计各相应信号源的功率，在无源定位和目标探测等实际应用中表现出突出的优势。