CN104730491B - 一种基于l型阵的虚拟阵列doa估计方法 - Google Patents

Abstract

本发明公开了一种基于L型阵列的虚拟阵列DOA估计方法，包括以下步骤：(1)基于移不变性质，将L型阵列的子阵Z_x,Z_y平移得到虚拟阵列Z_x',Z_y'，由于子阵的移不变性形成了两子信号的旋转不变性，虚拟子阵的信号等效于L型子阵Z_x,Z_y输入信号分别乘以旋转因子而得到；(2)将4个子阵的输出加以合并，构成虚拟阵列的输出信号Z(t)；(3)信号子空间和噪声子空间可用阵列输出的协方差矩阵的特征分解来描述，对阵列输出信号Z(t)进行互相关处理，得到R_zz，进行特征值分解得到信号子空间；(4)通过线性运算求解旋转因子，由其对角元素即可得到信号波达方向。本发明不需要计算谱函数，无需搜索峰值间接求解波达方向，降低了复杂度；减少了设备复杂性和成本；具有较高的定位精度。

Description

一种基于L型阵的虚拟阵列DOA估计方法

技术领域

本发明属于阵列信号处理二维DOA估计技术领域，尤其涉及的是一种基于L型阵列的虚拟阵元DOA估计方法。

背景技术

目标方位估计(DOA)估计是声呐、雷达、无线通信、医学成像、麦克风阵列处理等信号处理领域的一个重要分支，DOA解决的基本问题是确定同时处于某一空间领域中多个感兴趣的目标信号的空间位置，即各个目标信号到达传感器阵列的方向角。基于常规波束形成扫描的测向方法有固有的限定，受到“瑞利准则”的限制，估计的分辨率取决于阵列长度，仅当空间中两信源之间的分离角度大于阵列孔径的倒数时，才能被分辨。为了提高基阵的分辨率，当信号波长λ一定是，一般只有增加基阵的孔径长度，即增加阵元m或增大阵间距d，但是增加阵元个数会提高设备复杂性和成本，增大阵元间距又将引起次极大，出于实际情况的限制，阵列的尺度也不可能做得很大，因此仅仅依靠增加阵列孔径来达到提高分辨率的做法，实际工程中难以适用。为了克服这种局限性，设计了一种基于L型阵的虚拟阵列DOA估计方法，其基本思想是在有限尺度的阵列情况下，通过优化算法得到阵列的虚拟扩张，从而提高分辨率。

发明内容

本发明针对现有技术的不足，提出了一种基于L型阵的虚拟阵列DOA估计方法。

本发明的技术方案如下：

一种基于L型阵列的虚拟阵列DOA估计方法，其中，包括以下步骤：

步骤1：构造L型阵列，确定阵列接收的信号模型；

L型阵列中由2M-1个阵元接收声压时域信号，此L型阵列由x轴上阵元数为M的均匀线阵Z_x和y轴上阵元数为M的均匀线阵Z_y构成，其中，2M-1为L型阵列阵元数目，M为不小于2的整数，d为阵元间距。假设空间有K个信源入射到阵列上，其二维波达方向为 分别为第k个信号源的仰角和方位角。

假设入射到此阵列上的信号源数为K，则x轴、y轴上分别由M个阵元接收的信号分别为如下式(1)和式(2)：

x(t)＝A_xs(t)+n(t) (1)

y(t)＝A_ys(t)+n(t) (2)

式中s(t)为信号源矩阵，n(t)为噪声矩阵，A_x，A_y∈C^M×K，分别为L型阵列x轴、y轴上的方向矩阵，可表示为：

步骤2：借助旋转不变技术构造虚拟阵列获得输出信号矩阵Z；

借助旋转不变技术将L阵的子阵Z_x,Z_y进行虚拟扩张为子阵Z_x',Z_y'，由于子阵的移不变性形成了两个子阵信号的旋转不变性，即Z_x'的子阵信号为实际子阵Z_x的输入信号乘以旋转因子φ_x得到，Z_y'的子阵信号为实际子阵Z_y的输入信号乘以旋转因子φ_y得到，通过公式(1)与公式(2)先得到虚拟子阵的输出信号，然后将四个子阵输出加以合并，构成整个阵列的输出信号矩阵z(t)如下式(3)：

其中假设各信源的波达方向互不相同，则的列矢量之间线性独立，

并且

其中，矩阵φ_x，φ_y为K×K的对角矩阵，其对角元素为信号分别在Z_x，Z_y阵列上任意阵元之间的相位延迟，diag表示对角矩阵,即除了主对角线以外的元素均为零的方阵。

