CN111366893B - 一种均匀圆阵未知互耦条件下的非圆信号方位角估计方法 - Google Patents

Abstract

本发明提供了一种均匀圆阵未知互耦条件下的非圆信号方位角估计方法，构造均匀圆阵的互耦系数矩阵，获得输出信号模型，求四阶累计量矩阵，对四阶累计量矩阵进行奇异值分解，求得信号子空间和噪声子空间；对导向矢量做矩阵‑向量变换，最终利用导向矢量和噪声子空间之间的正交性得出空间谱，对空间谱进行降秩分解操作，对空间谱中的方位角进行一维谱峰搜索，找出来波方向角，达到测向的目的。本发明采用矩阵降秩分解变换操作，降低了搜索维度，减少了运算复杂度，有利于工程实现，实时性高。

Description

一种均匀圆阵未知互耦条件下的非圆信号方位角估计方法

技术领域

本发明涉及阵列信号处理领域，具体是一种非圆信号的估计方法。

背景技术

波达方向(direction-of-arrival,DOA)估计技术是阵列信号处理领域的一个重要研究方向。近几十年来，超分辨DOA估计技术因其优越的性能，在雷达、声纳以及通信等领域得到了广泛应用。随着现代无线通信系统的迅速发展，BPSK、PAM以及ASK等常见信号均具有非圆特性。因此，针对非圆信号的DOA测向算法在阵列信号处理领域受到了越来越多的关注，利用信号的非圆特性，阵列可以获得更大的虚拟孔径，从而能够检测更多的信号源，并提高测向精度。

然而，大多数超分辨DOA算法都需要精确的阵列流形，但在实际工程应用中，相邻天线阵元之间往往是存在互耦效应的，这使得真实阵列流形偏离理想(或标称)阵列流形，从而导致超分辨DOA算法性能退化，甚至失效。

由于阵元之间的互耦效应通常随着阵元间距的增大而快速减弱，因此，利用均匀圆阵的圆周对称特性，均匀圆阵的互耦系数矩阵可以表示为具有对称性的Toeplitz矩阵，现有存在对未知互耦影响下的均匀圆阵超分辨DOA算法主要利用这一点进行方位角估计(如文献：‘L.Xiang,Z.Ye,X.Xu,C.Chang,W.Xu,Y.S Hung.“Direction of arrivalestimation for uniform circular array based on fourth-order cumulants in thepresence of unknown mutual coupling”,IET Microw.Antennas Propag.,2008,Vol.2,No.3,pp.

281-287.’)。但是，针对互耦效应影响下的均匀圆阵方位角估计问题，此类方法仍存在如下不足：(1)没有充分利用信号的非圆特性，存在估计的信号源数目少和估计精度低等问题；(2)没有充分利用四阶累计量对于非高斯信号的统计作用，然而在实际场景下大多数信号是呈现非高斯分布的，工程应用会导致算法性能退化(3)没有充分利用扩展的四阶累计量矩阵对于提高角度分辨率的作用。

发明内容

为了克服现有技术的不足，解决现有技术在存在未知互耦影响的均匀圆阵中对方位角的估计问题上存在估计角度精度低、估计信号源数目少等缺陷，本发明提供一种均匀圆阵未知互耦条件下的非圆信号方位角估计方法。本发明通过将四维空间谱函数进行降秩分解，可以有效地降低计算复杂度，实现对来波方向角的高分辨估计；同时，利用信号的非圆特性，可以提高DOA的估计精度并增加估计的信号源数目。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1、利用均匀圆阵的圆周对称性以及阵元之间的互耦强度与距离呈反比的特性，构造均匀圆阵的互耦系数矩阵；基于互耦系数矩阵和均匀圆阵的接收信号模型，构建阵列的导向矢量，获得均匀圆阵的输出信号模型；

步骤1.1：求均匀圆阵的互耦系数矩阵；

假设均匀圆阵有N个阵元，因为均匀圆阵具有圆周对称性，因此互耦系数矩阵为圆周对称的Toeplitz矩阵，具体表达式如下：

