CN110749855B - 一种基于协方差域零化的均匀线阵波达方向估计方法 - Google Patents

Abstract

本发明提供了一种基于协方差域零化的均匀线阵波达方向估计方法，具体包括以下步骤：均匀线阵阵列数据接收；构造加权索引矩阵；构造协方差域的数据向量及对应的阵列残差的协方差矩阵；基于常规波束形成法的零化滤波器系数初估计；高斯牛顿算法迭代求解阵列残差和零化滤波器系数；根据零化滤波器系数计算所有信号的波达方向。该方法通过在协方差域进行零化，将阵列多快拍数据压缩为单快拍，相比于阵元域零化，无需进行奇异值分解，降低了计算量，且计算量不会随目标数增大而增大；本发明直接利用多块拍形成的阵列协方差矩阵，不需要对阵列数据进行主成分分析，可避免低信噪比或小快拍数条件下主成分分析性能差导致的波达方向估计精度低的问题。

Description

一种基于协方差域零化的均匀线阵波达方向估计方法

技术领域

本发明属于阵列信号处理领域，特别涉及对雷达、通信等电磁波和声学机械波信号波达方向的估计，具体是一种适用于均匀线阵的波达方向估计方法。

背景技术

波达方向估计属于阵列信号处理领域的一个重要分支，是用传感器阵列来估计信号源的来波到达方向，广泛应用于雷达、通信、声呐及医学诊断等军用和民用领域。为了对波达方向进行有效估计，研究者们提出了各种阵列构型及其对应波达方向估计算法。对于阵列构型，均匀线阵是阵元等间隔分布的直线阵，属于最基本的结构，是构成其他阵型的基础，有着广泛的应用。对于波达方向估计算法，人们希望其拥有高的估计精度的同时，能保持低的计算量，以便实时应用。

在时域中，零化滤波器理论可以用来对信号的频率进行高分辨估计，有着很好的性能。文献：Pan Y,Luo G Q,Jin H,et al.Direction-of-Arrival Estimation With ULA:A Spatial Annihilating Filter Reconstruction Perspective[J].IEEE Access.2018,6:23172-23179将时域的零化滤波器理论转化进了空域，实现了波达方向估计的功能，其估计性能优于传统的子空间类算法，且无需进行空域搜索。但将时域算法转至空域需要处理多块拍数据的利用问题。由于该算法是一种阵元域零化算法，所以其采取的是提取数据主成分的方法。具体来说是将阵列数据进行奇异值分解，取数量为目标数的最大奇异值对应的奇异值加权后的多个奇异向量作为新的阵列数据，此时快拍数压缩至目标数。但这样做会带来两个问题，一方面奇异值分解增大计算量且主成分分析会导致算法的计算量随目标数增加而增加，另一方面，由于主成分分析在低信噪比及小快拍数时性能较差，从而导致算法在对应场景下估计精度下降。

发明内容

针对现有方法的不足，本发明提出一种基于协方差域零化的均匀线阵波达方向估计方法。相对于阵元域零化算法，该方法可有效降低计算量，同时提高其在低信噪比及小快拍数场景下的估计精度。

本发明是通过以下技术方案实现的：一种基于协方差域零化的均匀线阵波达方向估计方法，所述方法包含以下步骤：

(S1)：K个远场窄带不相关信号入射到M阵元均匀线阵，K＜M，对各阵元输出进行匹配滤波得到N个快拍的阵列数据向量

n＝1,2,...,N。

(S2)：构造加权索引矩阵

(S3)：构造协方差域的数据向量y及对应的阵列残差的协方差矩阵Σ。

(S4)：对协方差域的数据向量y进行基于常规波束形成法的波达方向估计，并以此对零化滤波器系数

进行初估计。

(S5)：根据高斯牛顿算法迭代求解阵列残差ε和零化滤波器系数h。

(S6)：根据迭代求解得到的零化滤波器系数h计算所有信号的波达方向。

