CN113219399A - 基于全实值计算的远场窄带无线电信号波达方向估计方法 - Google Patents

Abstract

本发明涉及信号处理技术领域，具体地说是一种能够显著降低计算复杂度、提高精度的基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，获取辐射源信号后，计算信号协方差矩阵，构造实导向矢量后，构造实值多项式，计算多项式系数及根，解模糊并获得信号波达方向，本发明与现有技术相比有益效果为：实现了特征值分解和多项式求根的全实值计算，避免了进行复值多项式求根所需的大量计算，为波达方向的工程化实现提供了技术支持。

Description

基于全实值计算的远场窄带无线电信号波达方向估计方法

技术领域：

本发明涉及信号处理技术领域，具体地说是一种能够显著降低计算复杂度、 提高精度的基于全实值计算的远场窄带无线电信号波达方向估计方法。

背景技术：

信号的波达方向估计是雷达、声呐、无线通信和无源定位等应用中经常遇到 的重要研究课题，以多重信号分类和旋转不变子空间为代表的子空间类算法的提 出，实现了传统空间谱估计向超分辨测角的飞跃，但MUSIC算法庞大的计算量和 ESPRIT算法较低的估计精度阻碍了超分辨算法的工程化进度。针对此问题，促 成了酉MUSIC(Unitary MUSIC,U-MUSIC)算法的诞生。

酉求根MUSIC算法是最经典的实值波达方向(Direction of Arrival,DOA) 估计方法之一。与传统的求根MUSIC算法在特征值分解(Eigenvalue Decomposition,EVD)与多项式求根两个过程都需要进行复值运算不同的是，酉 求根算法利用酉变换、前向/后向(Forward/Backward,FB)平滑和坐标变换等 数学技术，实现了在对称实矩阵上进行实值特征值分解。

但是，与其它经典的酉算法相比，例如，酉MUSIC(U-MUSIC)，酉旋转不变 子空间法，酉求根MUSIC算法只能被视为半实值的估计方法，因为它只能在EVD 阶段实现有限的实值计算，而在多项式的计算上仍需要进行大量的复值计算。

发明内容：

本发明针对经典的酉求根MUSIC算法只能在EVD过程实现实值计算而在多 项式求根中仍为复值计算的问题，提出一种通过三角多倍角公式改写导向矢量， 以及变量替换实现实值多项式的构建，实现了在多项式求根过程的全实值计算， 显著降低了计算复杂度，为波达方向实际工程化提供重要技术支持的基于全实值 计算的远场窄带无线电信号波达方向估计方法。

本发明通过以下措施达到：

一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于， 获取辐射源信号后，计算信号协方差矩阵，构造实导向矢量后，构造实值多项式， 计算多项式系数及根，解模糊并获得信号波达方向。

本发明所述获取辐射源信号具体包括以下内容：

第一步，利用天线阵列接收辐射源信号，天线阵列具有L个相互独立的阵元， 以d等间距组成均匀线阵(ULA)，空间中存在K个远场窄带信号入射到阵列，其 中，假设K先验已知，d满足d≤λ/2以避免相位模糊，λ为窄带信号的波长， 则辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，s(t)为K×1维的入射信号矢量，n(t) 为L×1维的加性高斯白噪声矢量，a(θ)为A(θ)的列向量，表示为：

a(z)＝[1,z,z²,…,z^L-1]^T,

其中，z＝jφ，φ＝(2π/λ)dsinθ；

L×L维阵列协方差矩阵为：

其中R_ss＝E[s(t)s^H(t)]是K×K的信号协方差矩阵，理论上R_xx不可知，但 可以用T快拍数据对其估计：

其复值特征值分解可以表示为：

本发明所述计算信号实协方差矩阵具体为：

(1)用前后向平滑协方差矩阵代替

从而得到 对称实协方差阵：

其中U是酉矩阵，由L的奇 偶性其具有两种形式：

为了简化而不失一般性，文中假设L＝2M,实值协方差矩阵的特征值分解为：

其中，

是实信号矩阵，

是实噪声矩阵，因此波达 角可以通过低复杂度实值计算进行估计；

U-MUSIC的谱函数为：

其中a(φ)是2M×1的实导向矢量：

本发明所述构造实导向矢量具体为：

定义新的变量ψ＝φ/2，则用ψ表示的导向矢量为：

其中，b(ψ)和d(ψ)是两个M×1维的实向量：

再定义两个变量x＝cosψ和 y＝sinψ，通过三角多倍角公式可以表示：