其中c_i为第i个阵元与中心阵元之间的互耦系数，其中 为向下取整函数；

步骤1.2：根据均匀圆阵的几何模型，假设第m个信号源到均匀圆阵的方位角为θ_m，其中θ_m∈[0,2π]，第i个阵元与x轴的夹角为该处的位置矢量为当一个波数为/>的窄带平面波在方向为-r上传播，其中λ为波长，r为单位矢量，r＝(cosθ_m,sinθ_m,0)，故在原点与第i个阵元的相位差为/>即/>得出信号源的方位角为θ_m时的导向矢量a(θ_m)为：

对于K个信号源所构成的导向矢量矩阵为：

A＝[a(θ₁),...,a(θ_K)]

步骤1.3：考虑到阵元的输入信号为零均值的平稳的非高斯随机过程，记为s(t)，故s(t)＝[s₁(t),...,s_K(t)]^T，N个阵元输出信号记为x(t)＝[x₁(t),...,x_N(t)]^T，叠加在输出信号上的为复高斯白噪声，记为n(t)＝[n₁(t),...,n_N(t)]^T，故输出信号模型为：

x(t)＝CAs(t)+n(t)

步骤2、对阵列的输出信号x(t)求其3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，得到扩展的四阶累积量矩阵C_x；

步骤2.1：先求3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，假设接收信号分别为由四阶累计量的定义得到：

E{·}代表对任意一组随机变量求期望，cum(·)是对一组随机变量求四阶累计量；

因为输出信号为独立的随机信号，于是输出信号的四阶累积量表示为每个输出信号的四阶累计量的和，定义C_x,1为N²×N²的矩阵，其中第((k₁-1)N+k₂)行第((k₃-1)N+k₄)列的元素为

Ca(θ_i)表示互耦系数矩阵与导向矢量的乘积，其中互耦系数矩阵C与导向矢量a(θ_i)在上文中均已定义过，CA为互耦系数矩阵C与阵列流型矩阵A的乘积，θ_i为来波方向角；

定义C_S为K×K的矩阵，表达式为C_S＝diag{γ₁,γ₂,...,γ_k}，其中γ_i＝cum(s_0,i,s_0,i,s_0,i,s_0,i)；

其中H代表矩阵的共轭转置，代表Kronecker积，/>代表Khatri-Rao积；

同理C_x,2和C_x,3表示为：

步骤2.2：将步骤2.1得到的3个四阶累计量矩阵按如下形式合并为一个扩展的四阶累计量矩阵C_x：

步骤3、对步骤2.2扩展的四阶累计量矩阵进行奇异值分解，求得信号子空间和噪声子空间；对导向矢量做矩阵-向量变换，得到单独含有方位角的矩阵和单独含有互耦系数的向量；最终利用导向矢量和噪声子空间之间的正交性得出空间谱：

步骤3.1：对扩展的四阶累计量矩阵C_x作奇异值分解如下：

其中U_S,U_n分别为信号子空间和噪声子空间，C_x作如下矩阵分解：

记新的导向矢量

步骤3.2：根据互耦系数矩阵的性质，对导向矢量矩阵做如下的矩阵变换：

Ca(θ)＝T(θ)c

其中向量c满足如下表达式：

c₀＝1，c_i＝C_1i(i＝1,2,3,...,M)

其中T(θ)分解为如下表达式：

T＝T₁+T₂+T₃+T₄

定义T(θ)为N×M矩阵，步骤3.3：对非圆信号定义为做如下定义：

其中*取共轭；

将非圆信号写成如下公式：

s(t)＝Φs₀(t)

s₀(t)＝[s_0,1(t),...,s_0.K(t)]^T

其中s_0,i(t)(i＝1,2,...,K)是当s_i(t)的相位为0时的实值信号，对角矩阵Φ定义如下：

其中为s_i(t)的初始相位；

步骤3.4：得出空间谱为如下表达式：

其中||·||代表取二范数；

步骤4、对空间谱进行降秩分解操作，降低计算复杂度；具体步骤如下：

