CN116299150B - 一种均匀面阵中降维传播算子的二维doa估计方法 - Google Patents

Abstract

本发明公开了一种均匀面阵中降维传播算子的二维DOA估计方法，包括如下步骤：利用整个均匀面阵的接收信号构造协方差矩阵；对构造的协方差矩阵进行分块处理，得到各个分块矩阵，通过传播算子矩阵的旋转不变性进行初始估计；接着，通过降维处理，将二维优化问题转变为一维优化问题；然后，根据最小二乘思想直接求解目标信号的方位角信息作为初估计。最后，对初估计进行转换得到均匀面阵中的二维DOA估计值。本发明不需对接收信号的协方差矩阵进行特征值分解，只需一维局部谱峰搜索即可获得更加精确的估计值，以及二维估计参数的自动配对，所提方法的性能与传统均匀面阵二维PM算法相近。此外，该方法在保证估计性能的同时显著降低了计算复杂度。

Description

一种均匀面阵中降维传播算子的二维DOA估计方法

技术领域

本发明属于阵列信号处理DOA估计技术领域，尤其涉及一种均匀面阵中降维传播算子的二维DOA估计方法。

背景技术

阵列信号处理时信号处理的一个重要分支，其中的二维DOA估计是阵列信号处理的一个重要内容，也是阵列信号方向重要的研究领域。由于二维DOA估计可以更加准确的定位到某一确定空间中的多个感兴趣的空间信号的方向角与俯仰角，所获取的信息量更大，使得近年来关于二维DOA估计的相关内容越来越受到人们的关注。

现有的关于DOA估计的二维PM算法的角度估计精度很高，许多学者将二维PM算法运用于双L型阵列、立体十字型互耦阵列、双平行线阵等DOA估计中，得到了较好的DOA估计性能。但是传统二维PM算法进行DOA估计时需要进行全局谱峰搜索，这也带来了高维搜索问题，计算复杂度大大增加。为了解决上述问题，本发明提出了一种均匀面阵中降维传播算子的二维DOA估计方法。

发明内容

发明目的：为了解决现有技术存在的问题，本发明提供了一种均匀面阵中降维传播算子的二维DOA估计方法，在保证估计性能的同时显著降低了计算复杂度，易于实时处理。

技术方案：本发明所述的一种均匀面阵中降维传播算子的二维DOA估计方法，包括以下步骤：

一种均匀面阵中降维传播算子的二维DOA估计方法，包括以下步骤：

1)构造均匀面阵的信号模型，得到接收信号模型x(t)，构造所述接收信号模型x(t)的协方差矩阵

2)对构造的协方差矩阵进行分块处理，通过传播算子矩阵的旋转不变性进行初始估计，得到用于谱峰搜索的初始估计值/>

3)通过降维处理，将二维优化问题转变为一维优化问题，进行一维局部谱峰搜索获得俯仰角u的估计值

4)根据最小二乘思想获取方位角v的估计值

5)利用估计值和/>得到均匀面阵中的二维DOA精确估计值/>

优选的，步骤1)中：

所述均匀面阵共有M×N个均匀分布的阵元，M与N分别为x轴与y轴上的阵元个数，x轴与y轴上的相邻阵元间距均为d，且d≤λ/2，λ为波长，其中将每一列的阵元集合称作一个子阵，则有子阵1,2,…,n,…,N，共N个子阵；

假设某一空间有K个信号源以二维波达方向θ_k,入射到所述均匀面阵上，其中θ_k,/>分别表示第k(k＝1,2,…,K)个信号源的仰角与方位角，定义接收信号x(t)为：

其中：

其中：为Kronecker乘积，/>为方向矩阵A_y和A_x的Khatri-Rao乘积；

为信号源矢量，s_K(t)为第K个信号源矢量，(·)^T为转置；n(t)＝[n₁(t)^T,n₂(t)^T,…,n_N(t)^T]^T为所述均匀面阵的加性高斯白噪声，n_n(t)为子阵n的加性高斯白噪声；/>为x轴上的信号源方向向量，/>为y轴上的信号源方向向量；

