CN106019234A - L型天线阵的低计算复杂度二维波达方向估计方法 - Google Patents

Abstract

本发明公开了一种L型天线阵的低计算复杂度二维波达方向估计方法，主要解决现有L型天线阵二维波达方向估计方法计算复杂度高且在低信噪比、低快拍情况下估计性能差的问题，其实现步骤是：(1)接收信号；(2)构造互相关矩阵；(3)构造方位增广矩阵；(4)构造俯仰增广矩阵；(5)估计方位角；(6)估计俯仰角；(7)二维角度配对。本发明利用旋转不变子空间算法和对角度托普利兹矩阵进行特征值分解，以低计算复杂度、低信噪比、低快拍，实现了L型天线阵二维波达方向估计，可用于雷达、通信中的目标定位。

Description

L型天线阵的低计算复杂度二维波达方向估计方法

技术领域

本发明属于雷达技术领域，更进一步涉及雷达技术领域中的一种L型天线阵的低计算复杂度二维波达方向估计方法。本发明可用于雷达天线在低信噪比、低快拍情况下对目标信号的快速二维波达方向DOA(Direction of Arrival)的估计，提高雷达天线在低信噪比、低快拍情况下对目标信号的快速测向性能。

背景技术

二维DOA估计在雷达、通信等领域应用广泛。L型天线阵作为一种可用于二维DOA估计的阵列，已被证明具有很好的二维DOA估计性能。最近数十年以来，针对L型天线阵的二维DOA估计已经有了很多的研究成果。

N.Xi等人在其发表的论文“A Computationally Efficient Subspace Algorithmfor 2-D DOA Estimation with L-shaped Array”(《2014IEEE Signal ProcessingLetters》2014,pp:971-974)中提出了CESA方法。该方法的具体步骤是，第一步：采用2M个天线构成的L型天线阵接收目标信号，对每个天线接收到的目标信号分别按照奈奎斯特采样定理进行采样，得到L型天线阵数据矩阵；第二步：将L型天线阵数据矩阵中沿水平方向放置的天线的数据划分为x子阵数据矩阵，沿竖直方向放置的天线的数据划分为z子阵数据矩阵；第三步：利用L型天线阵的x子阵数据矩阵和z子阵数据矩阵，构造互相关矩阵；第四步：利用互相关矩阵，构造正交方位信号矩阵和正交俯仰信号矩阵；第五步：利用正交方位信号矩阵和正交俯仰信号矩阵，构造方位搜索函数和俯仰搜索函数，并对方位搜索函数和俯仰搜索函数进行谱峰搜索，得到方位角与俯仰角；第六步：利用最大互斥角度配对方法，得到二维角度配对的方位角和俯仰角。该方法虽然相对于已有方法运算量小、估计精度高，但是，该方法仍然存在的不足之处是，需要两次一维角度搜索，且在低信噪比、低快拍情况下的角度估计性能很差。

西安交通大学在其申请的专利“自动配准的二维波达方向估计装置及其方法”(申请号：CN201110077438；申请公布号：CN102142879A)公开了一种自动配准的二维波达方向估计装置及其方法。该方法的具体步骤是，第一步：采用2M个天线构成的L型天线阵接收目标信号，对每个天线接收到的目标信号分别按照奈奎斯特采样定理 进行采样，得到L型天线阵数据矩阵。第二步：将L型天线阵数据矩阵中沿水平方向放置的天线的数据划分为x子阵数据矩阵，沿竖直方向放置的天线的数据划分为z子阵数据矩阵；第三步：利用z子阵数据矩阵，构造俯仰搜索函数，并对俯仰搜索函数进行谱峰搜索，得到方位角；第四步，利用x子阵数据矩阵和z子阵数据矩阵，构造配对搜索函数，并对配对搜索函数进行谱峰搜索，得到俯仰角。该方法虽然能够实现自动配准的二维波达方向估计，但是，该方法仍然存在的不足之处是，需要多次大量的一维角度搜索。

发明内容

本发明的目的在于针对现有技术的L型天线阵二维DOA估计方法计算复杂度高且在低信噪比、低快拍情况下估计性能差的问题，提出一种可在低信噪比、低快拍情况下实现对目标信号的快速二维波达方向估计的方法。

