CN113325365B - 一种基于四元数的相干信号二维doa估计方法 - Google Patents

Abstract

本发明提供一种基于四元数的相干信号二维DOA估计方法，包括：建立观测信号复数模型，建立观测信号四元数模型，构造相关矩阵，四元数互相关向量化简，通过相关矩阵平滑实现信源解相干，利用传播算子方法对信号俯仰角进行估计，估计传播算子，估计入射信号水平方位角，重构阵列响应矩阵。本发明解决“如何在充分利用阵列数据的同时又能消除所引入的额外加性噪声”这一问题。

Description

一种基于四元数的相干信号二维DOA估计方法

技术领域

本发明涉及利用矢量传感器阵列对多个相干信号源的波达方向进行精确估计，尤其涉及一种基于四元数的相干信号二维DOA估计方法，属于阵列信号处理领域。

背景技术

相比传统标量传感器，声矢量传感器在测量声场中某一点声压的同时还可提供该点的三维粒子振速信息，从而可实现对声场信息更为全面、充分的获取，因此近些年来已在水声阵列信号处理领域扮演起重要角色。传统方法常将矢量水听器的输出信号在复数域内建模，并将信号各个分量拼接成“长矢量”形式，这种信号建模方式不但破环了矢量传感器输出信号各分量间的正交关系，也没有完整利用阵列的结构信息。因此近些年来，许多学者开始尝试用四元数对矢量信号进行建模，并已取得了众多成就。

对多个信号源的波达方向进行精确估计(DOA)是阵列信号处理领域所研究的主要内容之一，声纳、雷达、地震探测和通信***等领域都有着相关技术的身影，而此类技术的代表之一便是子空间方法。传统子空间方法都是通过对数据协方差矩阵进行特征分解来获得信号子空间或噪声子空间，再根据噪声子空间与信号子空间的正交性(MUSIC)或是在一定的阵列布放形式下信号子空间所具的旋转不变性(ESPRIT)来获得信号的DOA估计，其中，高维数据协方差矩阵的特征分解往往是费时费力的，为解决这一问题，传播算子方法应运而生。该方法首先对协方差矩阵进行子阵划分，随后通过线性操作便可获得噪声子空间或信号子空间，之后再以类似MUSIC或ESPRIT算法的方式获得信号的波达角，但由于传播算子方法未充分利用噪声的空域特性，导致低信噪比下算法的性能并不理想。因此，许多基于阵列互相关并应用传播算子技术的算法(如SUMWE、CODE)已被提出，上述算法通过阵列互相关成功抑制了加性噪声，从而为传播算子方法的优良性能提供了保证。

但通过阵列互相关去除加性噪声的缺点在于：(1)对阵列进行子阵划分这一步骤损失了阵列孔径；(2)只考虑各子阵之间的互相关而忽略子阵内部的自相关信息无法实现对阵列数据的高效利用。

因此，如何在充分利用阵列数据(即：既考虑阵列的互相关又能充分利用阵元的自相关信息)的同时又能消除所引入的额外加性噪声(由阵元自相关所引入)便成为了使此类算法性能获得进一步增益的关键突破口。

发明内容

为克服现有技术的不足，本发明提出一种在四元数域内对信号进行建模的全新思路，旨在解决“如何在充分利用阵列数据的同时又能消除所引入的额外加性噪声”这一问题。应当指出，所提方法具有一定的普适性，本专利在将所提方法与SUMWE算法相结合后(AQ-SUMWE)所获得的优良效果则进一步证明了该种方法的可行性、有效性。

本发明的目的是这样实现的：步骤如下：

步骤一：建立观测信号复数模型；

步骤二：建立观测信号四元数模型；构造如下两个四元数信号模型：

模型一：

z(t)＝p(t)+x(t)j+y(t)j

＝A(s(t)+Γ_xs(t)j+Γ_ys(t)j)+n(t)

模型二：

步骤三：构造相关矩阵：考虑第m个阵元输出信号的四元数模型与整个阵列输出信号的四元数模型之间的互相关，有：

将r_m展开，有如下形式：

步骤四：四元数互相关向量化简，r_m的最终形式为：

其中，

R_s＝E{s(t)s^H(t)}为信号协方差矩阵；

步骤五：通过相关矩阵平滑实现信源解相干；