其中α_m,k和β_m,k如下：

由此可以计算得到：

cos(3ψ)＝x(4x²-3)

sin(3ψ)＝y(4x²-1)

cos(5ψ)＝x(16x⁴-20x²+5)

sin(5ψ)＝y(16x⁴-12x²+1)；

用α_m,k和β_m,k表示b(ψ)和d(ψ)：

其中，

是M×1维范德蒙德向量：

并且

和

是两个由α_m,k和β_m,k表示的查询矩阵：

本发明所述构造实值多项式具体为：将噪声矩阵

分为四个子矩阵：

由此将U-MUSIC的功率谱表示为ψ的多项式:

再将其用x，

和y表示为：

其 中，

并且

和

都是2M-2阶多项式：

其中，γ_k，δ_k和ε_k都是待定系数，并且

和

都是M×M维实矩阵：

因为可以通过对

求根进行波达角估计，为了建立只有一个变量的多 项式，所以令

有如下关系：

其中

是2M-1阶的：

其中，ζ_k是待定系数，并且对y的表达式两边同时平方，并带入

得到只有

的4M-2阶多项式；

本发明所述计算多项式系数及根具体包括以下步骤：

(1)通过对

和

的观察，得到：

因此γ_k,k∈[0,2M-2]可以由此计算：

s.t.s+t-2＝k， 同样的，可以计算出：

s.t.s+t-2＝k

s.t.s+t-2＝k；

由以上结果可以确定ζ_k：

然后令

由此可以计算：

s.t.s+t＝k

s.t.s+t＝k，

最终可以确定系数：

计算出所有系数后，求解方程

确定方程 的根。

本发明中解模糊并获得信号波达方向包括以下内容：

因为根是最靠近实轴的，通过

得到可能的2K个根，进而估计出角度：

通过最大化

解模糊，得到K个真实的DOAs的估计值。

本发明与现有技术相比有益效果为：本发明实现了特征值分解和多项式求根 的全实值计算，避免了进行复值多项式求根所需的大量计算，为波达方向的工程 化实现提供了技术支持。

附图说明

图1是本发明的流程图。

图2是本发明与求根MUSIC算法和酉求根MUSIC算法根的分布图，其中L＝8， SNR＝10dB,T＝100,K＝2,θ₁＝6°,θ₂＝15°。

图3是本发明算法根的分布图，其中L＝8，SNR＝10dB,T＝100,K＝2,θ₁＝6°, θ₂＝15°。

图4是本发明与不同算法的RMSE随输入信噪比的变化情况，其中L＝8，SNR＝10dB,T＝100,K＝2,θ₁＝6°,θ₂＝15°。

图5是本发明与不同算法的RMSE随快拍数的变化情况，其中L＝8，SNR＝10dB, T＝100,K＝2,θ₁＝6°,θ₂＝15°。

图6是本发明与不同算法运行时间对比图。

具体实施方式

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

如附图1所示，本发明提出了一种基于全实值计算的多项式求根波达角估计 新方法，具体包括以下内容：

第一步，利用天线阵列接收辐射源信号：

(1)假设L个相互独立的阵元，以d等间距组成均匀线阵(ULA)，考虑空 间中存在K个远场窄带信号入射到阵列。其中，假设K先验已知，d满足d≤λ/2 以避免相位模糊，λ为窄带信号的波长。阵列接收辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，s(t)为K×1维的入射信号矢量，n(t) 为L×1维的加性高斯白噪声矢量，a(θ)为A(θ)的列向量，可以表示为：

a(z)＝[1,z,z²,…,z^L-1]^T,

其中，z＝jφ，φ＝(2π/λ)d sin θ；

(2)L×L维阵列协方差矩阵为：

其中R_ss＝E[s(t)s^H(t)]是K×K的信号协方差矩阵，理论上R_xx不可知，但 可以用T快拍数据对其估计：

其复值特征值分解可以表示为：

第二步，计算信号实协方差矩阵，所述第二步包括以下步骤：

(1)用前后向平滑协方差矩阵代替

从而得到对称实协方差矩阵：

其中U是酉矩阵，由L的奇偶性其具有两种形式：

为了简化而不失一般性，文中假设L＝2M,实值协方差矩阵的特征值分解为：

其中，

是实信号矩阵，

是实噪声矩阵，因此波达角可以通过低复杂度实 值计算进行估计；

(2)U-MUSIC的谱函数为：

其中a(φ)是2M×1的实导向矢量：