实现本发明的目的的思路是，首先利用L型天线阵接收目标信号，生成接收数据矩阵，并构造互相关矩阵。然后利用互相关矩阵分别构造方位增广矩阵和俯仰增广矩阵，并对方位增广矩阵和俯仰增广矩阵分别利用旋转不变子空间算法，得到方位增广矩阵的方位角和俯仰增广矩阵的俯仰角。最后利用最大互斥角度配对方法，将方位增广矩阵的方位角和俯仰增广矩阵的俯仰角代入角度配对代价函数，得到二维角度配对的L型天线阵接收目标信号的方位角和L型天线阵接收目标信号俯仰角。

本发明的具体实施步骤如下：

(1)接收信号：

(1a)由2M个天线构成的L型天线阵接收目标信号，对每个天线接收到的目标信号分别按照奈奎斯特采样定理进行采样，得到L型天线阵数据矩阵；

(1b)将L型天线阵数据矩阵中沿水平方向放置的天线的数据划分为x子阵数据矩阵，沿竖直方向放置的天线的数据划分为z子阵数据矩阵；

(2)构造互相关矩阵：

将x子阵数据矩阵和z子阵数据矩阵的共轭转置相乘，得到互相关矩阵；

(3)按照下式，构造方位增广矩阵：

其中，E表示方位增广矩阵，R表示互相关矩阵，J_M表示M×M维反对角线元 素全为1而其它元素为0的交换矩阵，*表示矩阵共轭操作，表示将两个矩阵按行组成一个新矩阵的操作；

(4)按照下式，构造俯仰增广矩阵：

其中，F表示俯仰增广矩阵，R表示互相关矩阵，J_M表示M×M维反对角线元素全为1而其它元素为0的交换矩阵，H表示共轭转置操作，T表示转置操作，表示将两个矩阵按行组成一个新矩阵的操作；

(5)估计方位角：

利用旋转不变子空间ESPRIT算法，估计方位增广矩阵的方位角；

(6)估计俯仰角：

利用旋转不变子空间ESPRIT算法，估计俯仰增广矩阵的俯仰角；

(7)二维角度配对：

(7a)利用互相关矩阵，构造角度相关向量；

(7b)利用角度相关向量，构造角度托普利兹矩阵；

(7c)按照下式，对角度托普利兹矩阵进行特征值分解，得到立体信号特征矩阵和立体信号特征值向量：

{U,h}＝EVD(Ξ)

其中，U表示立体信号特征矩阵，h表示立体信号特征值向量，EVD表示特征值分解操作，Ξ表示角度托普利兹矩阵；

(7d)按照下式，构造角度配对代价函数：

其中，f表示角度配对代价函数，a表示立体搜索导向矢量，H表示共轭转置操作，U表示立体信号特征矩阵；

(7e)利用最大互斥角度配对方法，将步骤(5)得到的方位增广矩阵的方位角和步骤(6)得到的俯仰增广矩阵的俯仰角代入角度配对代价函数，得到二维角度配 对的L型天线阵接收目标信号的方位角和L型天线阵接收目标信号俯仰角。

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

第一，由于本发明利用了旋转不变子空间ESPRIT算法，估计方位增广矩阵的方位角和俯仰增广矩阵的俯仰角，克服了现有技术因为角度搜索造成的复杂度高的问题，使得本发明可以在不进行角度搜索的情况下，以低计算复杂度，实现二维DOA估计。

第二，由于本发明利用了互相关矩阵的特点，构造了方位增广矩阵和俯仰增广矩阵，克服了现有技术二维DOA估计时可用数据少而造成的估计精度差的问题，使得本发明在可用数据少的情况下，以更高的估计精度，实现二维DOA估计。

第三，由于本发明对角度托普利兹矩阵进行特征值分解，得到立体信号特征矩阵和立体信号特征值向量，克服了现有技术信噪比低和快拍数低不能进行二维角度配对的问题，使得本发明对信噪比和快拍数的需求降低，以低信噪比、低快拍，实现二维DOA估计。

附图说明：

图1是本发明的流程图；

图2是分别利用本发明和现有技术的CESA算法得到的二维DOA估计的均方根误差随信噪比变化的曲线图；

图3是分别利用本发明和现有技术的CESA算法得到的二维DOA估计的均方根误差随快拍数变化的曲线图；