步骤六：利用传播算子方法对信号俯仰角进行估计：

步骤七：估计传播算子，将R进行如下划分：

其中R₁，R₂分别是R的前K行和后M-2K+1行，R₂＝P₀R₁

因此，传播算子可以表示为：

根据样本数据估计传播因子

俯仰角

可以从

的K个特征值

中获得：

步骤八：估计入射信号水平方位角，

将中给出的x(t)和y(t)级联，形成如下阵列输出向量z′(t)：

引入如下2M×M互相关矩阵R′：

其中

和n_p(t)在空间上彼此独立，即

因此在统计意义上，R′中的噪声已被消除；用R′_(1:M,:)、R′_(M+1:2M,:)分别表示R′的前M行和后M行，将

代入，有：

w_x,k,w_y,k,k＝1,2,...,K分别为W_x,W_y的第k行；由关系：u_k tanθ_k＝v_k，可知

tanθ_kw_x,k＝w_y,k

步骤九：重构阵列响应矩阵，参考A(φ)＝[a(φ₁),a(φ₂),...,a(φ_K)]并使用

k＝1,2,…,K来重构矩阵A，即

W_x,W_y的估计值为：

得第k个信号的方位角可以估计为：

至此，得到了入射信号的俯仰角和方位角，且两者之间自动匹配。

本发明还包括这样一些结构特征：

1.步骤一具体为：假设阵列为沿z轴按半波长等间距布放的矢量传感器线列阵，单矢量传感器包含一个声压传感器和两个分别指向x轴、y轴的质点振速传感器，则矢量传感器流形矢量可写为：

其中，θ_k和φ_k分别为第k个信号的水平角和俯仰角，u_k和v_k为方向余弦，分别对应于指向x轴和y轴的质点振速传感器；假设K个远场窄带信号入射到传感器阵列上，则可分别在复数域内建立声压传感器阵列、x轴指向的质点振速传感器阵列、y轴指向的质点振速传感器阵列的接收信号模型，依次为：

其中，A(φ)＝[a(φ₁),a(φ₂),...,a(φ_K)]为声压传感器阵列流形，

为对应第k个信号的导向矢量，

为传感器间空域相位因子，d为阵元间距，λ为信号波长，信号矢量s(t)＝[s₁(t),...,s_K(t)]^T为复高斯随机过程，β＝[β₁,…,β_K]^T为复衰减系数，其中β₁＝1且β_k≠0,k＝2,3,…K，n_l,l＝p,x,y为与信号不相关且彼此独立的噪声矢量，Γ_x＝diag{cosθ₁sinφ₁,…,cosθ_K sinφ_K}，Γ_y＝diag{sinθ₁sinφ₁,...,sinθ_K sinφ_K}。

2.步骤五具体为：考虑K个相干信号，我们通过对相关矩阵r_m进行平滑来构造矩阵

可以表示为：

其中，A₀包含了A的前M-K+1行，

是选择矩阵，

和

是两个零矩阵，Π_m＝[ψ_m,Φ¹ψ_m,...,Φ^K-1ψ_m],

可以发现，

的秩为K，因此

可以用来估计K个相干信号的方向；

定义如下的相关矩阵R以进行后续处理:

其中，Π＝[Π₁,Π₂,...,Π_M]。

3.步骤六具体为：传播算子方法基于对矩阵A₀的子阵划分，即

其中，因为A₁是一个K行的满秩矩阵，A₂为M-2K+1行的矩阵；A₀的前K行是线性无关的，则A₀的其余行可以表示为A₀前K行的线性组合；传播矩阵P₀定义为一个与矩阵A₁、A₂有关的线性算子：

P₀A₀＝A₁

则得：

其中，

用P_a,P_b分别表示P的前M-K行和后M-K行，用A_a,A_b表示A₀的前M-K行和后M-K行；有：

其中，P_a,P_b,Φ的关系为：

其中，

Φ＝diag{ξ₁,ξ₂,...,ξ_K}，Φ中的第k个对角线元素对应于第k个信号的传感器空间相位因子；

和Φ是相似矩阵，对

的特征分解得到的特征值等于对角元素ξ_k,k＝1,2,...,K；俯仰角估计结果φ_k从下式获得：