如式(3)z(t)包含x轴方向均匀线性子阵Z_x的输出信号x(t)_、y轴方向均匀线性子阵Z_y的输出信号y(t)、Z_x平移得到的虚拟子阵Z_x'的输出信号x'(t)、Z_y平移得到的虚拟子阵Z_y'的输出信号y'(t)，各阵列接收到的噪声相同，虚拟子阵Z_x'、Z_y'都为阵元数为M的均匀线阵。

步骤3：从阵列输出信号矩阵Z得到相关矩阵R_z，

信号子空间和噪声子空间可用阵列输出Z的协方差矩阵的特征分解得到，如式(4)所示：

式中R_s为信号的自相关矩阵，σ²为噪声方差，I为单位矩阵，式(4)中的E[.]，(.)^H分别表示为数学期望，共轭转置运算。

步骤4：将相关矩阵R_z做特征分解，估计信号个数；

阵列相关矩阵R_z可划分为两个空间，即K个的特征值对应的特征矢量E_s＝[s₁,s₂,...s_k]组成信号子空间，存在一个K×K的满秩矩阵T满足而且由于阵列的移不变特性E_s可分解为4部分，E_x,E_y,E_x',E_y'∈C^M×K，对应的子阵列分别为Z_x，Z_y，Z_x'，Z_y'，如式(5)所示，

步骤5：构造φ_x，φ_y的相似矩阵F，H；

由式(5)可推导出式(6)：

E_x'＝E_xT^-1φ_xT＝E_xT E_y'＝E_yT^-1φ_yT＝E_yH (6)

其中，F＝T^-1φ_xT,H＝T^-1φ_yT，T为满秩矩阵，因此F与φ_x，H与φ_y为相似矩阵，拥有相同的特征值，且其特征值为旋转因子φ_x，φ_y的对角元素。

步骤6：最小二乘法求解旋转因子φ_x，φ_y，计算波达方向

用最小二乘法解得旋转因子φ_x，φ_y如式(7)所示，便可从中得出信号的波达方向。

为E_x的伪逆，为E_y的伪逆，对F进行特征值分解得到 同时获得的估计值的u_k，对H进行特征值分解得到同时获得的估计值的v_k。可由式(8)估计出：

本发明具有如下优点：(1)与传统相干信号子空间算法相比，本发明不需要计算谱函数，无需搜索峰值间接求解波达方向，降低了复杂度；(2)在传感器数目确定的情况下，通过虚拟阵列增大阵列孔径，减少了设备复杂性和成本；(3)运用线性运算直接求解声源二维DOA的估计值，具有较高的定位精度。

附图说明

图1是L型阵列结构示意图；

图2是沿y轴方向平移的虚拟阵列Z_x'；

图3是沿x轴方向平移的虚拟阵列Z_y'。

图4是图4声源实际波达方向(/°)。

图5是三信源情况下MUSIC算法估计性能(SNR＝20dB)

图6是三信源情况下二维虚拟阵列算法估计性能(SNR＝20dB)

具体实施方式

下面结合附图和实施例对本发明作进一步的描述：

步骤1：构造L型阵列，确定阵列接收的信号模型；

如图1所示，L型阵列中由2M-1个阵元接收声压时域信号，此L型阵列由x轴上阵元数为M的均匀线阵Z_x和y轴上阵元数为M的均匀线阵Z_y构成，其中，2M-1为L型阵列阵元数目，M为不小于2的整数，d为阵元间距。假设空间有K个信源入射到阵列上，其二维波达方向为 分别为第k个信号源的仰角和方位角。

假设入射到此阵列上的信号源数为K，则x轴、y轴上分别由M个阵元接收的信号分别为如下式(1)和式(2)：

x(t)＝A_xs(t)+n(t) (1)

y(t)＝A_ys(t)+n(t) (2)

式中s(t)为信号源矩阵，n(t)为噪声矩阵，A_x，A_y∈C^M×K，分别为L型阵列x轴、y轴上的方向矩阵，可表示为：

步骤2：借助旋转不变技术构造虚拟阵列获得输出信号矩阵Z；

如图2、图3所示，借助旋转不变技术将L阵的子阵Z_x,Z_y进行虚拟扩张为子阵Z_x',Z_y'，由于子阵的移不变性形成了两个子阵信号的旋转不变性，即Z_x'的子阵信号为实际子阵Z_x的输入信号乘以旋转因子φ_x得到，Z_y'的子阵信号为实际子阵Z_y的输入信号乘以旋转因子φ_y得到，通过公式(1)与公式(2)先得到虚拟子阵的输出信号，然后将四个子阵输出加以合并，构成整个阵列的输出信号矩阵z(t)如下式(3)：

其假设各信源的波达方向互不相同，则的列矢量之间线性独立，

并且

其中，矩阵φ_x，φ_y为K×K的对角矩阵，其对角元素为信号分别在Z_x，Z_y阵列上任意阵元之间的相位延迟，diag表示对角矩阵,即除了主对角线以外的元素均为零的方阵。

