CN106569180B - 一种基于Prony方法的方位估计算法 - Google Patents

Abstract

本发明公开了一种基于Prony方法的方位估计算法，该方法适用于在雷达和声纳系统中对目标方位估计，其特征在于包括以下步骤：步骤1、获取DOA估计算法测量方程；步骤2、对于均匀线性阵列，直接利用Prony方法求解测量方程；步骤3、对于均匀圆形阵列，构造变换矩阵，再利用Prony方法求解测量方程；本发明的技术方案实现了在均匀线性阵列和均匀圆形阵列测量平台上对目标的DOA估计，在高信噪比条件下估计性能良好，无需进行搜索，不用进行特征分解从而减小了计算量，同时能够突破阵列的瑞利限，分辨率较高，并且在相干信源的条件下不用进行解相干的操作，为方位估计求解算法提供了一种新的方法。

Description

一种基于Prony方法的方位估计算法

技术领域

本发明属于雷达和声纳系统方位估计算法领域，涉及一种基于Prony方法的方位估计算法。

背景技术

波达方向(DOA：Direction of Arrival)估计是雷达和声纳系统中的关键问题之一，在目标探测、导航定位和无线通信等领域有着广泛的应用。最早的DOA估计算法称为常规波束形成(CBF)，对方位估计的角度分辨力受到瑞利限的限制；进而发展了一些高分辨DOA估计算法如线性预测算法、多重信号分类算法、最大似然算法、旋转不变子空间算法和压缩感知算法等，这些算法大多需要对谱峰进行搜索，或者构建特征多项式直接求解。在运算过程中，有的算法还存在一个特征值分解的过程，因此高分辨率算法的计算量较大；同时在相干信源的条件下，这些DOA估计算法效果不是很好，需要进行相应的解相干操作，进一步加大了算法的计算量。

发明内容

本发明所要解决的技术问题是，克服现有技术的缺点，提供一种基于Prony方法的方位估计算法，本发明方法基于Prony理论提出了一种新的DOA估计算法，该方位估计算法具有在高信噪比条件下估计性能良好，无需进行搜索，不用进行特征分解从而减小了计算量，同时能够突破阵列的瑞利限分辨率较高等优点，并且在相干信源的条件下不用进行解相干的操作，为方位估计求解算法提供了一种新的思路。

为了解决原有DOA估计算法的不足之处，本发明提供一种基于Prony方法的方位估计算法，包括如下步骤：

步骤(1)获取DOA估计算法测量方程：

在二维DOA方位测量上进行讨论，令阵元编号为k＝0,1,…,M-1，各信源方位为θ_i(i＝1,2,…,p)，信源为窄带信号，如式(1)所示：

其中，测量系统的采样频率为f_s，采样间隔Ts＝1/f_s，v(nT_s)为测量噪声。以信源i传播到阵元0上经过的传播时间为基准，τ_k(θ_i)是阵元k接收到信源i的时间相对于阵元0接收到信源i的时间差值；令信源的中心频率为f_c，取快拍数N满足式(2)：

则可以得到各阵元接收信号格点频率上的DFT值，如式(3)所示：

将式(1)代入式(3)中，忽略噪声v(nT_s)，可以得到式(4)：

于是可得测量方程，如式(5)所示：

步骤(2)对于均匀线性阵列，直接利用Prony方法求解测量方程：

测量方程(5)即为基于Prony方法的DOA估计算法测量模型，求解该方程即可以得到各信源的方位角θ_i，同时也能得出各信源的幅度参数a_i。

对于均匀线性阵列下的方位估计求解方法如下：

阵列是均匀线性阵列时，已知延时差如式(6)所示：

τ_k(θ_i)＝kdsinθ_i/c (6)

其中，d是阵元间距，c是传播速度，将式(6)代入式(5)，可以得出式(7)：

此时，DOA估计问题成为一个频率估计问题，若能解出z_i，则各信源的方位和幅度如式(8)所示：

θ_i＝arcsin(-arg(z_i)*λ/(2πd)) a_i＝|A_i| (8)

式(7)即为基本的Prony模型，可用Prony算法进行求解，X_k满足一阶AR模型(9)：

X_k+c₁X_k-1+…+c_pX_k-p＝0 (9)

可以根据式(10)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(11)所示：

求解方程(10)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(12)所示：

1+c₁z^-1+…+c_pz^-p＝0 (12)

式(12)的根z_i，也称其为Prony极点，于是式(5)模型可简化为参数A_i的线性方程，用矩阵形式表示即为式(13)所示：

ZA＝X (13)

其中极点矩阵Z和测量矩阵A如式(14)所示：