图4是分别利用本发明和现有技术的CESA算法统计平均运行时间随子阵阵元数M变化的曲线图。

具体实施方式：

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1，接收信号。

由2M个天线构成的L型天线阵接收目标信号，对每个天线接收到的目标信号分别按照奈奎斯特采样定理进行采样，得到L型天线阵数据矩阵。

将L型天线阵数据矩阵中沿水平方向放置的天线的数据划分为x子阵数据矩阵，沿竖直方向放置的天线的数据划分为z子阵数据矩阵，具体可以表示为：

X＝AS+N

Z＝BS+W

其中，X表示x子阵数据矩阵，Z表示z子阵数据矩阵，A表示x子阵流型矩阵，B表示z子阵流型矩阵，S表示信号矩阵，N表示x子阵噪声矩阵，W表示z子阵噪声矩阵。

所述的x子阵流型矩阵的具体形式如下：

其中，A表示x子阵流型矩阵，d表示L型天线阵的x子阵中相邻天线的间距，θ_q表示第q个目标信号方位角，q＝1,2,…,Q，Q表示目标信号数目，λ表示目标信号波长，m表示L型天线阵的x子阵中天线的序号，m＝1,2,…,M，M表示L型天线阵中x子阵的天线数。

z子阵流型矩阵具体形式如下：

其中，B表示z子阵流型矩阵，g表示L型天线阵的z子阵中相邻天线的间距，φ_q表示第q个目标信号俯仰角，q＝1,2,…,Q，Q表示目标信号数目，λ表示目标信号波长，n表示L型天线阵的z子阵中天线的序号，n＝1,2,…,N，N表示L型天线阵中z子阵的天线数。

在本发明的具体实施例中，要求L型天线阵中x子阵的天线数和L型天线阵中z子阵的天线数相等，L型天线阵的x子阵中相邻天线的间距和L型天线阵的z子阵中相邻天线的间距相等。

步骤2，构造互相关矩阵。

将x子阵数据矩阵和z子阵数据矩阵的共轭转置相乘，得到互相关矩阵：

R＝XZ^H

其中，R表示互相关矩阵，X表示x子阵数据矩阵，Z表示z子阵数据矩阵，H表示共轭转置操作。

步骤3，按照下式，构造方位增广矩阵：

其中，E表示方位增广矩阵，R表示互相关矩阵，J_M表示M×M维反对角线元素全为1而其它元素为0的交换矩阵，*表示矩阵共轭操作，表示将两个矩阵按行组成一个新矩阵的操作。

步骤4，按照下式，构造俯仰增广矩阵：

其中，F表示俯仰增广矩阵，R表示互相关矩阵，H表示共轭转置操作，J_M表示M×M维反对角线元素全为1而其它元素为0的交换矩阵，T表示转置操作，表示将两个矩阵按行组成一个新矩阵的操作。

步骤5，估计方位角。

利用旋转不变子空间ESPRIT算法，估计方位增广矩阵的方位角。具体实现如下：

按照下式，利用方位增广矩阵，得到方位增广协方差矩阵：

其中，表示方位增广协方差矩阵，E表示方位增广矩阵。

按照下式，对方位增广协方差矩阵进行特征值分解，得到方位增广特征矩阵和方位增广特征值向量：

其中，C表示方位增广特征矩阵，k表示方位增广特征值向量，EVD表示特征值分解操作，表示方位增广协方差矩阵。

按照下式，构造上选择矩阵：

其中，T表示上选择矩阵，I_(M-1)×(M-1)表示(M-1)×(M-1)维单位阵、0_(M-1)×(M-1)表示(M-1)×(M-1)维零矩阵，0_(M-1)×1表示(M-1)×1维零向量。

按照下式，构造下选择矩阵：

其中，Γ表示下选择矩阵，0_(M-1)×1表示(M-1)×1维零向量，I_(M-1)×(M-1)表示(M-1)×(M-1)维单位阵，0_(M-1)×(M-1)表示(M-1)×(M-1)维零矩阵。

按照下式，构造方位同征矩阵:

其中，Δ表示方位同征矩阵，T表示上选择矩阵，C表示方位增广特征矩阵， 表示伪逆操作，Γ表示下选择矩阵。

