CN113420411B - 用于无线信号的高分辨率窄带doa估计算法及实现方法 - Google Patents

Abstract

本发明公开了用于无线信号的高分辨率窄带DOA估计算法及实现方法，包括以下步骤：S1、建立二维天线阵列数据模型，在二维均匀矩形阵列模型中，仅在平面角俯仰角的回波方向范围内存在接收信号，对接收数据进行处理时，将回波方向均匀划分成K₁×K₂个网格点，每个网格点对应的平面角和俯仰角分别为和相邻阵元之间的间距为d，本发明SIAA算法针对自定义回波空间进行搜索，与IAA算法相比，进一步提高DOA估计的方向分辨率，而2‑DFFT与GS分解相结合的SIAA算法的快速实现方案则大大降低了运算复杂度，通过计算机仿真实验，与二维IAA算法采用相同网格数进行搜索的情况下，SIAA算法具有更高的角度分辨率，SIAA算法的快速实现大大降低了运算复杂度，增强了SIAA算法的实用性。

Description

用于无线信号的高分辨率窄带DOA估计算法及实现方法

技术领域

本发明涉及DOA估计算法技术领域，具体为用于无线信号的高分辨率窄带DOA估计算法及实现方法。

背景技术

作为信号处理的一个重要分支，阵列信号处理广泛应用在雷达、声呐、地震勘探和射电天文等领域，在阵列信号处理中，DOA估计是其中一个重要研究方向，现有的DOA估计方法大都基于多信号分类算法和子空间旋转不变等算法展开研究，这些算法很好的实现了高分辨率的DOA估计问题，然而当在少量快拍，甚至为单快拍、相干、非相干源或低信噪比等情况下，这些方法无法提供最优的估计结果。

基于WLS方法的参数化迭代自适应算法及其快速实现方案，其基本思路是将回波方向均匀网格化，应用WLS对整个空间进行网格搜索，通过不断更新权重矩阵的方式得到高分辨率DOA估计，然而，在具体情况下，探测方可以通过侦查和判断收发信号双方的相对地理位置等方案确定来波方向的大致范围，IAA算法由于无法利用这一先验信息，仍然需要对整个空间进行网格搜索，在一定程度上增加了运算复杂度

IAA算法是近年来提出的一种高分辨率DOA估计算法，其基本思想是将回波方向均匀网格化，应用最小二乘法对空间网格进行搜索，并通过迭代更新加权矩阵获得高分辨率频谱估计，然而在回波方向可确定情形下，无需对整个空间进行搜索，为进一步提高估计的分辨率，针对特定回波方向展开搜索的窄带IAA算法，同时，为了提高估计效率，将2-DFFT与GS分解方法相结合，提出了SIAA算法的快速实现方案。

发明内容：

本发明的目的就在于为了解决上述问题而提供用于无线信号的高分辨率窄带DOA估计算法及实现方法，解决了背景技术中提到的问题。

为了解决上述问题，本发明提供了一种技术方案：

用于无线信号的高分辨率窄带DOA估计算法，包括以下步骤：

S1、建立二维天线阵列数据模型，在二维均匀矩形阵列模型中，仅在平面角俯仰角/>的回波方向范围内存在接收信号，对接收数据进行处理时，将回波方向均匀划分成K₁×K₂个网格点，每个网格点对应的平面角和俯仰角分别为和/>相邻阵元之间的间距为d，二维天线阵列数据模型表示为：

S2、通过二维IAA算法对估计值进行初始化，然后通过矩阵R_MN和式(3)不断更新估计值并最终获得稳定的频谱估计，当接收信号的回波方向大致范围已知时，开始实施SIAA算法；

S3、通过GS分解和2-DFFT的方法降低SIAA算法运算复杂度；

S4、通过计算机仿真验证SIAA算法的高分辨率性能以及快速SIAA算法的有效性和运算效率；

S5、对SIAA快速算法与SIAA算法进行DOA估计所用时间进行对比；

所述S1中y_MN表示MN×1维的接收数据向量，表示信号振幅向量，表示MN×K₁K₂阶的阵列流形矩阵，q_MN表示独立的服从高斯分布的加性噪声，/>表示信源的导向矢量，且，/>表示Kronecker积，/> (·)^T表示转置。

所述S1中式(1)中的参数满足 的估计值，即，通过WLS方法得到：

且，干扰和噪声的协方差矩阵定义为(·)^H表示共轭转置，/>最小化公式(2)可以得到：