第三步，构造实导向矢量，所述第三步包括以下步骤：

(1)定义新的变量ψ＝φ/2，则用ψ表示的导向矢量为：

其中，b(ψ)和d(ψ)是两个M×1维的实向量：

再定义两个变量x＝cosψ和y＝sinψ，通过三角多倍角公式可以表示：

其中α_m,k和β_m,k如下：

由此可以计算得到：

cos(3ψ)＝x(4x²-3)

sin(3ψ)＝y(4x²-1)

cos(5ψ)＝x(16x⁴-20x²+5)

sin(5ψ)＝y(16x⁴-12x²+1)；

(2)用α_m,k和β_m,k表示b(ψ)和d(ψ)：

其中，

是M×1维范德蒙德向量：

并且

和

是两个由α_m,k和β_m,k表示的查询矩阵：

第四步，构造实值多项式，所述第四步包括以下步骤：

(1)将噪声矩阵

分为四个子矩阵：

由此可以将U-MUSIC的功率谱表示为ψ的多项式:

(2)再将其用x，

和y表示为：

其中，

并且

和

都是2M-2阶多项式：

其中，γ_k，δ_k和ε_k都是待定系数，并且

和

都是M×M维实矩阵：

因为可以通过对

求根进行波达角估计，为了建立只有一个变量的多 项式，所以令

有如下等式：

其中

是2M-1阶的：

其中，ζ_k是待定系数，并且对y的表达式两边同时平方，并带入

得到只有

的4M-2阶多项式。

第五步，计算多项式系数及根，所述第五步包括以下步骤：

(1)通过对

和

的观察，可以得到：

因此γ_k,k∈[0,2M-2]可以由此计算：

s.t.s+t-2＝k，

同样的，可以计算：

s.t.s+t-2＝k

s.t.s+t-2＝k；

(2)由以上结果可以确定ζ_k：

然后令

由此可以计算：

s.t.s+t＝k

s.t.s+t＝k，

最终可以确定系数：

(3)计算出所有系数后，求解方程

确 定方程的根。

第六步，解模糊并获得信号波达方向，所述第六步包括以下步骤：

(1)因为

因为根是最靠近实轴的，可以通过

得到可能的2K个根，进而估计出角度：

(2)通过最大化

解模糊，得到K个真实的DOAs的估计值。

本发明的性能可通过以下仿真说明：

1.仿真条件

假设采用12阵元的阵元间距为d＝λ/2的ULA阵型，两个入射信号的方向为 θ₁＝6°和θ₂＝15°。为了进一步评价本发明的性能，设置蒙特卡洛实验次数为1000， 将均方根误差(Root Mean Square Error,RMSE)作为评价指标，并引入克拉 美罗界(Cramér-Rao LowerBound,CRLB)作为对比。

2.仿真内容和结果

仿真1，设置阵元数L＝8，信源数K＝2，比较本发明与求根MUSIC算法和 酉求根MUSIC算法根的分布，其结果如图3所示。

由图2和图3可知，求根MUSIC算法和酉求根MUSIC算法需要计算14个根， 其中两对根贴近单位圆，同样本发明也需要计算14个根，其中两对模糊根贴近 实轴。这个结果证实了此前分析的正确性。

仿真2，设置T＝100，比较本发明与不同算法的RMSE随输入信噪比 (Signal-to-Noise Ratio,SNR)的变化情况，其结果如图4所示。

由图4可知，酉求根MUSIC算法在SNR≤-5dB的条件下的性能优于本发明， 但是随着SNR增加，本发明和其他两种算法性能接近。

仿真3，设置信噪比SNR＝0dB，比较本发明与不同算法的RMSE随快拍数 的变化情况，其结果如图5所示。

酉求根MUSIC算法在快拍数T≤50性能优于本发明，随着T的增加，本发明 的RMSE显著减少，最终本发明和其他两种算法性能接近，并且具有接近于CRLB 的估计精度。

仿真4，比较本发明与不同算法的计算效率，其中通过在Intel(R)Core(TM) DuoT5870 2.0GHz CPU处理器和1GB RAM内存的同一PC环境中运行MATLAB代 码，从CPU时间的角度对计算效率进行等效评估。

由图6可知，求根MUSIC算法的运行时间为1.7940e-4秒，酉求根MUSIC算 法的运行时间为1.7380e-4秒，最后本发明的运行时间为0.6350e-4秒，因此本 发明相比另外两种算法计算效率更高，并且可以提供接近CRLB的估计精度。