最后可以得出测量矩阵如式(15)所示：

A＝Z⁺X＝(Z^HZ)^-1Z^HX (15)

总结以上的推导，均匀线性阵列的Prony方位估计方法步骤如下：

(1)已知信源个数，根据式(3)得到各阵元接收信号的格点频率的DFT值X_k；

(2)利用式(11)计算样本函数r(m,n)并构造矩阵R，求解方程(10)得到系数c_i的估计值；

(3)求出特征多项式(12)的根z_i，并利用公式(8)解算出信源方位θ_i；

(4)利用式(13)～(15)解算出参数A_i，进而得出信源i的幅度a_i；

存在问题：根据式(8)，要求：

为了满足所有的角度，要求d/λ≥1/4，当d/λ＝1/2时，仅能够测量-60°～60°之间的声源。

步骤(3)对于均匀圆形阵列，构造变换矩阵，利用Prony方法求解测量方程。均匀圆阵下的方位估计求解方法如下：

阵列是均匀圆形阵列时，已知延时差如式(17)所示：

τ_k(θ_i)＝rcos(kθ_A-θ_i)/c (17)

其中r是圆阵半径，平均角度θ_A＝2π/M，M为阵元个数，可以得到格点频率上的DFT值如式(18)所示：

令β＝ωr/c，已知：

可以推出式(20)：

当|m|>β时，J_m(β)≈0，式(20)可化为式(21)

于是可以得到式(22)：

其中f_k(m)和变换矩阵Y如式(23)所示：

因此可以得到式(24)：

Y＝F⁺X＝(F^HF)^-1F^HX (24)

并且变换矩阵Y的各项如式(25)所示：

式(25)与式步骤(2)中的式(7)结构一致，同样的，DOA估计问题成为一个频率估计问题，若能估计出z_i和A_i，则各信源的方位和幅度如式(26)所示：

θ_i＝-arg(z_i) a_i＝|A_i| (26)

可用Prony算法进行求解，Y_k满足一阶AR模型(27)：

Y_k+c₁Y_k-1+…+c_pY_k-p＝0 (27)

可以根据式(28)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(29)所示：

求解方程(28)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(30)所示：

1+c₁z^-1+…+c_pz^-p＝0 (30)

式(30)的根z_i也是Prony极点，其中极点矩阵Z和测量矩阵A如式(31)所示：

最后可以得出测量矩阵如式(32)所示：

A＝Z⁺Y＝(ZHZ)^-1Z^HY (32)

因此，可以得出均匀圆形阵列的Prony方位估计算法步骤如下：

(1)已知信源个数，根据式(18)得到各阵元接收信号的格点频率的DFT值X_k；

(2)根据式(23)得到f_k(m)构造矩阵F，同时根据式(24)得到变换矩阵Y；

(3)利用式(29)计算样本函数r(m,n)并构造矩阵R，求解方程(28)得到系数c_i的估计值；

(4)求出特征多项式(30)的根z_i，并利用公式(26)解算出信源方位θ_i；

(5)利用式(31)～(32)解算出参数A_i，进而得出信源i的幅度a_i；

根据式(26)，该算法对于r/λ并无要求。

本发明的有益效果是：具有在高信噪比条件下估计性能良好，无需进行搜索，不用进行特征分解从而减小了计算量，同时能够突破阵列的瑞利限分辨率较高等优点，并且在相干信源的条件下不用进行解相干的操作。

附图说明

图1为本发明的原理流程图；

图2为本发明估计算法的测量模型建立流程图；

图3为仿真实验信源生成流程图；

图4为均匀线性阵列下的方位估计求解流程图；

图5为均匀圆形阵列下的方位估计求解流程图。

具体实施方式

实施例1

本实施例提供一种基于Prony方法的方位估计算法，原理如图1所示，包括如下步骤：

步骤(1)获取DOA估计算法测量方程：

在二维DOA方位测量上进行讨论，令阵元编号为k＝0,1,…,M-1，各信源方位为θ_i(i＝1,2,…,p)，信源为窄带信号，如式(1)所示：

其中，测量系统的采样频率为f_s，采样间隔Ts＝1/f_s，v(nT_s)为测量噪声。以信源i传播到阵元0上经过的传播时间为基准，τ_k(θ_i)是阵元k接收到信源i的时间相对于阵元0接收到信源i的时间差值；令信源的中心频率为f_c，取快拍数N满足式(2)：

则可以得到各阵元接收信号格点频率上的DFT值，如式(3)所示：

将式(1)代入式(3)中，忽略噪声v(nT_s)，可以得到式(4)：