所述s_mn的估计值，即，可以表示为：

且，X_MN表示窄带模型下的加权矩阵，根据加权WLS方法中加权矩阵的定义，X_MN表示为：

表示不存在回波的空间范围内噪声的平均功率，/>通过

进行求解，

由于X_MN中的元素是关于/>的隐函数，需要采用迭代的方法，第t次迭代的估计值/>表示为：

且/>可以通过式(5)进行更新；所述S3中2-DFFT的模型，/>

为实现SIAA的快速算法，需要满足H₁≥3M,H₂≥3N，假设相邻阵元之间的间距为一倍波长，则式(7)可重新表示为：

且， 将回波方向均匀划分成/>个网格点，且，/>式(1)中的阵列流形矩阵重新定义为/>提取/>波峰所在位置的/>和/>的值，此时接收到的信号的平面角/>俯仰角

所述S3中2-DFFT方法构建加权矩阵由(式10)可知，/>和/>的计算均需对/>矩阵求逆，根据(式5)中的定义，/>矩阵具有TBT结构，因此可以将/>矩阵用分块矩阵表示，即：

且，

根据(式11)和(式12)可知，表示为：

且，m＝-M+1,...,-1,0,1,...,M-1，n＝-N+1,...,-1,0,1,...,N-1；

由(式13)可知，元素构成的矩阵表示为：

且，表示对/>矩阵做长度为H₁×H₂的2-DFFT，表示计算过程中可忽略具体内容的矩阵，

0_mn表示由全0元素组成的m×n阶矩阵，而则表示由估计值的平方构成的阶的矩阵，将(式14)和(式15)代入式(11)，构建矩阵/>

作为优选，所述可分解成/>和/>两部分分别计算，将/>矩阵表示成分块矩阵的形式：

且，J_N表示副对角线上元素值为1，其余元素为0的N×N阶矩阵，引入：

根据GS分解[ii]可以将的逆矩阵表示为：

且，并且满足：

将上述矩阵分解成向量的形式， 和/>通过FFT/IFFT方法实现向量z_MN快速计算，然后对/>进行计算，对z_MN向量做长度为H₁、H₂的2-DFFT得到此时：/>

且，表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。

作为优选，所述的快速实现：

对于二维天线阵列模型，表示为：

根据(式22)，矩阵可以通过2-DFFT的方式计算得到，定义矩阵C：

且，0表示为构成H₁×H₂阶的C矩阵所需的由全0元素组成的矩阵，矩阵C₁，C₂，C₃和C₄由多项式系数表示：

将向量和/>表示成分块矩阵的形式，令MN×1阶的向量x_MN表示为x_M ^ij＝[x^1,ij,x^2,ij,...,x^M,ij]^T，矩阵C₁，C₂，C₃和C₄中元素计算过程如下：

且，

此时通过(式27)和(式28)得到的C矩阵中的元素值再根据(式25)和(式26)构建C矩阵，接下来根据(式23)可得此时：

且，/>表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。

本发明的有益效果是：SIAA算法针对自定义回波空间进行搜索，与IAA算法相比，进一步提高DOA估计的方向分辨率，而2-DFFT与GS分解相结合的SIAA算法的快速实现方案则大大降低了运算复杂度，通过计算机仿真实验，与二维IAA算法采用相同网格数进行搜索的情况下，SIAA算法具有更高的角度分辨率，SIAA算法的快速实现大大降低了运算复杂度，增强了SIAA算法的实用性。

附图说明：

为了易于说明，本发明由下述的具体实施及附图作以详细描述。

图1是本发明二维接收天线阵列示意图；

图2是本发明分辨率验证SIAA算法DOA估计频谱图；

图3是本发明分辨率验证IAA算法DOA估计频谱图；

图4是本发明有效性验证SIAA算法DOA估计频谱图；

图5是本发明有效性验证SIAA快速算法DOA估计频谱图；

图6是本发明SIAA算法时间展示图；

图7是本发明SIAA快速算法时间展示图。

具体实施方式：

如图1-7所示本具体实施方式采用以下技术方案：

实施例：

用于无线信号的高分辨率窄带DOA估计算法，包括以下步骤：

S1、建立二维天线阵列数据模型，在二维均匀矩形阵列模型中，仅在平面角俯仰角/>的回波方向范围内存在接收信号，对接收数据进行处理时，将回波方向均匀划分成K₁×K₂个网格点，每个网格点对应的平面角和俯仰角分别为和/>相邻阵元之间的间距为d，二维天线阵列数据模型表示为：

