一种基于平行互质阵列时空扩展的二维波达方向估计方法与流程

本发明涉及阵列信号处理领域，尤其指一种基于平行互质阵列时空扩展的二维波达方向估计方法。

背景技术：

阵列信号处理技术在众多领域已得到广泛应用，而阵列信号处理的基本问题之一是空间信号波达方向估计(doa估计)。这一技术应用广泛，在通信、雷达等领域发挥着越来越重要的作用。在过去的几十年里人们不断提出了一系列高效的

doa估计算法。近几年来，随着压缩感知理论的不断成熟和广泛应用，基于稀疏重构技术的doa估计算法受到更多的关注和研究。使用均匀线阵估计信号的波达方向时，阵列的自由度直接受限于阵元个数，所能估计的信号源数目有限。由于非均匀阵列能够获得更大的阵列孔径和更高的自由度，并且阵元的摆放位置也具有更大的灵活性，所以研究非均匀阵列的设计和相应的doa估计算法具有很好的实际意义和应用价值。互质并行阵列是一种非均匀阵列，由两个平行放置的互质子阵列组成，阵元间距为半波长的整数倍。通过矢量化子阵的互协方差矩阵，得到阵元数目更多的虚拟阵列，在实际物理阵列数一定的情况下，增大阵列的孔径，使得所能估计的信号源数目大于实际物理阵元数目，并且doa估计的准确率也进一步提高。在文献anefficientdictionarylearning-based2-ddoaestimationwithoutpairmatchingforco-primeparallelarrays,ieeeaccess，2018(以下简称文献[1])，提出了一种利用时间和空间信息的互质并行阵列扩展的方法，利用自相关函数的共轭对称性，构造出更多的虚拟阵元，进一步拓展了阵列孔径。

上述文献[1]中的方法虽然对互质阵列进行了虚拟扩展，但仅仅利用了部分空间信息来构造虚拟阵列，如果所有的时间-空间信息都能被利用，可以得到更多的虚拟阵元数目，那么阵列孔径和自由度将会进一步增加，doa估计也会更加准确。

针对上述问题，我们提出了一个基于co-prime并行阵列时空扩展的二维doa估计方法。首先，我们在互质并行阵列的基础上增加了两个额外的阵元，利用信号自相关函数的共轭对称性，得到不同时滞的自相关函数，将该自相关函数作为新的接收数据，并求新的数据的协方差矩阵，最后通过稀疏重构的方法还原信号。所提出的方法可以利用空间和时间信息并形成共轭增强时空虚拟阵列，该方法能在物理阵元数目一定的情况下，利用空间和时间信息完全的扩展互质并行阵列，最大限度的提高阵列孔径和自由度以及得到更加准确的doa估计值。

技术实现要素：

本发明的目的在于克服现有方法存在的各种不足，首先提供一种基于平行互质阵列时空扩展的波达方向估计方法，在文献[1]原有互质阵列的基础上，增加两个物理阵元，充分利用接收信号的时间和空间信息对物理阵列进行虚拟扩展，形成共轭增强的时空虚拟阵列，最大限度的提高阵列孔径和自由度，得到更加准确的doa估计值。

另外本发明的另一目的在于提供一种基于互质阵列扩展的阵列结构。

为了实现上述目的，本发明可以通过以下技术方案实现。

一种基于平行互质阵列时空扩展的二维波达方向估计方法，其包括：选取和增加参考阵元，使互质阵列的子阵分别对参考阵元求互相关，得到不同时滞的相关函数，并将该自相关函数作为新的接收数据，拓展出更多的虚拟阵元，由新的接收数据伪采样得到伪数据矩阵，然后对伪接收数据的协方差矩阵矢量化，最后通过稀疏重构的方法还原信号。

进一步地，基于互质阵列扩展的阵列结构包括一个平行互质阵列以及两个单独的阵元，互质并行阵列包含两个水平均匀线阵，均匀线阵分别有m1、m2个阵元，并且m1和m2互质；线阵1的阵元间距为线阵2的阵元间距为两个子阵的间距为以坐标系来说明各个阵元的相对位置，以线阵1第一个参考阵元为原点，线阵1位于坐标系水平轴上，线阵2位于线阵1上方(不仅限于上方，此仅举例说明本阵列的结构)，间距为在互质并行阵列基础上增加的两个单独的阵元分别位于坐标点(λ,0),