进一步地，步骤(S2)中构造加权索引矩阵

方法为：设置包含M个元素的列向量p，使p[m]＝m，其中p[m]为p的第m个元素，其余表示类似。令

其中(·)^T表示转置，vec(·)表示将矩阵的各列依次堆叠形成向量，e_M表示包含M个元素的全1列向量，其余表示类似。设置包含2M-1个元素的列向量p′，使p′[m]＝m-M，令

再将C中的非0元素改为0，0元素改为1,得到矩阵C′。最后将C′中的每行分别除以对应行的非零元素的个数得到加权索引矩阵B。

进一步地，步骤(S3)中构造协方差域的数据向量y及对应的阵列残差的协方差矩阵Σ的方法为：计算阵列协方差矩阵的估计值

并计算去噪后的向量化阵列协方差矩阵的估计值，得到

其中(·)^H表示共轭转置，I_M为有M个对角元的单位矩阵，

为阵列噪声方差

的估计值，通过取

的M-K个最小特征值的均值得到。则协方差域的数据向量

阵列残差的协方差矩阵

其中

表示Kronecker积。

进一步地，步骤(S4)中由基于常规波束形成法的波达方向估计结果，对零化滤波器系数

进行初估计的方法为：假设由常规波束形成法得到了K′个波达方向，设为

K′≤K。再定义函数

对该函数进行逆Z变换即得到零化滤波器系数h的初估计h⁽⁰⁾。其中

d为均匀线阵的阵元间距，λ为信号波长。

进一步地，步骤(S5)中迭代求解阵列残差ε和零化滤波器系数h的方法为：迭代之前需要进行初始化，阵列残差ε用零向量初始化，零化滤波器系数用步骤(S4)中的初估计h⁽⁰⁾初始化。迭代过程中ε与h的更新公式为：ε⁽ⁱ⁺¹⁾＝ε⁽ⁱ⁾+Δε，h⁽ⁱ⁺¹⁾＝h⁽ⁱ⁾+Δh，上标i表示迭代次数。其中Δh通过解方程组

得到，Δε通过

得到。其中，J₁＝-T(h⁽ⁱ⁾)，J₂＝L(y)-L(ε⁽ⁱ⁾)，T(·)与L(·)分别为两种Toeplitz算子。有

λ₂为辅助标量。迭代需设置最大迭代次数，每迭代一次检查是否达到最大迭代次数或者h是否收敛，若满足其中一个条件则停止迭代，若都不满足则继续迭代。

进一步地，步骤(S6)中根据迭代求解得到的零化滤波器系数h计算所有信号的波达方向的方法为：以h中元素为系数构造线性方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0，求其K个根，为

则K个信号的波达方向估计为

其中angle(·)表示求复数的幅角。

本发明相比现有技术的有益效果为：

首先，通过在协方差域进行零化，可将阵列多快拍数据压缩为单快拍，相比于阵元域零化，无需进行奇异值分解，降低了计算量，且计算量不会随信号源数增大而增大；其次，本发明直接利用多块拍形成的阵列协方差矩阵，不需要对阵列数据进行主成分分析，可避免低信噪比或小快拍数条件下主成分分析性能差导致的波达方向估计精度低的问题。

附图说明

图1为本发明的流程图；

图2为本发明方法与阵元域零化法在不同目标数下的计算耗时对比；

图3为本发明方法与阵元域零化法的估计均方根误差在不同信噪比下的对比；

图4为本发明方法与阵元域零化法的估计均方根误差在不同快拍数下的对比。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步详细说明。参照图1，本发明的具体实施步骤如下：

(S1)：K个远场窄带不相关信号入射到M阵元均匀线阵，K＜M，对各阵元输出进行匹配滤波得到N个快拍的阵列数据向量

n＝1,2,...,N。此时x(n)可建模为x(n)＝As(n)+υ(n) (1)