S2、通过二维IAA算法对估计值进行初始化，然后通过矩阵R_MN和式(3)不断更新估计值并最终获得稳定的频谱估计，当接收信号的回波方向大致范围已知时，开始实施SIAA算法；

S3、通过GS分解和2-DFFT的方法降低SIAA算法运算复杂度；

S4、通过计算机仿真验证SIAA算法的高分辨率性能以及快速SIAA算法的有效性和运算效率；

S5、对SIAA快速算法与SIAA算法进行DOA估计所用时间进行对比；

所述S1中y_MN表示MN×1维的接收数据向量，表示信号振幅向量，表示MN×K₁K₂阶的阵列流形矩阵，q_MN表示独立的服从高斯分布的加性噪声，/>表示信源的导向矢量，且，/>表示Kronecker积，/> (·)^T表示转置；

所述S1中式(1)中的参数满足 的估计值，即，通过WLS方法得到：

且，干扰和噪声的协方差矩阵定义为(·)^H表示共轭转置，/>最小化公式(2)可以得到：

所述s_mn的估计值，即，可以表示为：

且，X_MN表示窄带模型下的加权矩阵，根据加权WLS方法中加权矩阵的定义，X_MN表示为：

表示不存在回波的空间范围内噪声的平均功率，/>通过

进行求解，

由于X_MN中的元素是关于/>的隐函数，需要采用迭代的方法，第t次迭代的估计值/>表示为：

且/>可以通过式(5)进行更新；

所述S3中2-DFFT的模型，

为实现SIAA的快速算法，需要满足H₁≥3M,H₂≥3N，假设相邻阵元之间的间距为一倍波长，则式(7)可重新表示为：

且， 将回波方向均匀划分成/>个网格点，且，/>式(1)中的阵列流形矩阵/>重新定义为/>提取/>波峰所在位置的/>和/>的值，此时接收到的信号的平面角/>俯仰角/>

所述S3中2-DFFT方法构建加权矩阵由(式10)可知，/>和/>的计算均需对/>矩阵求逆，根据(式5)中的定义，/>矩阵具有TBT结构，因此可以将/>矩阵用分块矩阵表示，即：

且，

根据(式11)和(式12)可知，表示为：

且，m＝-M+1,...,-1,0,1,...,M-1，n＝-N+1,...,-1,0,1,...,N-1；

由(式13)可知，元素构成的矩阵表示为：

且，表示对/>矩阵做长度为H₁×H₂的2-DFFT，表示计算过程中可忽略具体内容的矩阵，

0_mn表示由全0元素组成的m×n阶矩阵，而则表示由估计值的平方构成的阶的矩阵，将(式14)和(式15)代入式(11)，构建矩阵/>

其中，所述可分解成/>和/>两部分分别计算，将/>矩阵表示成分块矩阵的形式：

且，J_N表示副对角线上元素值为1，其余元素为0的N×N阶矩阵，引入：

根据GS分解可以将的逆矩阵表示为：

且，并且满足：

将上述矩阵分解成向量的形式， 和/>通过FFT/IFFT方法实现向量z_MN快速计算，然后对/>进行计算，对z_MN向量做长度为H₁、H₂的2-DFFT得到此时：/>

且，表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。

其中，所述的快速实现：

对于二维天线阵列模型，表示为：

根据(式22)，矩阵可以通过2-DFFT的方式计算得到，定义矩阵C：

且，0表示为构成H₁×H₂阶的C矩阵所需的由全0元素组成的矩阵，矩阵C₁，C₂，C₃和C₄由多项式系数表示：

将向量和/>表示成分块矩阵的形式，令MN×1阶的向量x_MN表示为x_M ^ij＝[x^1,ij,x^2,ij,...,x^M,ij]^T，矩阵C₁，C₂，C₃和C₄中元素计算过程如下：

且，

此时通过(式27)和(式28)得到的C矩阵中的元素值再根据(式25)和(式26)构建C矩阵，接下来根据(式23)可得此时：

且，/>表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。