进一步地，所述的一种基于平行互质阵列时空扩展的波达方向估计方法，其特征在于包含如下步骤：

步骤一：子阵1信号接收x1(t)＝a1s(t)+z1(t)

其中，a1＝a1(α)＝[a1(α1),a1(α2),…,a1(αk)]表示子阵1的流形矩阵，k表示信号源的个数，α＝[α1,α2,…,αk]，αk表示第k个信号与y轴之间的夹角，表示第k个信号的方向向量，s(t)＝[s1(t),s2(t),…,sk(t)]^t表示信号矢量，表示第k个信号源，z1(t)表示子阵1接收到的零均值高斯白噪声矢量，方差为

步骤二：子阵2信号接收x2(t)＝a2φs(t)+z2(t)，子阵2信号相对于子阵1信号有φ的相移，φ中包含俯仰角度信息；

其中，a2＝a2(α)＝[a2(α1),a2(α2),…,a2(αk)]表示子阵2的流形矩阵，k表示信号源的个数，α＝[α1,α2,…,αk]，αk表示第k个信号与y轴之间的夹角，表示第k个信号的方向向量，s(t)＝[s1(t),s2(t),…,sk(t)]^t表示信号矢量，表示第k个信号源，z2(t)表示子阵2接收到的零均值高斯白噪声矢量，方差为

β＝[β1,β2,…,βk]，βk表示第k个信号与x轴之间的夹角。

步骤三：新增参考阵元信号接收，βi表示信号源俯仰角，

其中，y(t)表示阵元接收到的k个信号源数据之和，si(t)表示第i个信号源数据，i表示第i个信号源，z3(t)表示阵元接收到的噪声数据。w(t)表示阵元(λ,0)接收到的k个信号源数据之和，z4(t)表示阵元(λ,0)接收到的噪声数据。

由步骤一和步骤二得到的x1(t)和x2(t)分别对原有参考阵元(0,0),以及新增加的参考阵元(λ,0)的数据求互相关，以及利用相关函数的共轭对称性得到如下数据：

其中m1、m2为子阵1和子阵2的阵元数目，下角标m1表示阵列的第m1个阵元，下角标m1+1表示阵列的第m1+1个阵元，以此类推。rxx(τ)表示子阵1和子阵2的数据分别对参考阵元(0,0),的数据作相关运算；rxy(τ)表示子阵1的数据对参考阵元的数据作相关运算；rxw(τ)表示子阵2的数据对参考阵元的数据作相关运算，字母τ表示相关函数的时延，求相关函数的公式定义为：

其中，表示第k个信号源的相关函数。(an,bn)为参考阵元的坐标，(am,bm)为子阵中的阵元坐标，m,n表示阵元的位置。

步骤四：构建虚拟阵列，由步骤三求得的数据，得到虚拟阵列接收信号：

对r1(τ)和r2(τ)以不同时延τ采样得到伪数据矩阵：

其中s^eqv＝[rs(ts),rs(2ts),…,rs(npts)]，np是虚拟矩阵接收数据的伪快拍数，ts是虚拟矩阵接收数据的伪采样周期。

求r1和r2的互协方差矩阵：

rc＝e[r2r1^h]

互协方差矩阵可以表示为：

其中，r^eqv＝e[s^eqv(s^eqv)^h]是一个对角阵，该矩阵第k个元素为且φr^eqv也是对角阵。

矢量化协方差矩阵，则有：

其中，⊙表示矩阵的khatri-rao积，向量u包含对角矩阵φr^eqv中的对角元素。由于r^eqv是实值对角矩阵，φ中对角元素包含角度β的数据，而向量u中的元素包含φ的信息，所以一旦求得向量u，就能计算出角度β。

