CN108919189A - 一种频率和方位联合估计的阵列信号处理方法 - Google Patents

Abstract

本发明涉及一种频率和方位联合估计的阵列信号处理方法，首先在若干个假设频率的条件下，针对每一个假设的单频信号，建立了状态方程和观测方程。在此基础上，利用卡尔曼滤波器提高了信号的信噪比，进而再利用波束形成器获得目标方位。然后，建立与假设频率和目标方位有关的目标函数。最后，通过使得目标函数达到最大而获得最终的频率‑DOA联合估计。本发明使得频率估计和目标方位估计的均方根误差减小，实现低信噪比环境中频率‑DOA联合估计，提高了接收设备的作用距离，克服了现有方法在低信噪比下参数估计性能下降的问题。

Description

一种频率和方位联合估计的阵列信号处理方法

技术领域

本发明属于阵列信号处理领域，涉及一种频率和方位联合估计的阵列信号处理方法。

背景技术

参数估计是阵列信号处理中的主要任务之一，包括频率估计、目标方位估计(Direction of arrival，DOA)、频率-DOA联合估计等。比如在舰船和潜艇的辐射噪声中，存在单频线谱信号，因此它可以用于检测目标。但是在接收阵列工作环境中，当接收阵列与目标的距离较远时，则阵列接收的信噪比很低，工作距离有限。

频率估计是参数估计中的一个重要任务。低信噪比将导致传统的一些频率估计的方法性能下降，包括傅里叶变换(参见：Frequency estimation using wraped discreteFourier transform.Signal Processing,2003,83(8):1661-1671)、基于特征分析的频率估计方法(参见：Analysis of MUSIC and ESPRIT frequency estimates for sinusoidalsignals with lowpass envelopes.IEEE Transactions on Signal Processing,1996,44(9):2359-2364)以及时频分析方法(参见：Use of the Cross Wigner-VilleDistribution for Estimation of Instantaneous Frequency.IEEE Transactions onSignal Processing,1993,41(3):1439-1445)等等。

目标方位估计是阵列信号处理中的另一个主要任务。通常，利用子空间分解方法和波束形成方法来估计目标方位。基于子空间分解方法的算法最大的缺点是在低信噪比下，信号源数估计不正确，其性能可能会严重下降。波束形成方法包括延时求和波束形成方法(Delay-and-sum Beamforming，DAS)和最小方差无畸变响应波束形成方法(minimumvariance distortionless response,MVDR)等。当快拍数较少或者阵列存在位置误差和通道相幅误差时，MVDR波束形成方法的稳健性较差。相比之下，DAS波束形成方法具有较好的稳健性，得到了较为广泛的应用，但是该方法的噪声抑制能力受到了阵列孔径的限制。为了提高DAS波束形成方法的噪声抑制能力，进而提高目标方位估计的性能，国内外的学者研究了基于DAS波束形成方法的各种改进方法。有学者将协方差矩阵对角线置0，改善了DAS的性能(参见：Beamforming in Acoustic Testing.Berlin:2002:83-86.)。有学者提出了复杂噪声场条件下的对角减载技术，大大提高了噪声抑制的效果(参见：复杂噪声场下对角减载技术的原理及应用.物理学报,2017；66(1):152-161)。上述这些基于DAS波束形成方法的各种改进方法均是对协方差矩阵进行处理，没有充分挖掘噪声抑制能力。在低信噪比下，基于上述这些方法的方位估计性能下降。

在一些实际应用中，需要同时估计一些参数。ESPRIT方法(参见：Total leastsquares phased averaging and 3-D ESPRIT for joint azimuth-elevation-carrierestimation.IEEE Transactions on Signal Processing,2001,49(1):54–62)和最大似然方法(参见：Maximum likelihood angle-frequency estimation in partially knowncorrelated noise for low-elevation targets.IEEE Transactions on SignalProcessing,2005,53(8):3057–3064)是两种常见的频率-DOA联合估计方法。最大似然方法呈现出了较好的性能，但是其计算复杂度较高，需要进行多维搜索。相比之下，ESPRIT的计算复杂度较低，但是在低信噪比下，ESPRIT方法的性能下降。