按照下式，对方位同征矩阵进行特征值分解，得到方位特征矩阵与方位特征值向量：

{G,α}＝EVD(Δ)

其中，G表示方位特征矩阵，α表示方位特征值向量，EVD表示特征值分解操作，Δ表示方位同征矩阵。

按照下式，估计方位增广矩阵的方位角：

其中，表示第q个方位增广矩阵的方位角，q＝1,2,…,Q，Q表示目标信号数目，angle表示取相角操作，α表示方位特征值向量，α(q)表示方位特征值向量α的第q个元素，λ表示目标信号波长，d表示L型天线阵的x子阵中相邻天线的间距。

步骤6，估计俯仰角。

利用旋转不变子空间ESPRIT算法，估计俯仰增广矩阵的俯仰角。具体实现如下：

按照下式，利用俯仰增广矩阵，得到俯仰增广协方差矩阵：

其中，表示俯仰增广协方差矩阵，F表示俯仰增广矩阵。

按照下式，对俯仰增广协方差矩阵进行特征值分解，得到俯仰增广特征矩阵和俯仰增广特征值向量：

其中，P表示俯仰增广特征矩阵，t表示俯仰增广特征值向量，EVD表示特征值分解操作，表示俯仰增广协方差矩阵。

按照下式，构造俯仰同征矩阵:

其中，Π表示俯仰同征矩阵，T表示上选择矩阵，P表示俯仰增广特征矩阵， 表示伪逆操作，Γ表示下选择矩阵。

按照下式，对俯仰同征矩阵进行特征值分解，得到俯仰特征矩阵与俯仰特征值向量：

{H,β}＝EVD(Π)

其中，H表示俯仰特征矩阵，β表示俯仰特征值向量，EVD表示特征值分解操作，Π表示俯仰同征矩阵。

按照下式，估计俯仰增广矩阵的俯仰角：

其中，表示第q个俯仰增广矩阵的俯仰角，q＝1,2,…,Q，Q表示目标信号数目，angle表示取相角操作，β表示俯仰特征值向量，β(q)表示俯仰特征值向量β的第q个元素，λ表示目标信号波长，g表示L型天线阵的z子阵中相邻天线的间距。

步骤7，二维角度配对。

按照下式，利用互相关矩阵，构造角度相关向量：

r＝[R(1,1),R(2,2),…,R(m,m),…,R(M,M)]

其中，r表示角度相关向量，R表示互相关矩阵，R(m,m)表示互相关矩阵R的第 m行第m列元素，m＝1,2,…,M，M表示L型天线阵中x子阵的天线数。

按照下式，利用角度相关向量，构造角度托普利兹矩阵：

Ξ＝Toeplitz(r)

其中，Ξ表示角度托普利兹矩阵，Toeplitz(·)表示构造托普利兹矩阵操作，r表示角度相关向量。

按照下式，对角度托普利兹矩阵进行特征值分解，得到立体信号特征矩阵和立体信号特征值向量：

{U,h}＝EVD(Ξ)

其中，U表示立体信号特征矩阵，h表示立体信号特征值向量，EVD表示特征值分解操作，Ξ表示角度托普利兹矩阵。

按照下式，构造角度配对代价函数：

其中，f表示角度配对代价函数，a表示立体搜索导向矢量，H表示共轭转置操作，U表示立体信号特征矩阵。

立体搜索导向矢量a是按照如下公式构造的：

a＝[1,exp{j2πd(cosφ-cosθ)/λ},…,exp{j2π(M-1)d(cosφ-cosθ)/λ}]^T

其中，a表示立体搜索导向矢量，exp{·}表示指数函数操作，j表示虚数单位，d表示L型天线阵的x子阵中相邻天线的间距，φ表示任意的俯仰角，θ表示任意的方位角，λ表示目标信号波长，M表示L型天线阵中x子阵的天线数，T表示转置操作。

利用最大互斥角度配对方法，将步骤(5)得到的方位增广矩阵的方位角和步骤(6)得到的俯仰增广矩阵的俯仰角代入角度配对代价函数，并将角度配对代价函数值从大到小排列，选择前Q个角度配对代价函数值对应的方位增广矩阵的方位角和俯仰增广矩阵的俯仰角，作为二维角度配对的L型天线阵接收目标信号的方位角和L型天线阵接收目标信号俯仰角，其中Q表示目标信号数目。至此完成L型天线阵的低计算复杂度二维波达方向估计，得到二维角度配对的L型天线阵接收目标信号的方位角和L型天线阵接收目标信号俯仰角。