于是可得测量方程，如式(5)所示：

步骤(2)对于均匀线性阵列，直接利用Prony方法求解测量方程：

测量方程(5)即为基于Prony方法的DOA估计算法测量模型，求解该方程即可以得到各信源的方位角θ_i，同时也能得出各信源的幅度参数a_i。

对于均匀线性阵列下的方位估计求解方法如下：

阵列是均匀线性阵列时，已知延时差如式(6)所示：

τ_k(θ_i)＝kdsinθ_i/c (6)

其中，d是阵元间距，c是传播速度，将式(6)代入式(5)，可以得出式(7)：

此时，DOA估计问题成为一个频率估计问题，若能解出z_i，则各信源的方位和幅度如式(8)所示：

θ_i＝arcsin(-arg(z_i)*λ/(2πd)) a_i＝|A_i| (8)

式(7)即为基本的Prony模型，可用Prony算法进行求解，X_k满足一阶AR模型(9)：

X_k+c₁X_k-1+…+c_pX_k-p＝0 (9)

可以根据式(10)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(11)所示：

求解方程(10)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(12)所示：

1+c₁z^-1+…+c_pz^-p＝0 (12)

式(12)的根z_i，也称其为Prony极点，于是式(5)模型可简化为参数A_i的线性方程，用矩阵形式表示即为式(13)所示：

ZA＝X (13)

其中极点矩阵Z和测量矩阵A如式(14)所示：

最后可以得出测量矩阵如式(15)所示：

A＝Z⁺X＝(Z^HZ)^-1Z^HX (15)

根据式(8)，要求：

为了满足所有的角度，要求d/λ≥1/4，当d/λ＝1/2时，仅能够测量-60°～60°之间的声源。

步骤(3)对于均匀圆形阵列，构造变换矩阵，利用Prony方法求解测量方程。均匀圆阵下的方位估计求解方法如下：

阵列是均匀圆形阵列时，已知延时差如式(17)所示：

τ_k(θ_i)＝rcos(kθ_A-θ_i)/c (17)

其中r是圆阵半径，平均角度θ_A＝2π/M，M为阵元个数，可以得到格点频率上的DFT值如式(18)所示：

令β＝ωr/c，已知：

可以推出式(20)：

当|m|>β时，J_m(β)≈0，式(20)可化为式(21)

于是可以得到式(22)：

其中f_k(m)和变换矩阵Y如式(23)所示：

因此可以得到式(24)：

Y＝F⁺X＝(F^HF)^-1F^HX (24)

并且变换矩阵Y的各项如式(25)所示：

式(25)与式步骤(2)中的式(7)结构一致，同样的，DOA估计问题成为一个频率估计问题，若能估计出z_i和A_i，则各信源的方位和幅度如式(26)所示：

θ_i＝-arg(z_i) a_i＝|A_i| (26)

可用Prony算法进行求解，Y_k满足一阶AR模型(27)：

Y_k+c₁Y_k-1+…+c_pY_k-p＝0 (27)

可以根据式(28)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(29)所示：

求解方程(28)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(30)所示：

1+c₁z^-1+…+c_pz^-p＝0 (30)

式(30)的根z_i也是Prony极点，其中极点矩阵Z和测量矩阵A如式(31)所示：

最后可以得出测量矩阵如式(32)所示：

A＝Z⁺Y＝(Z^HZ)^-1Z^HY (32)

以下将结合附图1-4和仿真实验对本发明的技术方案进行详细描述。

图2为本发明估计算法的测量模型建立流程图，图3为仿真实验的信源设计流程，可分别按照均匀线性阵列和均匀圆形阵列阵元间的延时差τk(θi)进行生成，下面对仿真实验结果进行介绍。

考虑均匀线性阵列，数据快拍数N＝1000，阵元数M＝10，阵元间距d＝λ/2，已知瑞利限为10°，假设两个相关的声源，幅度为1，信噪声比30dB。使用prony算法测量1000次，按图4的求解流程，测量结果如表1所示：

表1两个相关声源在均匀线阵中的测量结果

可以看出，对于两个相关声源的情况，其方位分辨率约为3度，突破了瑞利门限。

然后考虑均匀圆形阵列，取快拍数N＝1000，阵元个数M＝32，r＝2λ，已知瑞利限约为30度。存在两个相关的声源，幅度为1，信噪比为30dB。使用prony算法测量1000次，按图5的求解流程，测量结果如表2所示：

表2两个相关声源在均匀圆阵中的测量结果

同样的，对于均匀圆形阵形而言，该算法也同样对相关信号突破了瑞利限。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (2)