因此，在低信噪比环境中实现参数估计是一个极具挑战而又亟待解决的问题。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种频率和方位联合估计的阵列信号处理方法，实现低信噪比下频率估计和方位估计。

技术方案

一种频率和方位联合估计的阵列信号处理方法，其特征在于步骤如下：

步骤1：以M个水听器组成的接收阵列接收来自空间的信号和噪声数据，表示为x(k)＝s(k)+n(k)，其中s(k)＝[s₁(k),s₂(k),…,s_M(k)]^T为阵列接收到的信号，第m号阵元接收到信号表示为s_m(k)＝A_msin[2πf(k-1)/f_s+φ_m]+v_m(k)，A_m为幅度，φ_m为相位，f_s为采样频率，v_m(k)表示信号畸变，是具有独立样本的高斯随机过程。n(k)＝[n₁(k),n₂(k),…,n_M(k)]^T为阵列接收到的高斯白噪声，协方差矩阵为Q_n；

步骤2：构建状态方程

所述B(f)为状态矩阵

所述K_Ts(k)为状态向量

其中

式中，表示Hilbert变换；

所述v(k)为信号畸变矢量，表示为其协方差矩阵为Q；

所述为状态方程中的驱动噪声矢量，表示为v′(k+1)＝v(k+1)-B(f)v(k)，其协

方差矩阵为Q′＝2Q；

步骤3：构建观测方程，表示为：

步骤4：将f_low≤f≤f_high频率范围按照等频率间隔，分成L个假设频率，记为f_l,l＝1,2,…,L，在每一个假设频率下计算得到状态矩阵B(f_l)；并根据步骤2和步骤3获得的状态方程和观测方程，利用Kalman滤波器获得每一个假设频率下的输出，具体递推式顺序计算如下：

1.初始化和(-1|-1)

2.预测：

3.最小预测MSE(Mean Squared Error)矩阵：

4.卡尔曼增益矩阵：K(k)＝M(k|k-1)D^T[Q_n+DM(k|k-1)D^T]^-1

5.修正：

6.最小MSE矩阵：M(k|k)＝[I-K(k)D]M(k|k-1)，I为单位矩阵；

重复2～6得到Kalman滤波器的输出将该输出乘以矩阵D，得到阵列的输出并利用Hilbert变换求取的解析信号形式，记为Z(k,f_l)；

步骤5：利用步骤4中的输出信号Z(k,f_l)，计算得到方位谱

P_l(θ)＝w^H(f_l,θ)R(f_l)w(f_l,θ)

R(f_l)为Z(k,f_l)的协方差矩阵，w(f_l,θ)为波束形成器的加权向量；上标H表示复共轭转置；

计算目标方位估计：

得到目标估计方位所在主瓣外的最大旁瓣为最大旁瓣级最大旁瓣级SLL_l；

步骤6：将阵列接收到的数据分成T段，对每一段数据，利用步骤4～步骤5计算得到最大旁瓣级SLL_l,t，t＝1,2,…,T和目标方位估计t＝1,2,…,T；

对每一个固定的l，计算得到的方差为同时计算得到平均最大旁瓣级为：

建立目标函数为：

通过使得目标函数最大化获得频率和方位的联合估计，表示为：

其中l_op是使得F(f_l)达到最大的f_l的下标，为频率估计，为方位估计。

有益效果

本发明提出的一种频率和方位联合估计的阵列信号处理方法，基于Kalman滤波器和波束形成器，提出了基于Kalman滤波器和波束形成器的频率-DOA联合估计方法，称为minimum variance and sidelobe level method(MVS)。首先在若干个假设频率的条件下，针对每一个假设的单频信号，建立了状态方程和观测方程。在此基础上，利用卡尔曼滤波器提高了信号的信噪比，进而再利用波束形成器获得目标方位。然后，建立与假设频率和目标方位有关的目标函数。最后，通过使得目标函数达到最大而获得最终的频率-DOA联合估计。本发明使得频率估计和目标方位估计的均方根误差减小，实现低信噪比环境中频率-DOA联合估计，提高了接收设备的作用距离，克服了现有方法在低信噪比下参数估计性能下降的问题。