下面结合仿真实验对本发明的效果做进一步说明。

1.仿真条件：

本发明的仿真实验中计算机配置环境为Intel(R)Core(i5-3470)3.20GHZ中央处理器、内存8G、WINDOWS 7操作系统，计算机仿真软件采用MATLAB R2013a软件。

本发明的仿真参数如下：假定L型天线阵中x子阵的天线数为8，L型天线阵中z子阵的天线数为8，L型天线阵的x子阵中相邻天线的间距为1，L型天线阵的z子阵中相邻天线的间距为1，目标信号波长为2，三个目标信号的方位角和俯仰角分别为：(80°,55°)、(90°,85°)、(100°,70°)。CESA算法角度搜索范围为0°～180°，角度搜索间隔为0.1^°。二维DOA估计均方根误差RMSE(Root Mean Square Error)表示为：

其中Ω表示蒙特卡洛仿真实验次数，Q表示目标信号数目，表示第η次蒙特卡洛仿真实验得到的第q个目标信号方位角，θ_q表示第q个目标信号方位角，表示第η次蒙特卡洛仿真实验得到的第q个目标信号俯仰角，φ_q表示第q个目标信号俯仰角。

2.仿真内容：

(1)当快拍数为512时，在不同信噪比条件下，采用本发明和现有技术CESA算法进行二维DOA估计，每个信噪比下分别进行500次蒙特卡洛仿真实验，得到二维DOA估计均方根误差，仿真结果如图2所示。

(2)当信噪比为10分贝时，在不同快拍数下，采用本发明和现有技术CESA算法进行二维DOA估计，每个快拍数下分别进行500次蒙特卡洛实验，得到二维DOA估计均方根误差，仿真结果如图3所示。

(3)当信噪比为0分贝且快拍数为512时，在不同x子阵天线数下，采用本发明和现有技术CESA算法进行二维DOA估计，统计两种算法在不同x子阵天线数下的平均运行时间，仿真结果如图4所示。

3.仿真结果分析：

图2是本发明和现有技术CESA算法的二维DOA估计均方根误差随信噪比变化的曲线图，图2中x轴表示信噪比，单位为分贝，图2中z轴表示二维DOA估 计均方根误差，单位为度，图2中以圆形标示的曲线表示本发明的二维DOA估计均方根误差随信噪比变化的曲线，图2中以三角形标示的曲线表示现有技术CESA算法的二维DOA估计均方根误差随信噪比变化的曲线。从图2可以看出，与现有技术CESA算法相比，本发明具有更小的二维DOA估计均方根误差，表明本发明可以在低信噪比情况下进行二维DOA估计，同时具有更高的二维DOA估计精度。

图3是本发明和现有技术CESA算法的二维DOA估计均方根误差随快拍数变化的曲线图，图3中x轴表示快拍数，图3中z轴表示二维DOA估计均方根误差，单位为度，图3中以圆形标示的曲线表示本发明的二维DOA估计均方根误差随快拍数变化的曲线，图3中以三角形标示的曲线表示现有技术CESA算法的二维DOA估计均方根误差随快拍数变化的曲线。从图3可以看出，与现有技术CESA算法相比，本发明具有更小的二维DOA估计均方根误差，表明本发明可以在低快拍情况下进行二维DOA估计，同时具有更高的二维DOA估计精度。

图4是本发明和现有技术CESA算法的平均运行时间随x子阵天线数变化的曲线图，图4中x轴表示x子阵天线数，图4中z轴表示平均运行时间，单位为秒，图4中以圆形标示的曲线表示本发明的平均运行时间随x子阵天线数变化的曲线，图4中以三角形标示的曲线表示现有技术CESA算法的平均运行时间随x子阵天线数变化的曲线。从图4可以看出，与现有技术CESA算法相比，本发明具有更小的平均运行时间，表明本发明能够以更低的计算复杂度实现二维DOA估计。

Claims (4)

1.一种L型天线阵的低计算复杂度二维波达方向估计方法，包括如下步骤：

