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

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

背景技术：

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

技术实现要素：

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

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

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

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

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

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

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

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

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

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

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

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

τk(θi)＝kdsinθi/c(6)

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

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

θi＝arcsin(-arg(zi)*λ/(2πd))ai＝|ai|(8)

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

xk+c1xk-1+…+cpxk-p＝0(9)

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

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

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

1+c1z^-1+…+cpz^-p＝0(12)

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

za＝x(13)

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

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

a＝z⁺x＝(z^hz)^-1z^hx(15)

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

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

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

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

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

存在问题：根据式(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|>β时，jm(β)≈0，式(20)可化为式(21)

于是可以得到式(22)：

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

因此可以得到式(24)：

y＝f⁺x＝(f^hf)^-1f^hx(24)

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

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

θi＝-arg(zi)ai＝|ai|(26)

可用prony算法进行求解，yk满足一阶ar模型(27)：

yk+c1yk-1+…+cpyk-p＝0(27)

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

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

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

1+c1z^-1+…+cpz^-p＝0(30)

式(30)的根zi也是prony极点，其中极点矩阵z和测量矩阵a如式(31)所示：

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

a＝z⁺y＝(zhz)^-1z^hy(32)

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

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

(2)根据式(23)得到fk(m)构造矩阵f，同时根据式(24)得到变换矩阵y；

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

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

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

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

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

附图说明

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

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

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

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

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

具体实施方式

实施例1

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

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

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

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

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

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

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

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

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

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

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

τk(θi)＝kdsinθi/c(6)

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

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

θi＝arcsin(-arg(zi)*λ/(2πd))ai＝|ai|(8)

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

xk+c1xk-1+…+cpxk-p＝0(9)

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

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

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

1+c1z^-1+…+cpz^-p＝0(12)

式(12)的根zi，也称其为prony极点，于是式(5)模型可简化为参数ai的线性方程，用矩阵形式表示即为式(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|>β时，jm(β)≈0，式(20)可化为式(21)

于是可以得到式(22)：

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

因此可以得到式(24)：

y＝f⁺x＝(f^hf)^-1f^hx(24)

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

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

θi＝-arg(zi)ai＝|ai|(26)

可用prony算法进行求解，yk满足一阶ar模型(27)：

yk+c1yk-1+…+cpyk-p＝0(27)

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

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

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

1+c1z^-1+…+cpz^-p＝0(30)

式(30)的根zi也是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两个相关声源在均匀圆阵中的测量结果

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

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