通过将新的导向矢量进行矩阵分解，实现方位角、互耦系数和非圆信初始相位的分离，从而使对使空间谱的四维搜索变成了一维搜索：

其中G(θ)和被定义为如下表达式：

对如下空间谱进行四维谱峰搜索：

因为并不恒为0，得到新的空间谱函数g(θ)为如下表达式：

只需对空间谱g(θ)中的方位角进行一维谱峰搜索，找出来波方向角θ，达到测向的目的，并且不需要思维搜索减少了计算量。

本发明的有益效果在于：

第一：采用矩阵降秩分解变换操作，降低了搜索维度，减少了运算复杂度，有利于工程实现，实时性高。

第二，相较于论文‘M.Lin,L.Yang,“Blind calibration and DOA estimationwith uniform circular arrays in the presence of mutual coupling”,IEEEAntennas Wireless Propag.Lett.,vol.5,pp.315–318,2006.’的方法，本发明提出的方法的优势是估计的信号源数目是其2倍以上而且角度分辨率更高在信号服从非高斯分布的情况下，在强互耦条件下本方法也能实现对来波方位角的高分辨估计。

第三，相较于论文‘L.Xiang,Z.Ye,X.Xu,C.Chang,W.Xu,Y.S Hung.“Direction ofarrival estimation for uniform circular array based on fourth-order cumulantsin the presence of unknown mutual coupling”,IET Microw.Antennas Propag.,2008,Vol.2,No.3,pp.281-287.’的方法，本发明充分利用了非圆信号信息，对于在均匀圆阵中对于非圆信号的互耦测向有着更高的角度分辨率和能估计更多的信号源数目。

附图说明

图1是本发明非圆信号方位角估计方法的流程图。

图2是本发明均匀圆阵的几何模型示意图。

图3是本发明与传统方法对于10个BPSK信号的空间谱。

图4是本发明与传统方法下RMSE随信噪比的变化对比图。

图5是本发明与传统方法下均方根误差(RMSE)对比图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

为使本发明的目的、技术方案和优点更加清楚，下面结合附图和实施例对本发明作进一步详细描述，本发明包括但不仅限于下述实施例。

1、构造均匀圆阵的输出信号模型，求出在互耦条件下的非圆信号激励下阵列的输出信号。

2、求输出信号的扩展四阶累计量矩阵，并对该矩阵进行奇异值分解得出信号子空间和噪声子空间

3、利用噪声子空间与阵列流型正交的关系的出空间谱函数，对空间谱实施变量分离得出只含有方位角的空间谱函数，对空间谱函数进行一维谱峰搜索，得出估计的来波方向角。

步骤1、利用均匀圆阵的圆周对称性以及阵元之间的互耦强度与距离呈反比的特性，构造均匀圆阵的互耦系数矩阵；基于互耦系数矩阵和均匀圆阵的接收信号模型，构建阵列的导向矢量，获得均匀圆阵的输出信号模型；

步骤1.1：求均匀圆阵的互耦系数矩阵；

假设均匀圆阵有N个阵元，因为均匀圆阵具有圆周对称性，因此互耦系数矩阵为圆周对称的Toeplitz矩阵，具体表达式如下：

其中c_i为第i个阵元与中心阵元之间的互耦系数，其中 为向下取整函数；

步骤1.2：根据图2均匀圆阵的几何模型，假设第m个信号源到均匀圆阵的方位角为θ_m，其中θ_m∈[0,2π]，第i个阵元与x轴的夹角为该处的位置矢量为当一个波数为/>的窄带平面波在方向为-r上传播，其中λ为波长，r为单位矢量，r＝(cosθ_m,sinθ_m,0)，故在原点与第i个阵元的相位差为/>即/>得出信号源的方位角为θ_m时的导向矢量a(θ_m)为：

对于K个信号源所构成的导向矢量矩阵为：

A＝[a(θ₁),...,a(θ_K)]