与现有技术相比，本发明的有益效果是：1.在四元数域内对矢量传感器输出信号进行建模保留了信号的矢量特性，也充分利用了信号各分量间的正交关系。相比复数域内的长矢量建模方式，四元数模型为多分量复信号提供了一种更为紧凑、方便的处理方式。2.可充分利用声矢量传感器阵列中相邻传感器间的空域相位信息与振速传感器内固有的方向余弦信息，对信息更为充分的利用实现了算法性能的增益。3.本方法不但考虑了阵元间的互相关信息，也考虑了自相关信息。传统方法中，由自相关引入的加性噪声无法去除，这将严重限制传播算子方法的性能，而本方法中，通过对数据在四元数域内进行处理，这种加性噪声已被统计消除，从而使传播算子方法得到了更为精确的DOA估计结果。

附图说明

图1为本发明流程框图；

图2为阵列布放示意图；

图3为不同信噪比下算法性能对比(图中，AQ-SUMWE为本专利所提方法，V-SUMWE与FBSS-ESPRIT为现有算法，下同)；

图4为不同快拍数下算法性能对比；

图5为不同阵元数下算法性能对比。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述。

结合图1至图5，本发明的步骤如下：

步骤一：建立观测信号复数模型：

假设阵列为沿z轴按半波长等间距布放的矢量传感器线列阵，单矢量传感器包含一个声压传感器和两个分别指向x轴、y轴的质点振速传感器，则矢量传感器流形矢量可写为：

其中，θ_k和φ_k分别为第k个信号的水平角和俯仰角，u_k和v_k为方向余弦，分别对应于指向x轴和y轴的质点振速传感器。假设K个远场窄带信号入射到传感器阵列上，则可分别在复数域内建立声压传感器阵列、x轴指向的质点振速传感器阵列、y轴指向的质点振速传感器阵列的接收信号模型，依次为：

其中，A(φ)＝[a(φ₁),a(φ₂),...,a(φ_K)]为声压传感器阵列流形，

为对应第k个信号的导向矢量，

为传感器间空域相位因子，d为阵元间距，λ为信号波长，信号矢量s(t)＝[s₁(t),...,s_K(t)]^T为复高斯随机过程，β＝[β₁,...,β_K]^T为复衰减系数，其中β₁＝1且β_k≠0,k＝2,3,…K，n_l,l＝p,x,y为与信号不相关且彼此独立的噪声矢量，Γ_x＝diag{cosθ₁sinφ₁,...,cosθ_K sinφ_K}，Γ_y＝diag{sinθ₁sinφ₁,...,sinθ_K sinφ_K}。

步骤二：建立观测信号四元数模型：

现根据以上信号模型构造如下两个四元数信号模型：

模型一：

模型二：

步骤三：构造相关矩阵：

考虑第m个阵元输出信号的四元数模型与整个阵列输出信号的四元数模型之间的互相关，有

将r_m展开，有如下形式：

步骤四：四元数互相关向量化简，

根据四元数运算规则：对于

满足jx＝x^*j，再考虑到信号为圆信号，我们有

E{s(t)js^H(t)}＝E{βs₁(t)s₁(t)β^Tj}＝0 (4.1)

所以，r_m可化简为如下形式：

其中，

当l≠m时，由于各噪声之间独立，故有

当l＝m且噪声为圆信号时，有

在上面的推导中，我们使用了以下两个假设：1)噪声是空间独立的，即

2)噪声向量n_p(t),n_x(t),n_y(t)皆为圆信号，即E{n_l,m(t)n_l,m(t)}＝0,l＝p,x,y。通过这种方法，我们成功地利用四元数代数理论在统计水平上消除了r_m中的加性噪声。

根据以上推导，r_m的最终形式为：

其中，

R_s＝E{s(t)s^H(t)}为信号协方差矩阵。

步骤五：通过相关矩阵平滑实现信源解相干：

考虑K个相干信号，我们通过对相关矩阵r_m进行平滑来构造矩阵

可以表示为

其中，A₀包含了A的前M-K+1行，

是选择矩阵，

和

是两个零矩阵，

可以发现，

的秩为K，因此

可以用来估计K个相干信号的方向。

最后，我们定义如下的相关矩阵R以进行后续处理:

其中，Π＝[Π₁,Π₂,...,Π_M]。

步骤六：利用传播算子方法对信号俯仰角进行估计：

传播算子方法基于对矩阵A₀的子阵划分，即

其中，因为A₁具有Vandermonde结构，所以是一个K行的满秩(非奇异)矩阵，A₂为M-2K+1行的矩阵。A₀的前K行是线性无关的，则A₀的其余行可以表示为A₀前K行的线性组合。传播矩阵P₀可以唯一定义为一个与矩阵A₁、A₂有关的线性算子:

P₀A₀＝A₁ (6.2)

由，得

其中，

用P_a,P_b分别表示P的前M-K行和后M-K行，同样，用A_a,A_b表示A₀的前M-K行和后M-K行。我们有

其中，P_a,P_b,Φ的关系为

这里，

Φ＝diag{ξ₁,ξ₂,...,ξ_K}，Φ中的第k个对角线元素对应于第k个信号的传感器空间相位因子。根据上式可知，

和Φ是相似矩阵，对

的特征分解得到的特征值等于对角元素ξ_k,k＝1,2,...,K。因此，俯仰角估计结果φ_k可以从下式获得

步骤七：估计传播算子：

为了估计传播算子，我们将R进行如下划分

其中R₁，R₂分别是R的前K行和后M-2K+1行，由有

R₂＝P₀R₁ (7.2)

因此，传播算子可以表示为

根据样本数据估计传播因子

俯仰角

可以从

的K个特征值

中获得：

步骤八：估计入射信号水平方位角，

这里，我们通过使用嵌入在声学粒子速度分量中的方向信息来获得方位角估计。首先，通过将中给出的x(t)和y(t)级联，可以形成如下阵列输出向量z′(t)：

然后，引入如下2M×M互相关矩阵R′：

注意，

和n_p(t)在空间上彼此独立，即

因此在统计意义上，R′中的噪声已被消除。我们用R′_(1:M,:)、R′_(M+1:2M,:)分别表示R′的前M行和后M行，将

代入，有

w_x,k,w_y,k,k＝1,2,...,K分别为W_x,W_y的第k行。由关系：u_k tanθ_k＝v_k，可知

tanθ_kw_x,k＝w_y,k (8.4)

步骤九：重构阵列响应矩阵，

为了估计阵列响应矩阵A，我们可以参考A(φ)＝[a(φ₁),a(φ₂),…,a(φ_K)]并使用

k＝1,2,…,K来重构矩阵A，即

因此，W_x,W_y的估计值为

然后，根据和，可得第k个信号的方位角可以估计为

至此，我们便获得了入射信号的俯仰角和方位角，且两者之间自动匹配，无需额外操作。图3、图4、图5中，我们分别对比了所提算法与已有算法(V-SUMWE、FBSS-ESPRIT)在不同参数(信噪比、快拍数、阵元数)下的性能优劣，仿真结果表明，所提出的算法均取得了更为精确的估计结果。

算法具体实施步骤见图1，应当指出，本专利所提出的“利用四元数消除加性噪声的方法”具有一定的普适性，可视为对数据的一种预处理，因此该方法还可与多种算法相结合(本专利以SUMWE为例)以使算法获得性能上的进一步提升。在不脱离本发明技术原理的前提下，本领域的技术人员还可对该方法做出若干变形(例如将该方法与其他已有算法相结合)，这些变形也应在本发明的保护范围之内。

Claims (3)

1.一种基于四元数的相干信号二维DOA估计方法，其特征在于：步骤如下：

步骤一：建立观测信号复数模型；

假设阵列为沿z轴按半波长等间距布放的矢量传感器线列阵，单矢量传感器包含一个声压传感器和两个分别指向x轴、y轴的质点振速传感器，则矢量传感器流形矢量可写为：

其中，θ_k和φ_k分别为第k个信号的水平角和俯仰角，u_k和v_k为方向余弦，分别对应于指向x轴和y轴的质点振速传感器；假设K个远场窄带信号入射到传感器阵列上，则可分别在复数域内建立声压传感器阵列、x轴指向的质点振速传感器阵列、y轴指向的质点振速传感器阵列的接收信号模型，依次为：

其中，A＝[a(φ₁),a(φ₂),...,a(φ_K)]为声压传感器阵列流形，