具体的，首先，将二维矩形阵列作为本文的数据模型，在平面角俯仰角/>的回波方向范围内接收信号，将回波方向均匀划分成K₁×K₂个网格点，每个网格点对应的平面角和俯仰角分别为/>和/>相邻阵元之间的间距为d，通过二维IAA算法实现过程需要先对估计值/>进行初始化，然后通过矩阵R_MN和式(3)不断更新估计值/>并最终获得稳定的频谱估计，当接收信号的回波方向大致范围已知时，原始IAA算法无法自适应地调整搜索空间，此时采用SIAA算法，a、通过式(1)中的定义构建阵列流形矩阵/>b、对估计值/>进行初始化；c、根据定义构建矩阵/>d、通过式(5)和式(6)计算矩阵X_MN；e、通过式(7)计算估计值/>循环步骤c-e，直到/>t表示迭代次数，同时通过GS分解和2-DFFT的方法降低运算复杂度，从而实现SIAA的快速算法，根据加权矩阵/>的特殊结构，利用2-DFFT对矩阵/>进行快速构建，然后，通过2-DFFT与GS分解相结合的方法分别对矩阵/>和/>进行处理，最后，将上述方案融合到每次迭代的过程中以实现SIAA的快速算法，二维SIAA快速算法迭代过程为，a1、根据定义构建阵列流形矩阵/>b1、通过式(15)对/>矩阵进行初始化；c1对初始化后的矩阵/>做2-D FFT，构造矩阵/>d1、通过迭代方式计算矩阵W和Q；e1、通过矩阵W和Q计算式(19)中的矩阵T₁，T₂，S₁和S₂；f1、通过式(21)和(28)计算矩阵/>和/>g1、通过式(11)得到/>循环步骤c1-g1直到/>t表示迭代次数，将SIAA快速算法的计算复杂度同SIAA算法相对比，设定H₁＝3M，H₂＝3N，SIAA算法中的协方差矩阵/>由矩阵和/>两部分构成，通过重构的方法计算/>矩阵，计算量为O(H₁H₂log₂(H₁H₂))，直接计算需要的计算量为/>和/>矩阵所需的计算量为O(M²N²+MNlog₂(MN))+O(H₁H₂log₂(H₁H₂))，直接计算/>和矩阵所需计算量为/>由上述分析可知，快速实现方案的提出从理论上降低了SIAA算法的运算复杂度，对所提SIAA算法的分辨率进行验证，设定坐标轴上相邻阵元之间的间距为接收信号波长的3倍，SIAA算法同IAA算法均采用8100个网格进行搜索，已知在方位角45°至55°，俯仰角55至65°的空间范围内存在多个来自不同方向的信号。通过SIAA算法与IAA算法处理后的频谱如图2、3所示，图3中，通过IAA算法处理后得到的频谱在来波密集的区域仅观察到一个波峰，而图2中通过SIAA算法处理后得到的频谱在相同区域观察到2～3个波峰，能够将来波密集区域的信号区分开，与IAA算法相比，SIAA算法具有更高的方向分辨率，对快速SIAA算法的有效性进行验证，设定坐标轴上相邻阵元之间的间距等于接收信号的波长，SIAA算法与快速SIAA算法具有相同的窄带搜索范围，令/>搜索范围0.45≤u≤0.6，0.6≤v≤0.8，均采用3600个网格数进行搜索，通过SIAA算法与快速SIAA算法处理后得到的频谱如图4、5所示，SIAA快速算法与SIAA算法处理后得到的频谱在相同方向均能观察到四个波峰，SIAA快速算法能够实现SIAA算法的频谱估计，具有同样的高分辨率，对SIAA快速算法与SIAA算法进行DOA估计所用时间进行对比，设定坐标轴上相邻阵元之间的间距等于接收信号的波长，SIAA算法与SIAA快速算法具有相同的窄带搜索范围，令/>搜索范围0.45≤u≤0.6，0.6≤v≤0.8，随着天线阵元数M和N的变化，SIAA算法与SIAA快速算法进行DOA估计所需时间对比如图6、7所示，时间取为10log(·)，与SIAA算法相比，SIAA快速算法随着天线阵元数目M和N的增加，进行DOA估计所用时间大幅减少，SIAA算法的快速实现大大减少了迭代运算过程的计算量，提高了运算效率。

在本发明的描述中，需要理解的是，术语“同轴”、“底部”、“一端”、“顶部”、“中部”、“另一端”、“上”、“一侧”、“顶部”、“内”、“前部”、“中央”、“两端”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (4)

1.用于无线信号的高分辨率窄带DOA估计算法，其特征在于，包括以下步骤：

S1、建立二维天线阵列数据模型，在二维均匀矩形阵列模型中，仅在平面角俯仰角/>的回波方向范围内存在接收信号，对接收数据进行处理时，将回波方向均匀划分成K₁×K₂个网格点，每个网格点对应的平面角和俯仰角分别为和/>相邻阵元之间的间距为d，二维天线阵列数据模型表示为：