稀疏重构算法求解。构造过完备字典θ来替代a，将空间划分为多个网格{θ1,θ2,…,θd}(d>>k)，真正的波达方向必定位于其中的网格附近。过完备字典由所有可能的方向{θ1,θ2,…,θd}(d>>k)构成。根据稀疏重构理论，将上述问题构建成一个稀疏优化问题，则有：

其中ρ＝[ρ1,ρ2,…,ρd]^t对应u中所有可能的方向。只要求解出ρ中的非零元素，就能求解出所有波达方向。具体来说，ρ中的非零元素所在的网格对应角度α，ρ中的非零元素的值对应于角度β，η是平衡稀疏度的正则化参数。

可以通过matlab软件的工具包求解得到ρ，ρ中非零元素的位置对应过完备基θ中网格的位置，此即为角度α的估计值。ρ中非零元素值为u中元素的估计值，用如下公式求解β：

βk＝arccos(-angle(uk)/π),k＝1,2,…,k,

在本式中，uk为ρ中第k个非零元素的值，βk表示信号第k信号的俯仰角，k表示信号源的个数。

与现有技术相比，本发明具有如下优点：本发明与现有互质并行阵列相比，充分的利用时间-空间信息扩展了虚拟阵列，在物理阵元数相同的条件下，能拓展出比原有互质并行阵列更多的虚拟阵元，具有更大的阵列孔径与自由度,均方误差更小，具有更高的实用性。

附图说明

图1是现有技术中互质并行阵列结构图。

图2是实施例中增加两个阵元的互质并行阵列结构图。

图3是实施例中互质并行阵列时空拓展后的阵元结构图。

图4是实施例中增加两个阵元的互质并行阵列时空拓展结构图。

图5是实施例中不同信噪比下各个算法的doa估计性能比较图。

图6是实施例中不同采样快拍数下各个算法的doa估计性能比较图。

图7是本发明基于平行互质阵列时空扩展的二维波达方向估计方法的流程图。

具体实施方式

下面结合实施例子及附图对本发明的实施作进一步详细的描述，但本发明的实施方式不限于此，需指出的是，以下若有未特别详细说明之过程或符号，均是本领域技术人员可参照现有技术理解或实现的。

1、本实施例的阵列结构

本发明提供的基于互质阵列扩展的阵列结构如图2所示，该阵列结构包括一个平行互质阵列以及两个单独的阵元。互质并行阵列包含两个水平均匀线阵，均匀线阵分别有m1、m2个阵元，且m1与m2互质。线阵1的阵元间距为λ表示信号的波长，线阵2的阵元间距为两个子阵的间距为以坐标系来说明各个阵元的相对位置，线阵1位于坐标系水平轴上，线阵2位于线阵1上方(不仅限于上方，此仅举例说明本阵列的结构)，间距为在互质并行阵列基础上增加的两个单独的阵元分别位于坐标点(λ,0),

2、基于平行互质阵列时空扩展的二维波达方向估计方法

本发明提供的基于互质并行阵列时空扩展的波达方向估计方法包含以下几个具体步骤：

步骤一：子阵1信号接收。k个窄带目标信号源分别为s1,s2,…,sk，信号源矢量记作s(t)＝[s1(t)s2(t)…sk(t)]^t。信号的波长为λ。这些目标声源对应于水平坐标轴的方向角分别为α＝[α1,α2,…,αk],对应于垂直坐标轴的方向角分别为β＝[β1,β2,…,βk]。αk表示第k个信号源的水平角，βk表示第k个信号源的俯仰角，则该阵元接收到的信号为：

其中，sk(t)表示第k个信号源，且有1≤k≤k。n1(t)表示第一个阵元上的噪声。

接收信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。那么第m个阵元在同一时刻接收到的信号为：

其中λ表示信号源波长，αk表示第k个信号源的水平角，nm(t)表示第m个阵元上的噪声。将各阵元的接收信号排列成列向量形式，则整个子阵1的接收信号可用以下矢量式子表示：

x1(t)＝a1s(t)+z1(t)(1)

其中，为m1×k的导