构造接收信号x(t)的协方差矩阵为：

其中，L为快拍数，()^H表示共轭转置。

优选的，步骤2)中：

定义参量u和v分别为将x轴与y轴上的信号源方向向量和/>转化为：

a_x(v)＝[1,e^j2dv/λ,…,e^{-j2πd(N-1)v/λ}]^T

a_y(u)＝[1,e^-j2πdu/λ,…,e^{-jπ2d(M1)u/λ}]^T

构造矩阵P_c为：

其中，P为传播算子，I_K为K阶单位矩阵，()^H为共轭转置；

将矩阵P_c划分为N块：

P_c＝P_c1,P_c2,…,P_cN ^T

其中，矩阵n＝1,…,N；

分别构造矩阵P₁和P₂：

P₁＝P_c1,P_c2,…,P_c(N-1) ^T

P₂＝P_c2,P_c3,…,P_cN ^T

得到：

P₂＝P₁T^-1ΦT

其中，为对角矩阵，u_k为第k个信号源的俯仰角，T为K×K阶的满秩矩阵，然后将P₁ ⁺·P₂进行特征值分解，记得到的第k个特征值为p_k，(·)⁺为求广义逆，·表示相乘，从而得出u_k的初始估计值/>为：

其中k＝1,2,,K。

优选的，步骤3)中：

定义矩阵Q：

其中，P为传播算子，I_MN-K为MN-K阶单位矩阵；

定义参数V(u,v)：

将V(u,v)变形为：

为了便于后续的计算，令参数I_M为M阶单位矩阵；

用e₁ ^Ha_x(v)＝1来消除a_x(v)＝0_M的平凡解，其中矩阵考虑到e₁ ^Ha_x(v)＝1，则将变形之后的V(u,v)重构为：

然后构造代价函数L(u,v)：

L(u,v)＝a_x(v)^HC(u)a_x(v)-ω(e₁ ^Ha_x(v)-1)

其中，ω为一个常数，则有：

则有a_x(v)＝μC(u)^-1e₁，由于e₁ ^Ha_x(v)＝1，则μ＝1/e₁ ^HC^-1(u)e₁为常量，所以a_x(v)表示为：

结合得到俯仰角u的估计值/>

将估计值用另一种表示方式为：

优选的，步骤4)中：

在区间内对参量u进行一维局部搜索，得到C^-1(u)中第(1,1)元素中最大的K个值，其中Δu→0；这最大的K个值的峰值(u₁,u₂,…,u_K)便与sinθ_ksinφ_k,(k＝1,2,…,K)的估计值相对应；

然后求出K个导向矢量

由于得出：

g_k＝-angle(a_x(v_k))＝0,2πdv_k/λ,…,2π(M-1)dv_k/λ^T＝v_kq

为了便于进一步推导，定义参数q＝0,2πd/λ,…,2π(M-1)d/λ^T；归一化再通过最小二乘法求v_K，最小二乘法则为：

其中：E＝[1_M,q]，

c_k1是v_k的估计值；则c_k的解为：

得到v_k的估计值

优选的，步骤5)中：

将和/>代入上式即得到均匀面阵中的二维DOA精确估计值/>

有益效果：与现有技术相比，本发明的有益效果：本发明所提算法不需对接收信号的协方差矩阵进行特征值分解，而是首先利用传播算子矩阵的旋转不变性进行初始估计，然后在通过降维处理，将二维优化问题转变为一维优化问题，仅需一维局部谱峰搜索即可获得更加精确的估计值，同时还可以实现参数估计角度的自动配对，角度估计性能与2D-ESPRIT算法估计性能相近；同时，与均匀面阵中传统的二维PM算法相比，所提算法极大的降低了计算复杂度；本发明可以在保证估计性能的情况下显著降低计算复杂度。

附图说明

图1为本发明的流程图；

图2为均匀面阵分布图；

图3为本发明所述方法与2D-PM算法计算复杂度随信源数K变化示意图；

图4为本发明所述方法与2D-PM算法计算复杂度随快拍数J变化示意图；