其中，s(n)为包含K个信号的列向量，K个信号互不相关；υ(n)为阵列噪声列向量，阵列噪声与来波信号互不相关，各阵元噪声满足独立复高斯分布，即

为噪声方差，I_M为有M个对角元素的单位阵；A＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流型，a(θ_k)为阵列流型向量，有

其中，(·)^T表示转置，d表示均匀线阵相邻阵元的间距，

(S2)：构造加权索引矩阵

设置包含M个元素的列向量p，使p[m]＝m，其中p[m]为p的第m个元素，其余表示类似。令

其中(·)^T表示转置，vec(·)表示将矩阵的各列依次堆叠形成向量，e_M表示包含M个元素的全1列向量，其余表示类似。设置包含2M-1个元素的列向量p′，使p′[m]＝m-M，令

再将C中的非0元素改为0，0元素改为1,得到矩阵C′。最后将C′中的每行分别除以对应行的非零元素的个数得到加权索引矩阵B。

(S3)：构造协方差域的数据向量y及对应的阵列残差的协方差矩阵Σ：计算阵列协方差矩阵的估计值

并计算去噪后的向量化阵列协方差矩阵的估计值，得到

其中(·)^H表示共轭转置，I_M为有M个对角元的单位矩阵，

为阵列噪声方差

的估计值，通过取

的M-K个最小特征值的均值得到。则协方差域的数据向量

阵列残差的协方差矩阵

其中

表示Kronecker积。

对于步骤(S2)和(S3)，我们一起做如下解释：

根据式(1)中的模型，我们可以得到去噪阵列协方差矩阵向量化后的真实值r为

其中，

(·)^*表示求共轭。

中的元素分别表示各个入射信号的功率。可见，由于快拍数有限，r与其估计值

存在误差Δr，即

这里将该误差称为阵列残差。这里阵列残差满足渐进高斯分布，即可近似表示为

我们知道，对于不相关信号入射到阵列上，阵列数据的协方差矩阵R的各个对角线上的元素相同，是一个Toeplitz矩阵。如果令矩阵

则R′同样为Toeplitz矩阵。于是我们将R′的各个对角线元素分别进行平均，形成一个列向量，R′右上角的元素为新的列向量的第一个元素，左下角的元素为新的列向量的最后一个元素。设该新向量为

则我们可以通过步骤(S2)构造B，使得

由于y是通过阵列协方差矩阵得到的，我们称其为协方差域的数据向量。这样我们用简单的矩阵相乘实现了上述平均Toeplitz矩阵对角元的操作。此处

为r的估计值，

为阵列噪声方差

的估计值，通过取

的M-K个最小特征值的均值得到。于是此时阵列残差相应会变为ε＝BΔr,根据式(5)，我们有

所以此时阵列残差的协方差矩阵为

(S4)：对协方差域的数据向量y进行基于常规波束形成法的波达方向估计，并以此对零化滤波器系数

进行初估计。假设由常规波束形成法得到了K′个波达方向，设为

k＝1,2,...,K′，K′≤K。再定义函数

对该函数进行逆Z变换即得到零化滤波器系数h的初估计h⁽⁰⁾。其中

d为均匀线阵的阵元间距，λ为信号波长。

这里由于我们设置

根据式(3)和(4)有y＝BA_vσ²+ε。根据A_v的表达式和B的构造方式，我们知道BA_vσ²的元素为等间隔采样的复指数的实数加权组合构成，即可认为BA_vσ²为一个新的均匀线阵的阵列输出数据。因此y-BΔr可被零化滤波器

零化。即

(y-ε)*h＝0 (6)

此处*表示卷积。对该零化滤波器进行z变换，有

由上式可知，若可以对θ_k,k＝1,2,...,K进行初估计，则可以用逆z变换对h进行初估计。这里我们取新的数据向量y，并采用计算量最小的常规波束形成波达方向估计法对θ_k,k＝1,2,...,K进行初估计。对于常规波束形成，由于其分辨率有限，如果存在两个目标的波达方向很接近，则常规波束形成的空间谱只能形成一个峰，所以会出现常规波束形成得到的波达方向个数K′≤K，但这不会影响算法最终的性能。