向矢量矩阵，为m1×1的接收信号矩阵，s(t)＝[s1(t),s2(t),…,sk(t)]^t为k×1的源信号矩阵，sk(t)表示第k个信号源，且有1≤k≤k，z1(t)为m1×1的噪声矩阵。

步骤二：子阵2信号接收，线阵2相对于线阵1有固定相位的延时，可以得到线阵2的信号接收模型为：

x2(t)＝a2φs(t)+z2(t)(2)

其中，为m2×k的导

向矢量矩阵，为m1×1的接收信号矩阵，s(t)＝[s1(t),s2(t),…,sk(t)]^t为k×1的源信号矩阵，sk(t)表示第k个信号源，且有1≤k≤k，z2(t)为m1×1的噪声矩阵，

步骤三：共轭增广时空处理，用坐标点(am,bm)，(an,bn)表示阵列中的两个不同的阵元，对应的接收信号表示为xm(t)和xn(t)，τ表示时延，其协方差矩阵表示为：

其中，表示为：

表示为：

因此，公式(3)可以表示为：

(1)对阵元(0,0)共轭增广

使(an,bn)＝(0,0)，则公式(3)可以表示为：

令对子阵1和子阵2分别有：

r⁽¹⁾(τ)＝a1rs(τ),r⁽²⁾(τ)＝a2φrs(τ)(5)

因为所以分别用和r^(1-)(τ)表示和r⁽¹⁾(-τ)^*的最后m1-1行，则有：

(2)对阵元共轭增广

使则公式(3)可以表示为：

令

则有：r^(1')(τ)＝a1φ^*rs(τ),r^(2')(τ)＝a2rs(τ)(8)

因为所以有

令r^(1'-)(τ)表示(r^(1')(-τ))^*的后面m1-1行，和r^(2'-)(τ)分别表示和(r^(2')(-τ))^*的后m2-1行。

(3)对阵元共轭增广，该阵元的接收信号为：

使则公式(3)可以表示为：

令则有：

r⁽³⁾(τ)＝a1φrs(τ)(11)

(4)对阵元(λ,0)共轭增广，该阵元的接收信号为：

使(an,bn)＝(λ,0)，则公式(3)可以表示为：

令则有：

因为所以令：

步骤四，构建虚拟阵列，由公式(5)(6)(8)(9)(11)(12)，则有：

r1(τ)和r2(τ)可以视为虚拟阵列的接收信号，由于虚拟阵列是通过不同时滞和空间位置的互相关获得的，因此虚拟阵列称为空间-时间虚拟阵列，如图4所示为时空拓展后的虚拟阵列结构。

如图3所示，在未添加两个单独阵元，仅对阵元(0,0)进行共轭增广的情况下，子阵1和子阵2分别扩展出两组虚拟阵元。如图4所示，在增加两个单独阵元，并对阵元(λ,0)，进行共轭增广的情况下，子阵1和子阵2分别再次拓展出一组阵元。比较图4和图3，与文献[1]中互质并行阵列的时空拓展相比，本发明方法拓展的虚拟阵列包含更多传感器元素和更大的孔径，可以增大自由度。

由r1(τ)和r2(τ)构造伪数据矩阵：

其中s^eqv＝[rs(ts),rs(ts),…,rs(npts)]，np是伪快拍数，ts是伪采样周期。

因此，求得r1和r2的互协方差矩阵：

互协方差矩阵可以表示为：

其中，r^eqv＝e[s^eqv(s^eqv)^h]，φr^eqv是对角阵。

矢量化协方差矩阵，则有：

其中，向量u包含对角矩阵φr^eqv中的对角元素。由于r^eqv是实值矩阵，φ中对角元素包含角度β，向量u中的元素包含φ，所以一旦求得u，就能计算出角度β。

步骤五，稀疏重构算法求解，在本步骤中，我们使用稀疏重构算法来求解u。首先构造过完备字典θ来替代a，将空间划分为多个网格{θ1,θ2,…,θd}(d>>k)，真正的波达方向必定位于其中的网格中，那么过完备字典由所有可能的方向{θ1,θ2,…,θd}(d>>k)构成。根据稀疏重构理论，将上述问题看成稀疏优化问题，则有：