步骤1.3：考虑到阵元的输入信号为零均值的平稳的非高斯随机过程，记为s(t)，故s(t)＝[s₁(t),...,s_K(t)]^T，N个阵元输出信号记为x(t)＝[x₁(t),...,x_N(t)]^T，叠加在输出信号上的为复高斯白噪声，记为n(t)＝[n₁(t),...,n_N(t)]^T，故输出信号模型为：

x(t)＝CAs(t)+n(t)

步骤2、对阵列的输出信号x(t)求其3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，得到扩展的四阶累积量矩阵C_x；

步骤2.1：先求3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，假设接收信号分别为由四阶累计量的定义得到：

E{·}代表对任意一组随机变量求期望，cum(·)是对一组随机变量求四阶累计量；

因为输出信号为独立的随机信号，于是输出信号的四阶累积量表示为每个输出信号的四阶累计量的和，定义C_x,1为N²×N²的矩阵，其中第((k₁-1)N+k₂)行第((k₃-1)N+k₄)列的元素为

Ca(θ_i)表示互耦系数矩阵与导向矢量的乘积，其中互耦系数矩阵C与导向矢量a(θ_i)在上文中均已定义过，CA为互耦系数矩阵C与阵列流型矩阵A的乘积，θ_i为来波方向角；

上式中，定义C_S为K×K的矩阵，表达式为C_S＝diag{γ₁,γ₂,...,γ_k}，其中γ_i＝cum(s_0,i,s_0,i,s_0,i,s_0,i)；

其中H代表矩阵的共轭转置，代表Kronecker积，/>代表Khatri-Rao积，；

同理C_x,2和C_x,3表示为：

步骤2.2：将步骤2.1得到的3个四阶累计量矩阵按如下形式合并为一个扩展的四阶累计量矩阵C_x：

步骤3、对步骤2.2扩展的四阶累计量矩阵进行奇异值分解，求得信号子空间和噪声子空间；对导向矢量做矩阵-向量变换，得到单独含有方位角的矩阵和单独含有互耦系数的向量；最终利用导向矢量和噪声子空间之间的正交性得出空间谱：

步骤3.1：对扩展的四阶累计量矩阵C_x作奇异值分解如下：

其中U_S,U_n分别为信号子空间和噪声子空间，C_x作如下矩阵分解：

记新的导向矢量

步骤3.2：根据互耦系数矩阵的性质，对导向矢量矩阵做如下的矩阵变换：

Ca(θ)＝T(θ)c

其中向量c满足如下表达式：

c₀＝1，c_i＝C_1i(i＝1,2,3,...,M)

其中T(θ)分解为如下表达式：

T＝T₁+T₂+T₃+T₄

定义T(θ)为N×M矩阵，

步骤3.3：对非圆信号定义为做如下定义：

其中*取共轭；

为了对非圆信号进行统计分析，将非圆信号写成如下公式：

s(t)＝Φs₀(t)

s₀(t)＝[s_0,1(t),...,s_0.K(t)]^T

其中s_0,i(t)(i＝1,2,...,K)是当s_i(t)的相位为0时的实值信号，对角矩阵Φ定义如下：

其中为s_i(t)的初始相位；

步骤3.4：得出空间谱为如下表达式：

其中||·||代表取二范数；

步骤4、对空间谱进行降秩分解操作，降低计算复杂度；具体步骤如下：

通过将新的导向矢量进行矩阵分解，实现方位角、互耦系数和非圆信初始相位的分离，从而使对使空间谱的四维搜索变成了一维搜索：

其中G(θ)和被定义为如下表达式：

对如下空间谱进行四维谱峰搜索：

因为并不恒为0，得到新的空间谱函数g(θ)为如下表达式：

只需对空间谱g(θ)中的方位角进行一维谱峰搜索，就能找出来波方向角θ，达到测向的目的并且不需要思维搜索减少了计算量。

实施例的具体步骤如下：

步骤一：求均匀圆阵的输出信号

(1)构造均匀圆阵的互耦系数矩阵C：

(2)求均匀圆阵的导向矢量a(θ_m)和K个信号源所构成的导向矢量矩阵A：