为对应第k个信号的导向矢量，

为传感器间空域相位因子，d为阵元间距，λ为信号波长，信号矢量s(t)＝[s₁(t),...,s_K(t)]^T为复高斯随机过程，β＝[β₁,...,β_K]^T为复衰减系数，其中β₁＝1且β_k≠0,k＝2,3,…K，n_l,l＝p,x,y为与信号不相关且彼此独立的噪声矢量，Γ_x＝diag{cosθ₁ sinφ₁,...,cosθ_K sinφ_K}，Γ_y＝diag{sinθ₁ sinφ₁,...,sinθ_K sinφ_K}；

步骤二：建立观测信号四元数模型；构造如下两个四元数信号模型：

模型一：

z(t)＝p(t)+x(t)j+y(t)j

＝A(s(t)+Γ_xs(t)j+Γ_ys(t)j)+n(t)

模型二：

步骤三：构造相关矩阵：考虑第m个阵元输出信号的四元数模型与整个阵列输出信号的四元数模型之间的互相关，有：

将r_m展开，有如下形式：

步骤四：四元数互相关向量化简，r_m的最终形式为：

其中，

R_s＝E{s(t)s^H(t)}为信号协方差矩阵；

步骤五：通过相关矩阵平滑实现信源解相干；

步骤六：利用传播算子方法对信号俯仰角进行估计：

步骤七：估计传播算子，将R进行如下划分：

其中R₁，R₂分别是R的前K行和后M-2K+1行，R₂＝P₀R₁

因此，传播算子表示为：

根据样本数据估计传播因子

俯仰角

从

的K个特征值

中获得：

步骤八：估计入射信号水平方位角，

将步骤一中给出的x(t)和y(t)级联，形成如下阵列输出向量z′(t)：

引入如下2M×M互相关矩阵R′：

其中

和n_p(t)在空间上彼此独立，即

因此在统计意义上，R′中的噪声已被消除；用R′_(1:M,:)、R′_(M+1:2M,:)分别表示R′的前M行和后M行，将

代入，有：

w_x,k,w_y,k,k＝1,2,...,K分别为W_x,W_y的第k行；由关系：u_k tanθ_k＝v_k，知

tanθ_kw_x,k＝w_y,k

步骤九：重构阵列响应矩阵，参考A＝[a(φ₁),a(φ₂),...,a(φ_K)]并使用

来重构矩阵A，即

W_x,W_y的估计值为：

得第k个信号的方位角估计为：

至此，得到了入射信号的俯仰角和方位角，且两者之间自动匹配。

2.根据权利要求1所述的一种基于四元数的相干信号二维DOA估计方法，其特征在于：步骤五具体为：考虑K个相干信号，我们通过对相关矩阵r_m进行平滑来构造矩阵

可以表示为：

其中，A₀包含了A的前M-K+1行，

是选择矩阵，

和

是两个零矩阵，Π_m＝[ψ_m,Φ¹ψ_m,...,Φ^K-1ψ_m],

可以发现，

的秩为K，因此

可以用来估计K个相干信号的方向；

定义如下的相关矩阵R以进行后续处理:

其中，Π＝[Π₁,Π₂,...,Π_M]。

3.根据权利要求2所述的一种基于四元数的相干信号二维DOA估计方法，其特征在于：步骤六具体为：传播算子方法基于对矩阵A₀的子阵划分，即

其中，因为A₁是一个K行的满秩矩阵，A₂为M-2K+1行的矩阵；A₀的前K行是线性无关的，则A₀的其余行可以表示为A₀前K行的线性组合；传播矩阵P₀定义为一个与矩阵A₁、A₂有关的线性算子：

P₀A₀＝A₁

则得：

其中，

用P_a,P_b分别表示P的前M-K行和后M-K行，用A_a,A_b表示A₀的前M-K行和后M-K行；有：

其中，P_a,P_b,Φ的关系为：

其中，

Φ＝diag{ξ₁,ξ₂,...,ξ_K}，Φ中的第k个对角线元素对应于第k个信号的传感器空间相位因子；

和Φ是相似矩阵，对

的特征分解得到的特征值等于对角元素ξ_k,k＝1,2,...,K；俯仰角估计结果φ_k从下式获得：

Images