图5为本发明所述方法与2D-PM算法计算复杂度随阵元数M变化示意图；

图6为本发明所述方法与2D-PM算法计算复杂度随阵元数N变化示意图；

图7为本发明所述方法的角度估计散点图；

图8为本发明所述方法与2D-ESPRIT算法、2D-MUSIC算法的方位角角度估计性能对比示意图；

图9为本发明所述方法与2D-ESPRIT算法、2D-MUSIC算法的俯仰角角度估计性能对比示意图。

具体实施方式

下面结合附图对本发明做进一步详细说明。

本发明提供一种均匀面阵中降维传播算子的二维DOA估计方法，如图1所示，具体包括以下步骤：

步骤1：构造如图2所示的均匀面阵的信号模型，得到接收信号模型x(t)，构造所述接收信号模型x(t)的协方差矩阵

均匀面阵共有M×N个均匀分布的阵元(M与N分别为x轴与y轴上的阵元个数)，且相邻阵元间距为d，(d≤λ/2，λ为波长)。假设某一空间有K个信号源以二维波达方向(θ_k,φ_k)入射到均匀面阵上，其中θ_k,φ_k分别表示第k个信号源的仰角与方位角。x轴与y轴上的信号源方向向量分别为

均匀面阵中子阵1的接收信号可以表示为x₁(t)＝A_xs(t)+n₁(t)；其中，为x轴上子阵1对应的方向矩阵，n₁(t)是子阵1的加性高斯白噪声，/>是信号源矢量。

根据均匀面阵的特性，子阵1的方向矩阵为A_x，而子阵2的方向矩阵就要考虑沿y轴的偏移情况，每个阵元相对于参考阵元(子阵1)的波程差就等于子阵1的阵元的波程差加上2πdsinφsin/λ，所以可得子阵n的接收信号为其中，为由矩阵的n行构造的一个对角矩阵，n_n(t)是子阵n的加性高斯白噪声；

因此，得到整个均匀面阵的接收信号为：

其中，x_i(t),i＝1,2,,N表示第i个子阵的接收信号，表示A_y和A_x的Khatri-Rao乘积，/>为y轴上子阵1对应的方向矩阵，n(t)＝[n₁(t)^T,n₂(t)^T,…,n_N(t)^T]^T为均匀面阵的加性高斯白噪声，(·)^T表示转置；

根据Khatri-Rao乘积的定义，接收信号x(t)可以转化为：

其中，表示Kronecker乘积，/>

对矩阵A_x分块得

在假设阵列无空间模糊(A_x列为满秩)的情况下，是非奇异矩阵，而可以由A_x1的线性变换得到。所以，该阵列方向矩阵A可以写成

其中，A_x2＝P^HA_x1，P是传播算子，是单位矩阵。定义矩阵Q为

可以得到Q^HA_x＝0_(MN-K)×K，定义协方差矩阵为其中，L为快拍数。

步骤2：对构造的协方差矩阵进行分块处理，得到各个分块矩阵，通过传播算子矩阵的旋转不变性进行初始估计，得到用于谱峰搜索的初始估计值/>

其次，定义参量u和v分别为则/>和/>可以转化为：

a_x(v)＝[1,e^-j2πdv/λ,…,e^{-j2πd(N1)v/λ}]^T

a_y(u)＝[1,e^-j2πdu/λ,…,e^{-j2πd(M1)u/λ}]^T

然后构造矩阵P_c为

其中，P为传播算子，I_K为一个K阶单位矩阵。

将矩阵P_c划分为N块P_c＝P_c1,P_c2,…,P_cN ^T，其中，n＝1,…,N；然后分别构造矩阵P₁和P₂

P₁＝P_c1,P_c2,…,P_c(N-1) ^T

P₂＝P_c2,P_c3,…,P_cN ^T

可以得到

P₂＝P₁T^-1ΦT

其中，为由K个互不相关的信号源到达阵元的波程差构成的对角矩阵，式中u_k表示第k个(k＝1,2,...,K)互不相关信号源的俯仰角构成的方向信息；T为一个K×K阶的满秩矩阵；然后将P₁+P₂进行特征值分解得到的第k个特征值记为p_k，()⁺表示求广义逆，从而得出u_k的初始估计值/>为：