A＝[a(θ₁),...,a(θ_K)]

(3)求出均匀圆阵的输出信号：

x(t)＝CAs(t)+n(t)

步骤二：求输出信号的扩展四阶累计量矩阵，并对该矩阵进行奇异值分解得出信号子空间和噪声子空间

(1)求3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3：

(2)求扩展的四阶累计量矩阵C_x：

(3)对扩展的四阶累计量矩阵C_x进行奇异值分解和矩阵分解：

步骤三：求降低计算复杂度后的空间谱函数

(1)求新的导向矢量

(2)得出空间谱函数

(3)得出最终只含有方位角的空间谱函数g(θ)：

仿真条件说明：本次仿真实验平台在Windows 10操作系统中的MATLAB R2018a进行。本次仿真实验是在含有5个阵元的均匀圆阵中进行，信噪比为10dB，快拍数为1024，互耦系数为c＝[1,0.37+j*0.42,0.09+j*0.21]^T，实际信源的波达方向为10°,48°,86°,123°,161°,199°,237°,274°,312°,350°。

从图2可以看出所得出的空间谱峰尖锐的指向真正的来波方向角，相比传统方法有着更高的角度分辨率，同时也可以看出估计的信号源数目也比传统方法要多。

图3是两种方法对于10个BPSK信号的空间谱，其中传统方法是指论文‘L.Xiang,Z.Ye,X.Xu,C.Chang,W.Xu,Y.S Hung.“Direction of arrival estimation for uniformcircular array based on fourth-order cumulants in the presence of unknownmutual coupling”,IET Microw.Antennas Propag.,2008,Vol.2,No.3,pp.281-287.’的方法，即为传统的估计方法，从图3可以看出，在快拍数固定为1024，信噪比从0dB到30dB变化时，本方法对于角度估计的均方根误差比传统方法要小很多，而且随着信噪比的提高，本方法的均方根误差显著下降。

图4是RMSE随信噪比的变化示意图，当快拍数一定时，随着信噪比从0dB到30dB变化时，本方法与传统估计方法估计角度的均方根误差(RMSE)对比如图4所示；图4可以看出，在信噪比固定为20dB时，本方法对于角度估计的均方根误差比传统方法要小很多，而且随着快拍数的提高，本方法的均方根误差显著下降。

图5是均方根误差(RMSE)对比图，当信噪比一定时，随着快拍数变化时，本方法与传统估计方法估计角度的均方根误差(RMSE)对比图如图5所示；

为了使本领域普通技术人员理解本发明，而对本发明进行了详细描述，但可以想到，在不脱离本发明的权利要求所涵盖的范围内还可以做出其他的变化和修改，这些变化和修改均在本发明的保护范围之内。

Claims (1)

1.一种均匀圆阵未知互耦条件下的非圆信号方位角估计方法，其特征在于包括下述步骤：

步骤1、利用均匀圆阵的圆周对称性以及阵元之间的互耦强度与距离呈反比的特性，构造均匀圆阵的互耦系数矩阵；基于互耦系数矩阵和均匀圆阵的接收信号模型，构建阵列的导向矢量，获得均匀圆阵的输出信号模型；

步骤1.1：求均匀圆阵的互耦系数矩阵；

假设均匀圆阵有N个阵元，因为均匀圆阵具有圆周对称性，因此互耦系数矩阵为圆周对称的Toeplitz矩阵，具体表达式如下：

其中c_i为第i个阵元与中心阵元之间的互耦系数，其中 为向下取整函数；

步骤1.2：根据均匀圆阵的几何模型，假设第m个信号源到均匀圆阵的方位角为θ_m，其中θ_m∈[0,2π]，第i个阵元与x轴的夹角为该处的位置矢量为/>当一个波数为/>的窄带平面波在方向为-r上传播，其中λ为波长，r为单位矢量，r＝(cosθ_m,sinθ_m,0)，故在原点与第i个阵元的相位差为/>即/>得出信号源的方位角为θ_m时的导向矢量a(θ_m)为：