本发明采建立了单频信号的状态方程和观测方程，并用Kalman滤波器增强了单频信号，提高了信噪比，然后利用本发明所提的MVS方法，获得了频率和方位的联合估计。本发明使得频率估计和目标方位估计的均方根误差减小，实现低信噪比环境中频率-DOA联合估计，提高了设备的作用距离，克服了现有方法在低信噪比下目标方位估计性能下降的问题。

附图说明

图1：参数估计均方根误差随输入信噪比的变化情况之方位估计的均方根误差；

图2：参数估计均方根误差随输入信噪比的变化情况之频率估计的均方根误差

具体实施方式

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

下面详细介绍本发明的具体操作步骤：

步骤1：20个水听器组成的接收阵列，阵元间距为0.75米，接收来自空间的信号和噪声数据，表示为x(k)＝s(k)+n(k)，其中s(k)＝[s₁(k),s₂(k),…,s_M(k)]^T为阵列接收到的信号，第m号阵元接收到信号表示为s_m(k)＝A_msin[2πf(k-1)/f_s+φ_m]+v_m(k)，A_m为4v，φ_m为0rad，f_s为8000Hz，f为1000Hz，v_m(k)的方差为0.001。n(k)＝[n₁(k),n₂(k),…,n_M(k)]^T为阵列接收到的高斯白噪声，协方差矩阵为Q_n。接收数据时间长度为12.5s，采样点数为100000；

步骤2：构建状态方程，表示为：

K_Ts(k+1)＝B(f)[K_Ts(k)-v(k)]+v(k+1)

＝B(f)K_Ts(k)+v′(k+1)

依次说明上式中各变量的计算公式，

1.状态矩阵B表示为：

2.K_Ts(k)为状态向量，表示为：

其中

式中，表示Hilbert变换；

3.v(k)为信号畸变矢量，表示为其协方差矩阵为Q；

4.为状态方程中的驱动噪声矢量，表示为v′(k+1)＝v(k+1)-B(f)v(k)，其协方差矩阵为Q′＝2Q。

步骤3：根据步骤1构建观测方程，表示为：

x(k)＝DK_Ts(k)+n(k)

其中，观测矩阵D表示为：

步骤4：将感兴趣的频率范围(700≤f≤1300Hz)按等频率间隔，分成301个假设频率，记为f_l,l＝1,2,…,301，在每一个假设频率下计算得到状态矩阵B(f_l)。并根据步骤2和步骤3获得的状态方程和观测方程，利用Kalman滤波器获得每一个假设频率下的输出，具体递推式顺序计算如下：

1.初始化

2.预测：

3.最小预测MSE(Mean Squared Error)矩阵：M(k|k-1)＝B(f_l)M(k-1|k-1)B^T(f_l)+Q′

4.卡尔曼增益矩阵：K(k)＝M(k|k-1)D^T[Q_n+DM(k|k-1)D^T]^-1

5.修正：

6.最小MSE矩阵：M(k|k)＝[I-K(k)D]M(k|k-1)，I为单位矩阵。

重复2～6即可得到Kalman滤波器的输出将该输出乘以观测矩阵D，得到阵列的输出并利用Hilbert变换求取的解析信号形式，记为Z(k,f_l)；

步骤5：利用步骤4中的输出信号Z(k,f_l)，计算得到方位谱

P_l(θ)＝w^H(f_l,θ)R(f_l)w(f_l,θ)，

R(f_l)为Z(k,f_l)的协方差矩阵，w(f_l,θ)为波束形成器的加权向量，上标H表示复共轭转置。进而计算得到目标方位估计和最大旁瓣级SLL_l(方位谱中，除目标估计方位所在主瓣外的最大旁瓣为最大旁瓣级)；

步骤6：将阵列接收到的数据分成100段，每一段时长0.125s，采样点数为1000。对每一段数据，利用步骤4～步骤5计算得到最大旁瓣级SLL_l,t，t＝1,2,…,100和目标方位估计t＝1,2,…,100。进而，对每一个固定的l，计算得到的方差为同时计算得到平均最大旁瓣级为：

建立目标函数F(f_l)为：