其中k＝1,2,,K。

步骤3：通过降维处理，将二维优化问题转变为一维优化问题，实现参数估计角度的自动配对，仅需一维局部谱峰搜索可获得俯仰角u的估计表达式

接着定义

上式还可变形为

其中，为了便于后续的计算，令I_M为一个M阶单位矩阵。上式为二次优化问题，用e₁ ^Ha_x(v)＝1来消除a_x(v)＝0_M的平凡解，其中考虑到e₁ ^Ha_x(v)＝1，则上述二次优化问题可以重构为

然后构造如下的代价函数

L(u,v)＝a_x(v)^HC(u)a_x(v)-λ(e₁ ^Ha_x(v)-1)

其中，λ为一个常数，则有

根据上式，可以得到a_x(v)＝μC(u)^-1e₁，由于e₁ ^Ha_x(v)＝1，则μ＝1/e₁ ^HC^-1(u)e₁为一个常量。所以a_x(v)可以表示为

结合和上式，可以得到u的估计表达式为

此外，还可以表示为

步骤4：根据最小二乘思想直接求解目标信号的方位角信息作为DOA的初估计值v_K。

通过对u在区间内进行一维局部搜索，从而得到C^-1(u)中第(1,1)元素中最大的K个值，其中Δu→0；这最大的K个值的峰值(u₁,u₂,…,u_K)便与sinθ_ksinφ_k,(k＝1,2,…,K)的估计相对应；

然后可以求出K个导向矢量

由可以得出

g_k＝-angle(a_x(v_k))＝0,2πdv_k/λ,…,2π(M-1)dv_k/λ^T＝v_kq

其中，为了便于进一步推导，令q＝0,2πd/λ,…,2π(M-1)d/λ^T；归一化再通过最小二乘法求v_K。其中最小二乘法则为/> E＝[1_M,q]，c_k1是v_k的估计值。则c_k的解可以表示为

因为前面已经表示出是一个2×1阶矩阵，根据最小二乘法则可以求出/>的表达式即为/>也是一个2×1阶矩阵，所以该矩阵的第一行第一列元素值为c_k0的结果，

根据上式即可得出v_k的估计值。

步骤5：利用初估计值和/>通过/>进行角度信息转换，得到均匀面阵中的二维DOA精确估计值。

其中分别为

将俯仰角初估计值和方位角初始估计值/>代入上式，即可得到均匀面阵中二维DOA估计的俯仰角和方位角信息。

本发明方法得到的计算的复杂度为O{JM²N²}，初始估计的复杂度为O{KM²N²+K²MN+2K³+2K²M(N-1)}，局部谱峰搜索需要的复杂度为O{n₁K[(M²N+M²)(MN-K)+M²]}，则降维PM算法的二维DOA估计所需的总复杂度为O{KM²N²+K²MN+2K³+2K²M(N-1)+n₁K[(M²N+M²)(MN-K)+M²]+JM²N²}，其中n₁是局部谱峰搜索的循环次数；而2D-PM算法的总复杂度为O{KM²N²+2K³+n₂[(MN+1)(MN-K)]+JM²N²}，其中n₂为全局谱峰搜索的循环次数，且n₁＜＜n₂；所以本文算法的计算复杂度远小于2D-PM算法。

图3为本发明所述方法与传统2D-PM的DOA估计方法的计算复杂度(复数乘法次数)随信源数K数目变化的示意图。仿真条件为：信噪比SNR＝50dB，阵元数M＝8，N＝8、快拍数J＝100。

图4为本发明所述方法与传统2D-PM的DOA估计方法的计算复杂度(复数乘法次数)随快拍数J数目变化的示意图。仿真条件为：信噪比SNR＝50dB，阵元数M＝8，N＝8、信源数K＝3。

图5为本发明所述方法与传统2D-PM的DOA估计方法的计算复杂度(复数乘法次数)随阵元数M变化的示意图。仿真条件为：信噪比SNR＝50dB，快拍数J＝100、信源数K＝3、阵元数N＝8。