对于K个信号源所构成的导向矢量矩阵为：

A＝[a(θ₁),...,a(θ_K)]

步骤1.3：考虑到阵元的输入信号为零均值的平稳的非高斯随机过程，记为s(t)，故s(t)＝[s₁(t),...,s_K(t)]^T，N个阵元输出信号记为x(t)＝[x₁(t),...,x_N(t)]^T，叠加在输出信号上的为复高斯白噪声，记为n(t)＝[n₁(t),...,n_N(t)]^T，故输出信号模型为：

x(t)＝CAs(t)+n(t)

步骤2、对阵列的输出信号x(t)求其3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，得到扩展的四阶累积量矩阵C_x；

对非圆信号定义为做如下定义：

其中*取共轭；

将非圆信号写成如下公式：

s(t)＝Φs₀(t)

s₀(t)＝[s_0,1(t),...,s_0.K(t)]^T

其中s_0,j'(t),j'＝1,2,...,K是当s_i(t)的相位为0时的实值信号，对角矩阵Φ定义如下：

其中为第i个阵元与x轴的夹角；

步骤2.1：先求3个不同的四阶累计量矩阵C_x,1，C_x,2和C_x,3，假设接收信号分别为由四阶累计量的定义得到：

E{·}代表对任意一组随机变量求期望，cum(·)是对一组随机变量求四阶累计量；

因为输出信号为独立的随机信号，于是输出信号的四阶累积量表示为每个输出信号的四阶累计量的和，定义C_x,1为N²×N²的矩阵，其中第((k₁-1)N+k₂)行第((k₃-1)N+k₄)列的元素为

Ca(θ_i)表示互耦系数矩阵与导向矢量的乘积，其中互耦系数矩阵C与导向矢量a(θ_i)在上文中均已定义过，CA为互耦系数矩阵C与阵列流型矩阵A的乘积，θ_i为来波方向角；

定义C_S为K×K的矩阵，表达式为C_S＝diag{γ₁,γ₂,...,γ_k}，其中γ_j'＝cum(s_0,j',s_0,j',s_0,j',s_0,j')；

其中H代表矩阵的共轭转置，代表Kronecker积，/>代表Khatri-Rao积；

同理C_x,2和C_x,3表示为：

步骤2.2：将步骤2.1得到的3个四阶累计量矩阵按如下形式合并为一个扩展的四阶累计量矩阵C_x：

步骤3、对步骤2.2扩展的四阶累计量矩阵进行奇异值分解，求得信号子空间和噪声子空间；对导向矢量做矩阵-向量变换，得到单独含有方位角的矩阵和单独含有互耦系数的向量；最终利用导向矢量和噪声子空间之间的正交性得出空间谱：

步骤3.1：对扩展的四阶累计量矩阵C_x作奇异值分解如下：

其中U_S,U_n分别为信号子空间和噪声子空间，C_x作如下矩阵分解：

记新的导向矢量

步骤3.2：根据互耦系数矩阵的性质，对导向矢量矩阵做如下的矩阵变换：

Ca(θ)＝T(θ)c

其中向量c满足如下表达式：

c₀＝1，c_i'＝C_1i'，i’＝1,2,3,...,M；

其中T(θ)分解为如下表达式：

T＝T₁+T₂+T₃+T₄

定义T(θ)为N×M矩阵，

步骤3.4：得出空间谱为如下表达式：

其中||·||代表取二范数；

步骤4、对空间谱进行降秩分解操作，降低计算复杂度；具体步骤如下：

通过将新的导向矢量进行矩阵分解，实现方位角、互耦系数和非圆信初始相位的分离，从而使对使空间谱的四维搜索变成了一维搜索：

其中G(θ)和被定义为如下表达式：

对如下空间谱进行四维谱峰搜索：

因为并不恒为0，得到新的空间谱函数g(θ)为如下表达式：

只需对空间谱g(θ)中的方位角进行一维谱峰搜索，找出来波方向角θ，达到测向的目的，并且不需要思维搜索减少了计算量。