CN109061551B - 一种基于多项式求根的无网格稀疏谱估计方法 - Google Patents

Abstract

本发明涉及一种基于多项式求根的无网格稀疏谱估计方法，通过求解阵列接收信号的协方差矩阵，将所有快拍数进行结合，使至方法在目标运动缓慢或静止的情况下的定位精度高于单快拍方法，并且能够适用于更低的信噪比情况下。将协方差矩阵重新表示在连续空间上并基于该模型建立连续空间的DOA估计问题，利用半定规划求解该问题从而使DOA估计转化为多项式求根，实现了连续空间上的DOA估计，避免了由于网格划分不够精细而带来的DOA估计误差。

Description

一种基于多项式求根的无网格稀疏谱估计方法

技术领域

本发明属于信号处理等领域，涉及一种基于多项式求根的无网格稀疏谱估计方法，通过将DOA估计问题建立在连续空间上，避免网格划分与信号方位不匹配而引起的DOA估计误差，在目标信号运动缓慢或静止的情况下实现高精度定位。

背景技术

阵列信号处理在雷达、声纳等领域上具有广泛的应用，目标方位(Direction ofarrival,DOA)估计是阵列信号处理一大主要任务。常规的DOA估计方法有最小方差无失真响应方法(Minimum variance distortion response beamforming,MVDR)和多重信号分类方法(Multiple signal classification,MUSIC)，这些方法均能够实现高分辨性，但对快拍数要求较高。稀疏信号处理类的DOA估计方法是近十年发展起来的DOA估计方法，这种方法可用于小快拍和低信噪比情况下，其性能远远优于常规的DOA估计方法。

常见的稀疏信号处理方法在进行DOA估计前一般先将观测空间划分为离散的网格，在网格上建立接收信号模型并进行DOA估计。这类方法有一大缺陷，即当信号的真实方位与所划分的网格点不一致时，DOA估计的结果会存在一定误差。虽然增加空间网格的划分精度在一定程度上能够减轻目标方位与网格点不匹配的问题，但却大大增加了运算量。

国外学者Xenaki等人(Angeliki Xenaki,Peter Gerstoft.Grid-freecompressive beamforming[J].The Journal of the Acoustical Society of America,2015,137(4):1923-1935)将稀疏信号处理方法表示在连续空间上，提出了无网格的压缩感知(Grid-free compressive sensing,GFCS)方法，通过求解半正定规划问题将DOA估计转化为多项式求根，以实现连续空间上的DOA估计。尽管该方法一定程度上减弱了目标与网格点不匹配的问题，但仅适用于单快拍情况下。而当目标信号所处背景环境为慢变的，且目标运动较为缓慢或者静止时，通过将多个快拍进行结合通常能够给出精度更高的DOA估计结果，并且可应用于信噪比更低的环境下。除此之外，采用GFCS在噪声环境下进行DOA估计需选择一个超参数，该参数与阵列接收噪声幅值有关。由于环境噪声一般为随机信号，其幅值变化没有规律，使得该参数的选取较为困难，不合适的参数选取将大大影响GFCS的性能。因此，需要建立一个合适的多快拍模型，并将其表示在连续空间上，实现多快拍情况下的无网格DOA估计。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于多项式求根的无网格稀疏谱估计方法，在慢变环境下实现对静止或运动缓慢的目标信号更为精准的DOA估计，同时解决GFCS中超参数难以选取的问题，解决多快拍情况下的无网格DOA估计问题。

技术方案

一种基于多项式求根的无网格稀疏谱估计方法，其特征在于估计步骤如下：

步骤1：采用阵元间距为半波长的M元均匀线列阵接收窄带信号,均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),1≤n≤N,i＝1,...,M；

将观测空间[-90°,90°]划分为Q个网格；所述90°为端射方向，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]；

在该网格上，阵列的接收信号模型表示为x(n)＝A(Θ)s(n)+e(n),n＝1,...,N

其中：

和

分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，a(θ)＝[1e^{-j2πdsin(θ)f/c}...e^{-j2πd(M-1)sin(θ)f/c}]^T，f为窄带信号的中心频率，c为声速；

阵列接收信号的采样协方差矩阵为R＝E{x(n)x^H(n)}＝A(Θ)R_sA^H(Θ)+R_e

其中：E{·}为期望算子；R_s和R_e分别为信号协方差矩阵和噪声协方差矩阵；上标“H”为共轭转置符号；协方差矩阵由采样协方差矩阵