S2、通过二维IAA算法对估计值进行初始化，然后通过矩阵R_MN和式(3)不断更新估计值并最终获得稳定的频谱估计，当接收信号的回波方向范围已知时，开始实施SIAA算法；

S3、通过GS分解和2-DFFT的方法降低SIAA算法运算复杂度；

S4、通过计算机仿真验证SIAA算法的高分辨率性能以及快速SIAA算法的有效性和运算效率；

S5、对SIAA快速算法与SIAA算法进行DOA估计所用时间进行对比；

所述S1中y_MN表示MN×1维的接收数据向量，表示信号振幅向量，表示MN×K₁K₂阶的阵列流形矩阵，q_MN表示独立的服从高斯分布的加性噪声，/>表示信源的导向矢量，且，/>表示Kronecker积，/> (·)^T表示转置；

所述S1中式(1)中的参数满足 的估计值，即，通过WLS方法得到：

且，干扰和噪声的协方差矩阵定义为(·)^H表示共轭转置，最小化公式(2)可以得到：

所述s_mn的估计值，即，可以表示为：

且，X_MN表示窄带模型下的加权矩阵，根据加权WLS方法中加权矩阵的定义，X_MN表示为：

表示不存在回波的空间范围内噪声的平均功率，/>通过

进行求解，

由于X_MN中的元素是关于/>的隐函数，需要采用迭代的方法，第t次迭代的估计值/>表示为：

且/>可以通过式(5)进行更新；

所述S3中2-DFFT的模型，

为实现SIAA的快速算法，需要满足H₁≥3M,H₂≥3N，假设相邻阵元之间的间距为一倍波长，则式(7)可重新表示为：

且， 将回波方向均匀划分成/>个网格点，且，/>式(1)中的阵列流形矩阵重新定义为/>提取/>波峰所在位置的/>和/>的值，此时接收到的信号的平面角/>俯仰角

2.根据权利要求1所述的用于无线信号的高分辨率窄带DOA估计算法，其特征在于，所述S3中2-DFFT方法构建加权矩阵由(式10)可知，/>和/>的计算均需对矩阵求逆，根据(式5)中的定义，/>矩阵具有TBT结构，因此可以将/>矩阵用分块矩阵表示，即：

且，

根据(式11)和(式12)可知，表示为：

且，m＝-M+1,...,-1,0,1,...,M-1，n＝-N+1,...,-1,0,1,...,N-1；

由(式13)可知，元素构成的矩阵表示为：

且，表示对/>矩阵做长度为H₁×H₂的2-DFFT，表示计算过程中可忽略具体内容的矩阵，

0_mn表示由全0元素组成的m×n阶矩阵，而则表示由估计值的平方构成的/>阶的矩阵，将(式14)和(式15)代入式(11)，构建矩阵/>

3.根据权利要求2所述的用于无线信号的高分辨率窄带DOA估计算法，其特征在于，所述可分解成/>和/>两部分分别计算，将/>矩阵表示成分块矩阵的形式：

且，J_N表示副对角线上元素值为1，其余元素为0的N×N阶矩阵，引入：

根据GS分解[i]可以将的逆矩阵表示为：

且，并且满足：

将上述矩阵分解成向量的形式， 和/>通过FFT/IFFT方法实现向量z_MN快速计算，然后对/>进行计算，对z_MN向量做长度为H₁、H₂的2-DFFT得到/>此时：/>

且，表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。

4.根据权利要求3所述的用于无线信号的高分辨率窄带DOA估计算法，其特征在于，所述的快速实现：

对于二维天线阵列模型，表示为：

根据(式22)，矩阵可以通过2-DFFT的方式计算得到，定义矩阵C：

且，0表示为构成H₁×H₂阶的C矩阵所需的由全0元素组成的矩阵，矩阵C₁，C₂，C₃和C₄由多项式系数/>表示：

将向量和/>表示成分块矩阵的形式，令MN×1阶的向量x_MN表示为x_M ^ij＝[x^1,ij,x^2,ij,...,x^M,ij]^T，矩阵C₁，C₂，C₃和C₄中元素计算过程如下：

且，

此时通过(式27)和(式28)得到的C矩阵中的元素值再根据(式25)和(式26)构建C矩阵，接下来根据(式23)可得此时：

且，/>表示/>矩阵中h₁₁至h₁₂行、h₂₁至h₂₂列的元素构成的/>阶矩阵。