(1)接收信号：

(1a)由2M个天线构成的L型天线阵接收目标信号，对每个天线接收到的目标信号分别按照奈奎斯特采样定理进行采样，得到L型天线阵数据矩阵；

(1b)将L型天线阵数据矩阵中沿水平方向放置的天线的数据划分为x子阵数据矩阵，沿竖直方向放置的天线的数据划分为z子阵数据矩阵；

(2)构造互相关矩阵：

将x子阵数据矩阵和z子阵数据矩阵的共轭转置相乘，得到互相关矩阵；

(3)按照下式，构造方位增广矩阵：

$E=\left[\begin{array}{c}R\\ J_M{R}^{*}{J}_M\end{array}\right]$

其中，E表示方位增广矩阵，R表示互相关矩阵，J_M表示M×M维反对角线元素全为1而其它元素为0的交换矩阵，*表示矩阵共轭操作，表示将两个矩阵按行组成一个新矩阵的操作；

(4)按照下式，构造俯仰增广矩阵：

$F=\left[\begin{array}{c}{R}^H\\ J_M{R}^T{J}_M\end{array}\right]$

其中，F表示俯仰增广矩阵，R表示互相关矩阵，J_M表示M×M维反对角线元素全为1而其它元素为0的交换矩阵，H表示共轭转置操作，T表示转置操作，表示将两个矩阵按行组成一个新矩阵的操作；

(5)估计方位角：

利用旋转不变子空间ESPRIT算法，估计方位增广矩阵的方位角；

(6)估计俯仰角：

利用旋转不变子空间ESPRIT算法，估计俯仰增广矩阵的俯仰角；

(7)二维角度配对：

(7a)利用互相关矩阵，构造角度相关向量；

(7b)利用角度相关向量，构造角度托普利兹矩阵；

(7c)按照下式，对角度托普利兹矩阵进行特征值分解，得到立体信号特征矩阵和立体信号特征值向量：

{U,h}＝EVD(Ξ)

其中，U表示立体信号特征矩阵，h表示立体信号特征值向量，EVD表示特征值分解操作，Ξ表示角度托普利兹矩阵；

(7d)按照下式，构造角度配对代价函数：

$f=\frac{1}{a^{H}a-a^H{\mathrm{UU}}^{H}a}$

其中，f表示角度配对代价函数，a表示立体搜索导向矢量，H表示共轭转置操作，U表示立体信号特征矩阵；

(7e)利用最大互斥角度配对方法，将步骤(5)得到的方位增广矩阵的方位角和步骤(6)得到的俯仰增广矩阵的俯仰角代入角度配对代价函数，得到二维角度配对的L型天线阵接收目标信号的方位角和L型天线阵接收目标信号俯仰角。

2.根据权利要求1所述的L型天线阵的低计算复杂度二维波达方向估计方法，其特征在于：步骤(7a)中所述的角度相关向量是按照如下公式构造的：

r＝[R(1,1),R(2,2),…,R(m,m),…,R(M,M)]

其中，r表示角度相关向量，R表示互相关矩阵，R(m,m)表示互相关矩阵R的第m行第m列元素，m＝1,2,…,M，M表示L型天线阵中x子阵的天线数。

3.根据权利要求1所述的L型天线阵的低计算复杂度二维波达方向估计方法，其特征在于：步骤(7b)中所述的角度托普利兹矩阵是按照如下公式构造的：

Ξ＝Toeplitz(r)

其中，Ξ表示角度托普利兹矩阵，Toeplitz(·)表示构造托普利兹矩阵操作，r表示角度相关向量。

4.根据权利要求1所述的L型天线阵的低计算复杂度二维波达方向估计方法，其特征在于：步骤(7d)中所述的立体搜索导向矢量a是按照如下公式构造的：

a＝[1,exp{j2πd(cosφ-cosθ)/λ},…,exp{j2π(M-1)d(cosφ-cosθ)/λ}]^T

其中，a表示立体搜索导向矢量，exp{·}表示指数函数操作，j表示虚数单位，d表示L型天线阵的x子阵中相邻天线的间距，φ表示任意的俯仰角，θ表示任意的方位角，λ表示目标信号波长，M表示L型天线阵中x子阵的天线数，T表示转置操作。