图6为本发明所述方法与传统2D-PM的DOA估计方法的计算复杂度(复数乘法次数)随阵元数N变化的示意图。仿真条件为：信噪比SNR＝50dB，快拍数J＝100、信源数K＝3、阵元数M＝8。

从图3至图6中可以看出，本发明所述的降维PM方法无需谱峰搜索，计算复杂度远小于2D-PM算法。

本发明性能估计标准为均方根误差(root mean square error,RMSE)，定义RMSE、方位角的RMSE和俯仰角的RMSE分别为

其中，分别为第l次蒙特卡洛仿真时θ_k、φ_k的估计值。假设某一空间中有K＝3个信号源，方位角和俯仰角分别为/> 阵元间距为d＝λ/2；M与N分别代表x轴与y轴上的阵元数量，J代表快拍数。

图7为本发明所述方法的角度估计散点图。仿真条件为：有3个信源，它们的方位角与俯仰角分别为均匀面阵M＝8，N＝10，快拍数J为100，信噪比SNR为20dB，仿真1000次。从图7可以看出，本发明可以较准确的对信源的俯仰角进行估计。而且，在大信噪比情况下，算法的角度估计性能更好。

图8为本发明所述方法与2D-ESPRIT算法、2D-MUSIC算法的方位角角度估计性能对比。仿真条件为：有3个信源，它们的方位角与俯仰角分别为 均匀面阵M＝8，N＝8，快拍数J为100，信噪比SNR为50dB，仿真1000次。

图9为本发明所述方法与2D-ESPRIT算法、2D-MUSIC算法的俯仰角角度估计性能对比。仿真条件为：有3个信源，它们的方位角与俯仰角分别为 均匀面阵M＝8，N＝8，快拍数J为100，信噪比SNR为50dB，仿真1000次。

从图8和图9可以看出，本文所提降维PM(RD-PM)方法的角度估计性能与2D-EPRIT算法性能相近，均优于2D-MUSIC算法的估计性能。

综上所述，从仿真效果图的分析可知，本发明提出的一种均匀面阵中降维传播算子的二维DOA估计方法不需要进行二维全局谱峰搜索，而是通过一维局部谱峰搜索得到方向角与俯仰角的估计值，成功减少了二维全局谱峰搜索带来的巨大计算量，显著降低了计算复杂度，而且还可以实现角度估计参数的自动配对。同时，通过仿真验证了本文所提方法的角度估计性能与2D-ESPRIT算法估计性能相近。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (4)

1.一种均匀面阵中降维传播算子的二维DOA估计方法，其特征在于，包括以下步骤：

1)构造均匀面阵的信号模型，得到接收信号模型x(t)，构造所述接收信号模型x(t)的协方差矩阵

2)对构造的协方差矩阵进行分块处理，通过传播算子矩阵的旋转不变性进行初始估计，得到用于谱峰搜索的初始估计值/>

3)通过降维处理，将二维优化问题转变为一维优化问题，进行一维局部谱峰搜索获得俯仰角u的估计值

4)根据最小二乘思想获取方位角v的估计值

5)利用估计值和/>得到均匀面阵中的二维DOA精确估计值/>

步骤1)中：

所述均匀面阵共有M×N个均匀分布的阵元，M与N分别为x轴与y轴上的阵元个数，x轴与y轴上的相邻阵元间距均为d，且d≤λ/2，λ为波长，其中将每一列的阵元集合称作一个子阵，则有子阵1,2,…,n,…,N，共N个子阵；

假设某一空间有K个信号源以二维波达方向θ_k,入射到所述均匀面阵上，其中θ_k,/>分别表示第k(k＝1,2,…,K)个信号源的仰角与方位角，定义接收信号x(t)为：

其中：

其中：为Kronecker乘积，/>为方向矩阵A_y和A_x的Khatri-Rao乘积；为信号源矢量，s_K(t)为第K个信号源矢量，(·)^T为转置；n(t)＝n₁(t)^T,n₂(t)^T,…,n_N(t)^T T为所述均匀面阵的加性高斯白噪声，n_n(t)为子阵n的加性高斯白噪声；/>为x轴上的信号源方向向量，/>为y轴上的信号源方向向量；