其中ρ＝[ρ1,ρ2,…,ρd]^t对应u中所有可能的方向。只要求解出ρ中的非零元素，就能求解出所有波达方向。具体来说，ρ中的非零元素所在的网格对应角度α，ρ中的非零元素的值对应于角度β。η是平衡稀疏度的正则化参数。

公式(19)可以通cvx工具包求解ρ，ρ中非零元素的位置对应过完备基θ中网格的位置，此即为α的估计值。ρ中非零元素值为u的估计值，用如下公式求解β：

βk＝arccos(-angle(uk)/π),k＝1,2,…,k,(20)

uk为ρ中第k个非零元素的值。

以上算法流程图可以用图7表示。

本发明的主要工作步骤具体如下：

步骤1：信号源数目设定为k＝5，信号源与水平轴夹角为25°,36°,47°,58°,69°，与竖直轴水平夹角为83°,71°,59°,47°,35°，子阵1的阵元数目m1为4，子阵2的阵元数目m2为3。

步骤2：子阵1接收到的k个信号源的数据x1(t)＝a1s(t)+z1(t)，子阵2接收到的k个信号源数据为x2(t)＝a2φs(t)+n(t)，两个单独阵元接收到的k个信号源数据为以下又分为四个子步骤：

(1)子阵1对阵元(0,0)求自相关得到r⁽¹⁾(τ)，对r⁽¹⁾(τ)求共轭转置得到r^(1-)(τ)。子阵2对阵元(0,0)求自相关得到r⁽²⁾(τ)。

(2)子阵1对阵元求自相关并共轭转置得到r^(1'-)(τ)，子阵2对阵元求自相关得到r^(2')(τ)，对r^(2')(τ)求共轭转置得到r^(2'-)(τ)。

(3)子阵1对求自相关得到r⁽³⁾(τ)。

(4)子阵2对(λ,0)求自相关得到r^(4-)(τ)。

得到虚拟阵列的接收数据：

步骤3：对虚拟数据矩阵进行采样，得到r1＝[r1(ts),r1(ts),…,r1(npts)]，r2＝[r2(ts),r2(ts),…,r2(npts)]，其中伪快拍数设置为np＝500.ts＝1/4000s。求虚拟阵列接收数据的互协方差矢量化r＝vec(rc)。

步骤4：稀疏重构求解。构造过完备基θ，设置网格为[0:0.2:180]，则有

θ＝aa1*⊙aa1，求解稀疏优化问题可以用cvx工具包求解ρ，ρ中非零元素的位置对应过完备基θ中网格的位置，此即为α的估计值。ρ中非零元素值为u的估计值，用如下公式求解β：

βk＝arccos(-angle(uk)/π),k＝1,2,…,k,

在本式中，uk为ρ中第k个非零元素的值，βk表示信号第k信号的俯仰角，k表示信号源的个数。

步骤5，性能比较。实施例中用均方根误差(rmse)来评估各个算法的性能，其定义为

如图5所示，采样快拍数为500，参考信噪比snr从-10db到20db变化，间隔为5db，每个信噪比进行1000次蒙特卡洛实验。“methodin[1]m1＝4m2＝3”为文献[1]中的阵列结构下(如图1)，阵元数m1＝4,m2＝3时的均方误差，“methodin[1]m1＝4m2＝5”为文献[1]中的阵列结构下(如图1)，m1＝4，m2＝5时的均方误差，“ourmethodm1＝4m2＝3”为本发明所提供的阵列结构及扩展方法下(如图2)，增加两个单独阵元，且m1＝4，m2＝3时的均方误差。可以看出，本发明提供的算法性能优于文献[1]中算法，并且在阵元数目相等的情况下，性能同样优于原有算法。如图6所示，参考信噪比snr为10db，采样快拍数从100到1000变化，间隔为100，每个快拍数进行1000次蒙特卡洛实验，可以看出本发明所提供的算法同样性能较好。本发明在增加两个单独阵元的条件下，最大限度的对原有阵列进行了时间-空间拓展，得到更多的虚拟阵元，提高阵列孔径，具有更小的均方误差。