如式(3)z(t)包含x轴方向均匀线性子阵Z_x的输出信号x(t)_、y轴方向均匀线性子阵Z_y的输出信号y(t)、Z_x平移得到的虚拟子阵Z_x'的输出信号x'(t)、Z_y平移得到的虚拟子阵Z_y'的输出信号y'(t)，各阵列接收到的噪声相同，虚拟子阵Z_x'、Z_y'都为阵元数为M的均匀线阵。

步骤3：从阵列输出信号矩阵Z得到相关矩阵R_z，

信号子空间和噪声子空间可用阵列输出Z的协方差矩阵的特征分解得到，如式(4)所示：

式中R_s为信号的自相关矩阵，σ²为噪声方差，I为单位矩阵，式(4)中的E[.]，(.)^H分别表示为数学期望，共轭转置运算。

步骤4：将相关矩阵R_z做特征分解，估计信号个数；

阵列相关矩阵R_z可划分为两个空间，即K个的特征值对应的特征矢量E_s＝[s₁,s₂,...s_k]组成信号子空间，存在一个K×K的满秩矩阵T满足而且由于阵列的移不变特性E_s可分解为4部分，E_x,E_y,E_x',E_y'∈C^M×K，对应的子阵列分别为Z_x，Z_y，Z_x'，Z_y'，如式(5)所示，

步骤5：构造φ_x，φ_y的相似矩阵F，H；

由式(5)可推导出式(6)：

E_x'＝E_xT^-1φ_xT＝E_xT E_y'＝E_yT^-1φ_yT＝E_yH (6)

其中，F＝T^-1φ_xT,H＝T^-1φ_yT，T为满秩矩阵，因此F与φ_x，H与φ_y为相似矩阵，拥有相同的特征值，且其特征值为旋转因子φ_x，φ_y的对角元素。步骤6：最小二乘法求解旋转因子φ_x，φ_y，计算波达方向

用最小二乘法解得旋转因子φ_x，φ_y如式(7)所示，便可从中得出信号的波达方向。

为E_x的伪逆，为E_y的伪逆，对F进行特征值分解得到 同时获得的估计值的u_k，对H进行特征值分解得到同时获得的估计值的v_k。可由式(8)估计出：

步骤7：基于L型阵的虚拟阵列DOA估计方法的运行仿真结果；

仿真条件为：采用的L型阵列纵轴、横轴都为8个全向阵元组成的均匀线阵，声音在空气中的传播速度为c＝340 m/s，声源频率为f＝3000 Hz，阵元间距取d＝λ/2，即d＝5.7cm，信噪比SNR＝20 dB。对比基于MUSIC算法与虚拟阵列算法的DOA估计性能，仿真模拟三个声源，实际波达方向为(10°,15°)、(30°,25°)、(50°,35°)，如图4所示，图中黑点即为声源真实的入射角度。运用MUSIC算法通过搜索峰峰值估计声源的波达方向，如图5所示，只能估计出1个声源的波达方向，DOA为(28°，24°)，旁瓣较多，定位精度较模糊，不能准确识别出声源的准确位置。通过声源位置与麦克风之间的关系，运用基于L型阵的虚拟阵列DOA估计方法直接计算出三个声源的波达估计值为(10.01°,15.07°)、(30°,25.01°)、(49.98°,34.99°)，如图6所示，图中黑点即为估计的DOA，与真实的波达方向相比，计算得到的DOA在±1°之间变化，误差较小，基于L型阵虚拟阵列算法的DOA估计方法定位精度相对较高。

结果中可以看出，与MUSIC算法相比，虚拟阵列算法不需要计算谱函数，无需搜索峰值间接求解波达方向，降低了复杂度；在传感器数目确定的情况下，通过虚拟阵列增大阵列孔径，减少了设备复杂性和成本；运用线性运算直接求解声源二维DOA的估计值，具有较高的定位精度。

Claims (1)

1.一种基于L型阵列的虚拟阵列DOA估计方法，其特征在于：该方法包括以下步骤：

步骤1：构造L型阵列，确定阵列接收的信号模型；

L型阵列中由2M-1个阵元接收声压时域信号，此L型阵列由x轴上阵元数为M的均匀线阵Z_x和y轴上阵元数为M的均匀线阵Z_y构成，其中，2M-1为L型阵列阵元数目，M为不小于2的整数，d为阵元间距；假设空间有K个信源入射到阵列上，其二维波达方向为θ_k,分别为第k个信号源的仰角和方位角；

假设入射到此阵列上的信号源数为K，则x轴、y轴上分别由M个阵元接收的信号分别为如下式(1)和式(2)：