构造接收信号x(t)的协方差矩阵为：

其中，L为快拍数，(·)^H表示共轭转置；

步骤2)中：

定义参量u和v分别为将x轴与y轴上的信号源方向向量和/>转化为：

a_x(v)＝[1,e^-j2πdv/λ,…,e^{-j2πd(N-1)v/λ}]^T

a_y(u)＝[1,e^-j2πdu/λ,…,e^{-j2πd(M-1)u/λ}]^T

构造矩阵P_c为：

其中，P为传播算子，I_K为K阶单位矩阵，(·)^H为共轭转置；

将矩阵P_c划分为N块：

P_c＝[P_c1,P_c2,…,P_cN]^T

其中，矩阵n＝1,…,N；

分别构造矩阵P₁和P₂：

P₁＝[P_c1,P_c2,…,P_c(N-1)]^T

P₂＝[P_c2,P_c3,…,P_cN]^T

得到：

P₂＝P₁T^-1ΦT

其中，为对角矩阵，u_k为第k个信号源的俯仰角，T为K×K阶的满秩矩阵，然后将P₁ ⁺·P₂进行特征值分解，记得到的第k个特征值为p_k，(·)⁺为求广义逆，·表示相乘，从而得出u_k的初始估计值/>为：

其中k＝1,2,…,K。

2.如权利要求1所述的一种均匀面阵中降维传播算子的二维DOA估计方法，其特征在于，步骤3)中：

定义矩阵Q：

其中，P为传播算子，I_MN-K为MN-K阶单位矩阵；

定义参数V(u,v)：

V(u,v)＝[a_y(u) a_x(v)]^HQQ^H[a_y(u) a_x(v)]

将V(u,v)变形为：

V(u,v)＝a_x(v)^H[a_y(u) I_M]^HQQ^H[a_y(u) I_M]a_x(v)＝a_x(v)^HC(u)a_x(v)

为了便于后续的计算，令参数C(u)＝[a_y(u) I_M]^HQQ^H[a_y(u) I_M]，I_M为M阶单位矩阵；

用e₁ ^Ha_x(v)＝1来消除a_x(v)＝0_M的平凡解，其中矩阵考虑到e₁ ^Ha_x(v)＝1，则将变形之后的V(u,v)重构为：

然后构造代价函数L(u,v)：

L(u,v)＝a_x(v)^HC(u)a_x(v)-ω(e₁ ^Ha_x(v)-1)

其中，ω为一个常数，则有：

则有a_x(v)＝μC(u)^-1e₁，由于e₁ ^Ha_x(v)＝1，则μ＝1/{e₁ ^HC^-1(u)e₁}为常量，所以a_x(v)表示为：

结合得到俯仰角u的估计值/>

将估计值用另一种表示方式为：

3.如权利要求2所述的一种均匀面阵中降维传播算子的二维DOA估计方法，其特征在于，步骤4)中：

在区间内对参量u进行一维局部搜索，得到C^-1(u)中第(1,1)元素中最大的K个值，其中Δu→0；这最大的K个值的峰值(u₁,u₂,…,u_K)便与sinθ_ksinφ_k,(k＝1,2,…,K)的估计值相对应；

然后求出K个导向矢量

由于得出：

g_k＝-angle(a_x(v_k))＝[0,2πdv_k/λ,…,2π(M-1)dv_k/λ]^T＝v_kq

为了便于进一步推导，定义参数q＝[0,2πd/λ,…,2π(M-1)d/λ]^T；归一化再通过最小二乘法求v_K，最小二乘法则为：

其中：E＝[1_M,q]，

c_k1是v_k的估计值；则c_k的解为：

得到v_k的估计值

4.如权利要求3所述的一种均匀面阵中降维传播算子的二维DOA估计方法，其特征在于，步骤5)中：

将和/>代入上式即得到均匀面阵中的二维DOA精确估计值/>