(S5)：根据高斯牛顿算法迭代求解阵列残差ε和零化滤波器系数h：迭代之前需要进行初始化，阵列残差ε用零向量初始化，零化滤波器系数用步骤(S4)中的初估计h⁽⁰⁾初始化。迭代过程中ε与h的更新公式为：ε⁽ⁱ⁺¹⁾＝ε⁽ⁱ⁾+Δε，h⁽ⁱ⁺¹⁾＝h⁽ⁱ⁾+Δh，上标i表示迭代次数。

其中Δh通过解方程组

得到，Δε通过

得到。其中，J₁＝-T(h⁽ⁱ⁾)，J₂＝L(y)-L(ε⁽ⁱ⁾)，T(·)与L(·)分别为两种Toeplitz算子。有

λ₂为辅助标量。迭代需设置最大迭代次数，每迭代一次检查是否达到最大迭代次数或者h是否收敛，若满足其中一个条件则停止迭代，若都不满足则继续迭代。

式(6)的零化关系可以写成矩阵相乘的形式L(y-ε)h＝0或T(h)(y-ε)＝0。于是为了解未知参数ε与h，我们建立如下优化问题：

上述优化问题可以采用高斯牛顿的思想转化为

式中的参数定义见步骤(S5)。解上式可用拉格朗日乘的方法，最终可得到步骤(S5)中的更新公式。

(S6)：根据迭代求解得到的零化滤波器系数h计算所有信号的波达方向：以h中元素为系数构造线性方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0，求其K个根，为

则K个信号的波达方向估计为

其中angle(·)表示求复数的幅角。

这一步中的求根操作是基于公式(7)的原理。这里最终的测向只需要对迭代求解后的零化系数h构成的方程求根，再加少量的运算即可完成。

为了验证本发明提出的基于协方差域零化的均匀线阵波达方向估计方法的正确性和相对于现有技术的优越性，做以下仿真实验。

考虑均匀现在M＝8，d＝λ/2的情形，于是此时共有8个阵元。对于本发明提出的方法的(S5)步骤中，最大迭代次数设为100，当||h⁽ⁱ⁺¹⁾-h⁽ⁱ⁾||₂/||h⁽ⁱ⁾||₂≤10^-8则认为h已收敛，其中上标i表示第i次迭代。采用的性能对比算法为背景技术中所引文献中的阵元域零化法。所有仿真结果都基于500次蒙特卡洛实验取平均。

实验一：

设置信噪比为20dB，快拍数为500个。目标数从1扫描到6，目标的波达方向集设为{-60°,-35°,-10°,10°,35°,60°}。当目标数为1时，取第一个波达方向，当目标数为2时，取前两个波达方向，以此类推。每运行一次算法，统计一下运行时间。计算机配置为：i7-7700CPU,16GB内存。阵元域零化算法和本发明算法的计算耗时对比如图2所示。可见，在1个目标的情况下，阵元域零化算法的计算耗时大于本发明算法，这是由于阵元域零化算法需要进行奇异值分解，且我们还发现阵元域零化算法迭代收敛算法稍慢于本发明算法。另外阵元域零化算法的计算耗时随着目标数增加而增加，几乎线性增长。这是由于阵元域零化算法的主成分分析提取的主成分数量等于目标数。本发明算法的计算耗时无上述现象，其计算量低于阵元域算法，大约只需1.5ms即可完成一次波达方向估计。

实验二：

设置快拍数为100，信噪比从-5dB扫描至15dB。设置两个目标，其来波方向分别为-4°和5°。运行算法，并统计波达方向估计的均方根误差。仿真结果如图3所示。可以发现，在高信噪比下阵元域零化算法的估计均方根误差与本发明所想趋同，但在低信噪比下，其均方根误差大于本发明算法。这是由于阵元域零化算法在低信噪比下主成分分析性能较差导致。