通过使得目标函数最大化获得频率和方位的联合估计，表示为：

其中l_op是使得F(f_l)达到最大的f_l的下标，为频率估计，为方位估计。

上述过程可以获得频率和方位的联合估计，为了评估MVS方法在频率-DOA联合估计中的性能，采用蒙特卡洛实验的方法，实验次数为100次，频率估计和方位估计的均方根误差(RMSE)定义为：

其中，θ_s为真实目标方位，和分别为一次仿真实验的方位估计和频率估计。得到方位估计的RMSE和频率估计的RMSE如图1所示。为了对比方位估计的性能，DAS方法的方位估计结果是在频率已知的条件下计算得到的。可以发现，当信噪比为-30dB时，MVS方法的方位估计和频率估计的均方根误差依然很小。相比之下，ESPRIT方法的频率和方位估计的均方根误差在低信噪比下较大。此外，DAS方法的方位估计的均方根误差在信噪比低于-20dB时便开始增大。充分说明了本发明提出的MVS方法在低信噪比下的有效性和优越性。

Claims (1)

1.一种频率和方位联合估计的阵列信号处理方法，其特征在于步骤如下：

步骤1：以M个水听器组成的接收阵列接收来自空间的信号和噪声数据，表示为x(k)＝s(k)+n(k)，其中s(k)＝[s₁(k),s₂(k),…,s_M(k)]^T为阵列接收到的信号，第m号阵元接收到信号表示为s_m(k)＝A_msin[2πf(k-1)/f_s+φ_m]+v_m(k)，A_m为幅度，φ_m为相位，f_s为采样频率，v_m(k)表示信号畸变，是具有独立样本的高斯随机过程。n(k)＝[n₁(k),n₂(k),…,n_M(k)]^T为阵列接收到的高斯白噪声，协方差矩阵为Q_n；

步骤2：构建状态方程

所述B(f)为状态矩阵所述K_Ts(k)为状态向量

其中

式中，表示Hilbert变换；

所述v(k)为信号畸变矢量，表示为其协方差矩阵为Q；所述为状态方程中的驱动噪声矢量，表示为v′(k+1)＝v(k+1)-B(f)v(k)，其协方差矩阵为Q′＝2Q；

步骤3：构建观测方程，表示为：

步骤4：将f_low≤f≤f_high频率范围按照等频率间隔，分成L个假设频率，记为fl,l＝1,2,…,L，在每一个假设频率下计算得到状态矩阵B(f_l)；并根据步骤2和步骤3获得的状态方程和观测方程，利用Kalman滤波器获得每一个假设频率下的输出，具体递推式顺序计算如下：

1.初始化和M(-1|-1)

2.预测：

3.最小预测MSE(Mean Squared Error)矩阵：

M(k|k-1)＝B(f_l)M(k-1|k-1)B^T(f_l)+Q′

4.卡尔曼增益矩阵：K(k)＝M(k|k-1)D^T[Q_n+DM(k|k-1)D^T]^-1

5.修正：

6.最小MSE矩阵：I为单位矩阵；

重复2～6得到Kalman滤波器的输出将该输出乘以矩阵D，得到阵列的输出并利用Hilbert变换求取的解析信号形式，记为Z(k,f_l)；

步骤5：利用步骤4中的输出信号Z(k,f_l)，计算得到方位谱

P_l(θ)＝w^H(f_l,θ)R(f_l)w(f_l,θ)

R(f_l)为Z(k,f_l)的协方差矩阵，w(f_l,θ)为波束形成器的加权向量；上标H表示复共轭转置；

计算目标方位估计：

得到目标估计方位所在主瓣外的最大旁瓣为最大旁瓣级最大旁瓣级SLL_l；

步骤6：将阵列接收到的数据分成T段，对每一段数据，利用步骤4～步骤5计算得到最大旁瓣级SLL_l,t，t＝1,2,…,T和目标方位估计

对每一个固定的l，计算得到的方差为同时计算得到平均最大旁瓣级为：

建立目标函数为：

通过使得目标函数最大化获得频率和方位的联合估计，表示为：

其中l_op是使得F(f_l)达到最大的f_l的下标，为频率估计，为方位估计。