本发明与求根MUSIC算法和酉求根MUSIC算法计算复杂度，如表1所示：

表1

其中O(L³)表示对一个L×L的实矩阵进行特征值分解和对L阶实值多项式进 行求根的计算复杂度，因为求根MUSIC算法是复值算法，所以在EVD和多项式求 根过程需要4×O(L³)和32×O(L³)的计算复杂度，酉求根MUSIC算法因为在EVD 过程只涉及实值计算因此相比与求根MUSIC算法在该过程减少到了四分之一的 计算量，对于本发明，只涉及实值计算，所以在EVD和多项式求根过程分别只需 要O(L³)和8×O(L³)的计算量，但是对于系数的计算需要额外的步骤，这部分的 计算复杂度约为5/4×O(L³)，而对矩阵的计算也同样需要1/4×O(L³)的计算复 杂度。由表1可知，本发明相对其他两种算法显著降低了计算复杂度。

Claims (7)

1.一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，获取辐射源信号后，计算信号协方差矩阵，构造实导向矢量后，构造实值多项式，计算多项式系数及根，解模糊并获得信号波达方向。

2.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，获取辐射源信号具体包括以下内容：

利用天线阵列接收辐射源信号，天线阵列具有L个相互独立的阵元，以d等间距组成均匀线阵(ULA)，空间中存在K个远场窄带信号入射到阵列，其中，假设K先验已知，d满足d≤λ/2以避免相位模糊，λ为窄带信号的波长，则辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，s(t)为K×1维的入射信号矢量，n(t)为L×1维的加性高斯白噪声矢量，a(θ)为A(θ)的列向量，表示为：

a(z)＝[1,z,z²,…,z^L-1]^T,其中，z＝jφ，φ＝(2π/λ)d sinθ；

L×L维阵列协方差矩阵为：

其中R_ss＝E[s(t)s^H(t)]是K×K的信号协方差矩阵，理论上R_xx不可知，但可以用T快拍数据对其估计：

其复值特征值分解可以表示为：

3.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，计算信号实协方差矩阵具体为：

用前后向平滑协方差矩阵代替

得到对称实协方差阵：

其中U是酉矩阵，由L的奇偶性其具有两种形式：

为了简化而不失一般性，假设L＝2M,实值协方差矩阵的特征值分解为：

其中，

是实信号矩阵，

是实噪声矩阵，因此波达角通过低复杂度实值计算进行估计；

U-MUSIC的谱函数为：

其中a(φ)是2M×1的实导向矢量：

4.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，所述构造实导向矢量具体为：

定义变量ψ＝φ/2，则用ψ表示的导向矢量为：

其中，b(ψ)和d(ψ)是两个M×1维的实向量：

再定义两个变量x＝cosψ和y＝sinψ，通过三角多倍角公式表示：

其中α_m,k和β_m,k如下：

由此可以计算得到：

cos(3ψ)＝x(4x²-3)

sin(3ψ)＝y(4x²-1)

cos(5ψ)＝x(16x⁴-20x²+5)

sin(5ψ)＝y(16x⁴-12x²+1)；

用α_m,k和β_m,k表示b(ψ)和d(ψ)：

其中，

是M×1维范德蒙德向量：

并且

和

是两个由α_m,k和β_m,k表示的查询矩阵：

5.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，构造实值多项式具体为：将噪声矩阵

分为四个子矩阵：

由此将U-MUSIC的功率谱表示为ψ的多项式:

再将其用x，

和y表示为：

其中，

并且

和

都是2M-2阶多项式：

其中，γ_k，δ_k和ε_k都是待定系数，并且

和

都是M×M维实矩阵：

通过对

求根进行波达角估计，为了建立只有一个变量的多项式，所以令

有如下关系：

其中

是2M-1阶的：

其中，ζ_k是待定系数，并且对y的表达式两边同时平方，并带入

得到只有

的4M-2阶多项式。

6.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，计算多项式系数及根具体包括以下步骤：

通过

和

得到：

因此γ_k,k∈[0,2M-2]由此计算：

同样的，计算出：

由以上结果确定ζ_k：

然后令

由此计算：

最终确定系数：

计算出所有系数后，求解方程

确定方程的根。

7.根据权利要求1所述的一种基于全实值计算的远场窄带无线电信号波达方向估计方法，其特征在于，解模糊并获得信号波达方向包括以下内容：

因为根是最靠近实轴的，通过

得到可能的2K个根，进而估计出角度：

通过最大化

解模糊，得到K个真实的DOAs的估计值。

Images