实验三：

设置快拍数为0dB，快拍数从20个扫描至200个，信号源设置同上。仿真结果如图4所示。可以发现，本发明的算法的均方根误差低于阵元域零化算法。特别是当快拍数小于60个时，即小快拍数下条件下，本发明算法的波达方向估计精度高出阵元域零化算法较多。

以上所述仅为本发名的较佳实施范例，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于协方差域零化的均匀线阵波达方向估计方法，其特征在于，包括以下步骤：

(S1)：K个远场窄带不相关信号入射到M阵元均匀线阵，K＜M，对各阵元输出进行匹配滤波得到N个快拍的阵列数据向量

(S2)：构造加权索引矩阵

所述步骤(S2)中构造加权索引矩阵

方法为：设置包含M个元素的列向量p，使p[m]＝m，其中p[m]为p的第m个元素，其余表示类似；令

其中(·)^T表示转置，vec(·)表示将矩阵的各列依次堆叠形成向量，e_M表示包含M个元素的全1列向量，其余表示类似；设置包含2M-1个元素的列向量p′，使p′[m]＝m-M，令

再将C中的非0元素改为0，0元素改为1,得到矩阵C′；最后将C′中的每行分别除以对应行的非零元素的个数得到加权索引矩阵B；

(S3)：构造协方差域的数据向量y及对应的阵列残差的协方差矩阵Σ；

所述步骤(S3)中构造协方差域的数据向量y及对应的阵列残差的协方差矩阵Σ的方法为：计算阵列协方差矩阵的估计值

并计算去噪后的向量化阵列协方差矩阵的估计值，得到

其中(·)^H表示共轭转置，I_M为有M个对角元的单位矩阵，

为阵列噪声方差

的估计值，通过取

的M-K个最小特征值的均值得到；则协方差域的数据向量

阵列残差的协方差矩阵

其中

表示Kronecker积；

(S4)：对协方差域的数据向量y进行基于常规波束形成法的波达方向估计，并以此对零化滤波器系数

进行初估计；

所述步骤(S4)中由基于常规波束形成法的波达方向估计结果，对零化滤波器系数

进行初估计的方法为：假设由常规波束形成法得到了K′个波达方向，设为

再定义函数

对该函数进行逆Z变换即得到零化滤波器系数h的初估计h⁽⁰⁾；其中

d为均匀线阵的阵元间距，λ为信号波长；

(S5)：根据高斯牛顿算法迭代求解阵列残差ε和零化滤波器系数h；

所述步骤(S5)中迭代求解阵列残差ε和零化滤波器系数h的方法为：迭代之前需要进行初始化，阵列残差ε用零向量初始化，零化滤波器系数用步骤(S4)中的初估计h⁽⁰⁾初始化；迭代过程中ε与h的更新公式为：ε⁽ⁱ⁺¹⁾＝ε⁽ⁱ⁾+Δε，h⁽ⁱ⁺¹⁾＝h⁽ⁱ⁾+Δh，上标i表示迭代次数；其中Δh通过解方程组

得到，Δε通过

得到；其中，J₁＝-T(h⁽ⁱ⁾)，J₂＝L(y)-L(ε⁽ⁱ⁾)，T(·)与L(·)分别为两种Toeplitz算子；有

λ₂为辅助标量；迭代需设置最大迭代次数，每迭代一次检查是否达到最大迭代次数或者h是否收敛，若满足其中一个条件则停止迭代，若都不满足则继续迭代；

(S6)：根据迭代求解得到的零化滤波器系数h计算所有信号的波达方向；

所述步骤(S6)中根据迭代求解得到的零化滤波器系数h计算所有信号的波达方向的方法为：以h中元素为系数构造线性方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0，求其K个根，为

则K个信号的波达方向估计为

其中angle(·)表示求复数的幅角。

Images