1.一种基于Prony方法的方位估计算法，其特征在于，包括如下步骤：

步骤(1)获取DOA估计算法测量方程：

基于二维DOA方位测量，令阵元编号为k＝0,1,…,M-1，各信源方位为θ_i(i＝1,2,…,p)，信源为窄带信号，如式(1)所示：

其中，测量系统的采样频率为f_s，采样间隔Ts＝1/f_s，v(nT_s)为测量噪声；

以信源i传播到阵元0上经过的传播时间为基准，τ_k(θ_i)是阵元k接收到信源i的时间相对于阵元0接收到信源i的时间差值；

令信源的中心频率为f_c，取快拍数N满足式(2)：

则可得到各阵元接收信号格点频率上的DFT值，将X_k(Nf_c/f_s)简写定义X_k，如式(3)所示：

将式(1)代入式(3)中，忽略噪声v(nT_s)，可以得到式(4)：

于是可得测量方程，如式(5)所示：

步骤(2)对于均匀线性阵列，直接利用Prony方法求解测量方程：

测量方程(5)即为基于Prony方法的DOA估计算法测量模型，求解该方程即可以得到各信源的方位角θ_i，同时也能得出各信源的幅度参数a_i；

对于均匀线性阵列下的方位估计求解方法如下：

阵列是均匀线性阵列时，已知延时差如式(6)所示：

τ_k(θ_i)＝kd sinθ_i/c (6)

其中，d是阵元间距，c是传播速度，将式(6)代入式(5)，可以得出式(7)：

此时，DOA估计问题成为一个频率估计问题，若能解出z_i，则各信源的方位和幅度如式(8)所示：

θ_i＝arcsin(-arg(z_i)*λ/(2πd)) a_i＝|A_i| (8)

式(7)即为基本的Prony模型，可用Prony算法进行求解，X_k满足一阶AR模型(9)：

X_k+c₁X_k-1+…+c_pX_k-p＝0 (9)

可以根据式(10)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(11)所示：

求解方程(10)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(12)所示：

1+c₁z^-1+…+c_pz^-p＝0 (12)

式(12)的根z_i，也称其为Prony极点，于是式(5)模型可简化为参数A_i的线性方程，用矩阵形式表示即为式(13)所示：

ZA＝X (13)

其中极点矩阵Z和测量矩阵A如式(14)所示：

最后可以得出测量矩阵如式(15)所示：

A＝Z⁺X＝(Z^HZ)^-1Z^HX (15)

步骤3、对于均匀圆形阵列，构造变换矩阵，利用Prony方法求解测量方程，均匀圆阵下的方位估计求解方法如下：

阵列是均匀圆形阵列时，已知延时差如式(16)所示：

τ_k(θ_i)＝r cos(kθ_A-θ_i)/c (16)

其中r是圆阵半径，平均角度θ_A＝2π/M，M为阵元个数，可以得到格点频率上的DFT值如式(17)所示：

令β＝ωr/c，已知：

可以推出式(19)：

当|m|>β时，J_m(β)≈0，式(19)可化为式(20)

于是可以得到式(21)：

其中f_k(m)和变换矩阵Y如式(22)所示：

因此可以得到式(23)：

Y＝F⁺X＝(F^HF)^-1F^HX (23)

并且变换矩阵Y的各项如式(24)所示：

式(24)与式步骤(2)中的式(7)结构一致，同样的，DOA估计问题成为一个频率估计问题，若能估计出z_i和A_i，则各信源的方位和幅度如式(25)所示：

θ_i＝-arg(z_i) a_i＝|A_i| (25)

式(7)即为基本的Prony模型，可用Prony算法进行求解，Y_k满足一阶AR模型(26)：

Y_k+c₁Y_k-1+…+c_pY_k-p＝0 (26)

可以根据式(27)求解系数c_i：

其中e_p是AR过程中的最小误差，并且相关函数如式(28)所示：

求解方程(27)，可得到系数c_i和最小误差e_p的估计值，即可建立特征多项式如式(29)所示：

1+c₁z^-1+…+c_pz^-p＝0 (29)

式(29)的根z_i也是Prony极点，其中极点矩阵Z和测量矩阵A如式(30)所示：

最后可以得出测量矩阵如式(31)所示：

A＝Z⁺Y＝(Z^HZ)^-1Z^HY (31)。

2.根据权利要求1所述的基于Prony方法的方位估计算法，其特征在于，所述步骤(2)中根据式(32)，要求：

为了满足所有的角度，要求d/λ≥1/4，当d/λ＝1/2时，仅能够测量-60°～60°之间的声源。

Images