x(t)＝A_xs(t)+n(t) (1)

y(t)＝A_ys(t)+n(t) (2)

式中s(t)为信号源矩阵，n(t)为噪声矩阵，A_x，A_y∈C^M×K，分别为L型阵列x轴、y轴上的方向矩阵，可表示为：

步骤2：借助旋转不变技术构造虚拟阵列获得输出信号矩阵Z；在传感器数目确定的情况下，通过虚拟阵列增大阵列孔径，减少了设备复杂性和成本；

借助旋转不变技术将L阵的子阵Z_x,Z_y进行虚拟扩张为子阵Z_x',Z_y'，由于子阵的移不变性形成了两个子阵信号的旋转不变性，即Z_x'的子阵信号为实际子阵Z_x的输入信号乘以旋转因子φ_x得到，Z_y'的子阵信号为实际子阵Z_y的输入信号乘以旋转因子φ_y得到，通过公式(1)与公式(2)先得到虚拟子阵的输出信号，然后将四个子阵输出加以合并，构成整个阵列的输出信号矩阵z(t)如下式(3)：

$z\left(t\right)={\[x\left(t\right),y\left(t\right),x^{\′}\left(t\right),y^{\′}\left(t\right)\]}^T=\stackrel{\‾}{A}s\left(t\right)+n\left(t\right)---\left(3\right)$

其中假设各信源的波达方向互不相同，则的列矢量之间线性独立，

并且

其中，矩阵φ_x，φ_y为K×K的对角矩阵，其对角元素为信号分别在Z_x，Z_y阵列上任意阵元之间的相位延迟，diag表示对角矩阵,即除了主对角线以外的元素均为零的方阵；

如式(3)z(t)包含x轴方向均匀线性子阵Z_x的输出信号x(t)、y轴方向均匀线性子阵Z_y的输出信号y(t)、Z_x平移得到的虚拟子阵Z_x'的输出信号x'(t)、Z_y平移得到的虚拟子阵Z_y'的输出信号y'(t)，各阵列接收到的噪声相同，虚拟子阵Z_x'、Z_y'都为阵元数为M的均匀线阵；

步骤3：从阵列输出信号矩阵Z得到相关矩阵R_z，

信号子空间和噪声子空间可用阵列输出Z的协方差矩阵的特征分解得到，如式(4)所示：

$R_z=E\[z\left(t\right)z^H\left(t\right)\]=\stackrel{\‾}{A}{R}_s{\stackrel{\‾}{A}}^H+{\mathrm{\σ}}^{2}I---\left(4\right)$

式中R_s为信号的自相关矩阵，σ²为噪声方差，I为单位矩阵，式(4)中的E[.]，(.)^H分别表示为数学期望，共轭转置运算；

步骤4：将相关矩阵R_z做特征分解，估计信号个数；

阵列相关矩阵R_z可划分为两个空间，即K个的特征值对应的特征矢量E_s＝[s₁,s₂,...s_k]组成信号子空间，存在一个K×K的满秩矩阵T满足而且由于阵列的移不变特性E_s可分解为4部分，E_x,E_y,E_x',E_y'∈C^M×K，对应的子阵列分别为Z_x，Z_y，Z_x'，Z_y'，如式(5)所示，

$E_s=\left[\begin{array}{c}{E}_x\\ E_y\\ E_{x^{\′}}\\ E_{y^{\′}}\end{array}\right]=\left[\begin{array}{c}{A}_{x}T\\ A_{y}T\\ A_x{\mathrm{\φ}}_{x}T\\ A_y{\mathrm{\φ}}_{y}T\end{array}\right]---\left(5\right)$

步骤5：构造φ_x，φ_y的相似矩阵F，H；

由式(5)可推导出式(6)：

E_x'＝E_xT^-1φ_xT＝E_xT E_y'＝E_yT^-1φ_yT＝E_yH (6)

其中，F＝T^-1φ_xT,H＝T^-1φ_yT，T为满秩矩阵，因此F与φ_x，H与φ_y为相似矩阵，拥有相同的特征值，且其特征值为旋转因子φ_x，φ_y的对角元素；

步骤6：最小二乘法求解旋转因子φ_x，φ_y，计算波达方向

用最小二乘法解得旋转因子φ_x，φ_y如式(7)所示，便可从中得出信号的波达方向；

$\begin{array}{cc}\hat{F}=E_x^{+}{E}_{x^{\′}}& \hat{H}=E_y^{+}{E}_{y^{\′}}\end{array}---\left(7\right)$

为E_x的伪逆，为E_y的伪逆，对F进行特征值分解得到同时获得sinθ_k的估计值的u_k，对H进行特征值分解得到同时获得sinθ_k的估计值的v_k；θ_k,可由式(8)估计出：