所代替；

对协方差矩阵两端进行向量化得

其中：

vec(·)为矩阵向量化算子；矩阵

的第i列为

上标“*”表示求共轭，

表示Kronecker积；p＝diag(R_s)为R_s对角线元素组成的向量，diag(·)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵；p_e＝vec(R_e)；

保留矩阵

中元素不同的(2M-1)行，并将每列元素按照

的顺序进行排列，组成新的矩阵

采用同样的方式保留并重新排列r和p_e中的元素，得到新的向量r′和p′_e，则新的协方差矩阵模型表示为

步骤2：将p表示在连续空间上为

其中：t＝sinθ∈[-1,1]；t_i＝sinθ_i,i＝1,2,...,K为真实信号方位所对应的正弦值；p_i为相应的信号功率；δ(t-t_i)为在t_i上的狄拉克函数；

将协方差矩阵模型表示在连续空间上为

其中：p为连续空间上的信号功率；

表示为傅里叶变换算子；r′中第m个元素为

将稀疏谱估计的优化问题表示为：

||p′_e||₂≤ε

||p′_e||₂≤ε

其中：

为原子范数，p的原子范数表示为

||·||₂为l₂范数；

trace(·)表示矩阵求迹；

构造优化问题的对偶问题为：

其中：

为关于

的对偶因子；Re[·]为求实部；

表示半正定符号；

步骤3、构造求根多项式：P(z)＝e^{j2πfd/ct2(M-1)}(1-|H(z)|²)

其中：z＝e^j2πfd/ct；

|·|为求模符号；

求解多项式的根，并选择多项式根中模值为1的根

计算目标信号的方位角：

其中：angle(·)表示求相位符号；

步骤4：计算出所有DOA估计值

重新构造维数为M×K的矩阵

所述

的第k列

得到估计的稀疏谱为：

其中：上标“+”为求Moore-Penrose逆。

有益效果

本发明提出的一种基于多项式求根的无网格稀疏谱估计方法，通过求解阵列接收信号的协方差矩阵，将所有快拍数进行结合，使至方法在目标运动缓慢或静止的情况下的定位精度高于单快拍方法，并且能够适用于更低的信噪比情况下。将协方差矩阵重新表示在连续空间上并基于该模型建立连续空间的DOA估计问题，利用半定规划求解该问题从而使DOA估计转化为多项式求根，实现了连续空间上的DOA估计，避免了由于网格划分不够精细而带来的DOA估计误差。

附图说明

图1：基于多项式求根的无网格稀疏谱估计方法流程

图2：GFSSE进行DOA估计的结果，其中“*”表示真实信号方位

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案包括以下步骤：

1)建立接收信号协方差矩阵模型

采用阵元间距为半波长的M元均匀线列阵作为接收阵列，接收窄带信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),1≤n≤N,i＝1,...,M。将空间[-90°,90°](其中90°为端射方向)划分为Q个网格，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]。在该离散网格上，阵列的接收信号可表示为x(n)＝A(Θ)s(n)+e(n),n＝1,...,N，其中

和

分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，对于均匀线列阵来说，a(θ)＝[1e^{-j2πdsin(θ)f/c}...e^{-j2πd(M-1)sin(θ)f/c}]^T，f为信号频率，d为阵元间距，c为声速。根据阵列接收信号模型，接收信号的协方差矩阵可表示为R＝E{x(n)x^H(n)}＝A(Θ)R_sA^H(Θ)+R_e，将等式两端向量化，可得：

其中r＝vec(R)，vec(·)为矩阵向量化算子；矩阵

的第i列为

上标“*”表示求共轭，

表示Kronecker积；p＝diag(R_s)为R_s对角线元素组成的向量，diag(·)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵；p_e＝vec(R_e)。一般地，协方差矩阵由采样协方差矩阵

所代替。

中仅存在(2M-1)个不同的元素，它们可表示为

当矩阵

中任意两行的元素完全相同时，这两行所提供的信号信息相同，则可舍弃其中一行。保留矩阵

中元素不同的任意(2M-1)行，并将每列元素按照

的顺序进行排列，组成新的矩阵

采用同样的方式保留并重新排列r和p_e中的元素，得到新的向量r′和p′_e，则新的协方差矩阵模型可写为：

2)构造优化问题

为将DOA估计问题转化至连续空间，将p重新表示在连续空间上：

其中t＝sinθ∈[-1,1]，t_i＝sinθ_i,i＝1,2,...,K为真实信号方位所对应的正弦值，p_i为相应的信号功率；δ(t-t_i)为在t_i上的狄拉克函数。r′中第m个元素可表示为：

r′_m和p′_em分别为r′和p′_e中第m个元素。根据式(4)可将式(2)表示在连续空间上，得：

其中

表示为傅里叶变换算子。

定义p的原子范数

为

则GFSSE可表示为：

其中||·||₂为l₂范数；参数ε与噪声功率σ²有关，满足ε≥σ²，用

对ε的值进行估计，trace(·)表示矩阵求迹。

将式(6)的优化问题转化为相应的对偶问题：

其中

和

分别是关于

的对偶因子和关于||p′_e||-ε≤0的拉格朗日乘子；||·||_∞代表无穷范数。

对g(c,μ)关于μ和p′_e求导并置零，可得如下关系式：

将式(8)带入至(7)中，优化问题转化为：

是一个多项式，可写为

该多项式的模值等于1时所对应的t值即为信号方位角的正弦值。从(9)的约束条件可以看出，该多项式具有上限，因此该约束条件可表示为有限维线性矩阵不等式，优化问题可重新写为：

其中Re[·]为求实部；

表示半正定符号。式(10)中待优化变量为有限维的向量，同时约束条件也为有限维的矩阵不等式，因此式(10)可由优化工具箱进行高效求解。

3)求根多项式的构造及DOA估计

当求解完最优变量c后，为方便求解使

模值等于1所对应的t值，构造求根多项式

P(z)＝e^{j2πfd/ct2(M-1)}(1-|H(z)|²) (11)

其中z＝e^j2πfd/ct；|·|为求模符号。计算式(11)所示的多项式的零点，并选择模值为1的零点

则目标信号的方位角可由下式计算：

其中angle(·)表示求相位符号。

4)稀疏谱估计

按照等式(12)计算出所有DOA估计值

重新构造维数为M×K的矩阵

其中

的第k列

则信号的稀疏谱可由下式计算得到：

上标“+”表示求Moore-Penrose逆。

具体实施流程：

基于多项式求根的无网格稀疏谱估计方法流程图见图1，主要过程如下：

1)采用阵元间距为半波长的M元均匀线列阵作为接收阵列，接收窄带信号。均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),1≤n≤N,i＝1,...,M。将空间[-90°,90°](其中90°为端射方向)划分为Q个网格，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]。在该离散网格上，阵列的接收信号可表示为x(n)＝A(Θ)s(n)+e(n),n＝1,...,N。根据阵列接收信号模型，接收信号的协方差矩阵可表示为R＝E{x(n)x^H(n)}＝A(Θ)R_sA^H(Θ)+R_e，将等式两端向量化，可得

一般地，协方差矩阵由采样协方差矩阵

来代替。保留矩阵

中元素不同的任意(2M-1)行，并将每列元素按照

的顺序进行排列，得到新的矩阵

采用同样的方式保留并重新排列r和p_e中的元素，得到新的向量r′和p′_e，新的协方差矩阵的关系式可写为

2)为将DOA估计问题转化至连续空间，将p表示在连续网格上

其中t＝sinθ∈[-1,1]，t_i＝sinθ_i,i＝1,2,...,K为真实的信号方位所对应的正弦值，p_i为相应的信号功率，δ(t-t_i)为在t_i上的狄拉克函数。向量r′中第m个元素可表示为

r′_m和p′_em分别为r′和p′_e的第m个元素，在连续空间上阵列接收信号的协方差矩阵模型可写为

定义p的原子范数

为

则GFSSE可表示为：

||p′_e||₂≤ε

trace(·)表示矩阵求迹。

将该问题转化为相应的对偶问题：

其中

和

分别是关于

的对偶因子和关于||p′_e||-ε≤0的拉格朗日乘子；||·||_∞代表无穷范数。

对g(c,μ)关于μ和p′_e求导并置零，可分别得p′_eo＝c/(2μ_o)和μ_o＝||c||₂/(2ε)，并把它们带入到对偶问题中得到新的优化问题：

是一个多项式，可写为

从约束条件可以看出，该多项式具有上限，因此该约束条件可表示为有限维线性矩阵不等式，则优化问题被重新写为：

该优化问题中待优化变量为有限维的向量，同时约束条件也为有限维的矩阵不等式，因此可由优化工具箱进行高效求解。

3)当求解完最优变量c后，构造多项式P(z)＝e^{j2πfd/ct2(M-1)}(1-|H(z)|²)，其中z＝e^{j2 πfd/ct}。计算该多项式的零点，并寻找模值为1的零点

则目标信号的方位角可表示为

其中angle(·)表示求相位符号。

4)计算出所有的DOA估计值

重新构造维数为M×K的矩阵

其中

的第k列

则信号功率的估计值可计算为

两个不相关的窄带信号分别从-5.6°和-8.6°入射至阵元间距为4m的32元均匀线列阵上,信号频率均为187.5Hz，信噪比均为-16dB，快拍数为5000。采用GFSSE对目标方位进行估计，DOA估计结果见图2。

从仿真结果中可以看出，当信号处于静止状态时，GFSSE通过将多个快拍进行结合可在一个较低的信噪比环境下实现高分辨性能，并且该方法的估计精度不受网格误差的影响，此时的估计误差为0.26°(若采用以1°为间隔的网格进行DOA估计，则估计误差至少为0.4°)。

Claims (1)

1.一种基于多项式求根的无网格稀疏谱估计方法，其特征在于估计步骤如下：

步骤1：采用阵元间距为半波长的M元均匀线列阵接收窄带信号,均匀线列阵上各个传感器将接收到的水声信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号x_i(n),1≤n≤N,i＝1,...,M；

将观测空间[-90°,90°]划分为Q个网格；所述90°为端射方向，各网格点所代表的方向角组成的向量记为Θ，Θ＝[θ₁,θ₂,...,θ_Q]；

在该网格上，阵列的接收信号模型表示为x(n)＝A(Θ)s(n)+e(n),n＝1,...,N

其中：

和

分别为各阵元上的接收信号、信号源以及各阵元上接收的环境噪声所组成的向量，上标“T”表示为转置符号；

为阵列流形矩阵，a(θ)＝[1e^{-j2πdsin(θ)f/c}...e^{-j2πd(M-1)sin(θ)f/c}]^T，f为窄带信号的中心频率，d为阵元间距，c为声速；

阵列接收信号的采样协方差矩阵为R＝E{x(n)x^H(n)}＝A(Θ)R_sA^H(Θ)+R_e

其中：E{·}为期望算子；R_s和R_e分别为信号协方差矩阵和噪声协方差矩阵；上标“H”为共轭转置符号；协方差矩阵由采样协方差矩阵

所代替；

对协方差矩阵两端进行向量化得

其中：

vec(·)为矩阵向量化算子；矩阵

的第i列为

上标“*”表示求共轭，

表示Kronecker积；p＝diag(R_s)为R_s对角线元素组成的向量，diag(·)表示由矩阵主对角元素组成的向量或以向量元素为主对角线的对角矩阵；p_e＝vec(R_e)；

保留矩阵

中元素不同的(2M-1)行，并将每列元素按照

的顺序进行排列，组成新的矩阵

采用同样的方式保留并重新排列r和p_e中的元素，得到新的向量r'和p'_e，则新的协方差矩阵模型表示为

步骤2：将p表示在连续空间上为

其中：t＝sinθ∈[-1,1]；t_i＝sinθ_i,i＝1,2,...,K为真实信号方位所对应的正弦值；p_i为相应的信号功率；δ(t-t_i)为在t_i上的狄拉克函数；

将协方差矩阵模型表示在连续空间上为

其中：p为连续空间上的信号功率；

表示为傅里叶变换算子；r'中第m个元素为

将稀疏谱估计的优化问题表示为：

||p′_e||₂≤ε

其中：

为原子范数，p的原子范数表示为

||·||₂为l₂范数；

trace(·)表示矩阵求迹；

构造优化问题的对偶问题为：

其中：

为关于

的对偶因子；Re[·]为求实部；

表示半正定符号；

步骤3、构造求根多项式：P(z)＝e^{j2πfd/ct2(M-1)}(1-|H(z)|²)

其中：z＝e^j2πfd/ct；

|·|为求模符号；

求解多项式的根，并选择多项式根中模值为1的根

计算目标信号的方位角：

其中：angle(·)表示求相位符号；

步骤4：计算出所有DOA估计值

重新构造维数为M×K的矩阵

所述

的第k列

得到估计的稀疏谱为：

其中：上标“+”为求Moore-Penrose逆。

Images