CN109188344A - 脉冲噪声环境下基于互循环相关music算法信源个数与来波方向角估计方法 - Google Patents

Abstract

本发明公开了脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，属于阵列信号处理领域。首先构建包括M个阵元和N个远场窄带信号的远场窄带信号源模型，计算第m个阵元t时刻的接收数据x_m(t)。然后针对t时刻，利用任意两个阵元的接收数据，计算这两个阵元的互循环协方差矩阵R(m,k)；将M×M个互循环相关协方差矩阵R(m,k)求和取均值，得到一个互循环相关协方差矩阵R_mutual。最后对互循环相关协方差矩阵R_mutual进行奇异值分解，并计算其功率谱密度P(θ)：搜索功率谱密度P(θ)的局部谱峰即得远场窄带信号的来波方向角估计值。本发明利用信号的循环平稳特性，互相关在自相关的基础上进一步地去掉随机噪声和杂波成分，提升了强脉冲噪声背景下来波方向角估计性能。

Description

脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波 方向角估计方法

技术领域

本发明属于阵列信号处理领域，具体是一种脉冲噪声环境下基于互循环相关MUSIC(多重信号分类，Multiple Signal Classification)算法信源个数与来波方向角估计方法。

背景技术

阵列信号处理领域的空间谱估计方向从子空间分解类算法的处理方式上看，可以分为两类，一是以MUSIC算法为代表的噪声子空间类算法，另一类是以旋转不变子空间(ESPRIT)为代表的信号子空间类算法。以MUSIC算法为代表的算法包括特征矢量法、MUSIC、求根MUSIC及MNM(最小范数，Minimum Norm)等；以ESPRIT为代表的算法主要有TAM(Toeplitz近似，Toeplitz Approximation)、LS-ESPRIT(最小二乘-旋转不变子空间，LeastSquares ESPRIT)及TLS-ESPRIT(总体最小二乘-旋转不变子空间，Total Least SquaresESPRIT)等。

实际监测环境中常见的是脉冲噪声而不是高斯噪声，脉冲噪声也称为非高斯噪声，在非高斯背景下，基于傅里叶分析的相关研究方法的效果将发生明显的下降。Gardner等人通过对人造信号的研究，发现大部分人造信号具有循环平稳特性，相同循环频率的信号循环相关，不同循环频率的信号循环互相关为零，将这一特性引入DOA(波束到达方向，Direction of Arrival)估计中，用阵列接收数据的循环自相关矩阵代替数据的自相关矩阵，通过信号子空间拟合等方法估计信号的方位。然而，自相关熵的方法引入了非常明显的低频噪声。

阵列信号的实际接收环境中，存在着诸如大气环境噪声、海杂波、地杂波、雷达反向散射回波及电磁噪声等大量的脉冲噪声。在非高斯背景下，以传统MUSIC算法为代表的噪声子空间相关研究方法的效果将发生明显的下降，针对此类噪声，可以利用接收信号的循环谱特性，消除干扰和背景噪声。用阵列接收数据的循环自相关矩阵代替数据的自相关矩阵，再进行SVD(奇异值分解，Singular Value Decomposition)分解与谱峰搜索，与相关协方差矩阵相比，阵列接收数据的循环自相关矩阵可以更好地抑制脉冲噪声，提高DOA估计算法的鲁棒性。

发明内容

为解决自相关引入的低频噪声影响准确性的问题，本发明提出了一种脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，在高脉冲噪声时依旧保持较高准确率，是一种高鲁棒性的角度估计算法。

具体步骤如下：

步骤一、构建包括M个阵元和N个远场窄带信号的远场窄带信号源模型；

所述信号源模型如下：

M个阵元组成均匀等距的阵列天线，阵元间距为d，N个远场窄带信号入射到该阵列上，各阵元接收到信号后经各自的传输通道送到处理器，处理器处理来自M个通道的数据。

步骤二、计算第m个阵元t时刻的接收数据x_m(t)；

输出的表达式为：

s_i(t)为第i个远场窄带信号的复包络；θ_i为第i个远场窄带信号入射到阵列天线上时的入射角度；λ为远场窄带信号的波长；n_m(t)表示服从α稳定分布的第m个阵元在t时刻的加性噪声；

步骤三、针对t时刻，利用任意两个阵元的接收数据，计算这两个阵元的互循环协方差矩阵R(m,k)；

针对第m个阵元与第k个阵元的接收数据，互循环协方差矩阵定义式如下：

R(m,k)＝<E[x_m(t)x_m(t+τ)](E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})])^H>_t

其中，m＝1,2,...,M，k＝1,2,...,M，<·>_t表示时间平均，(·)^H表示共轭转置，x_m(t+τ)表示第m个阵元t+τ时刻的接收信号，τ表示时延；x_k(t)表示第k个阵元t时刻的接收信号，ε为循环频率；E[x_m(t)x_m(t+τ)]是x_m(t)的自相关函数，E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})]是x_k(t)的循环自相关矩阵。

步骤四、将M×M个互循环相关协方差矩阵R(m,k)求和取均值，得到一个互循环相关协方差矩阵R_mutual；

步骤五、对互循环相关协方差矩阵R_mutual进行奇异值分解；

互循环相关协方差矩阵R_mutual大小为M×M，经过奇异值分解之后，U，∑和V均为M×M的方阵；U的元素向量称为左奇异向量，∑是实数对角矩阵，对角线外的元素都为0，对角线上的元素是按从大到小排列的奇异值；V的元素向量被称为右奇异向量；s代表信号子空间，n代表噪声子空间。U_s为信号子空间的左奇异向量方阵；U_n为噪声子空间的左奇异向量方阵；∑_s为信号子空间的实数对角矩阵；∑_n为噪声子空间的实数对角矩阵；V_s为信号子空间的右奇异向量方阵；V_n为噪声子空间的右奇异向量方阵；

步骤六、对奇异值分解后的互循环相关协方差矩阵R_mutual计算功率谱密度P(θ)：

a(θ)＝[1,exp(-j(2π/λ)dsinθ),...,exp(-j(2π/λ)(M-1)dsinθ)]^T；

步骤七、搜索功率谱密度P(θ)的局部谱峰即得远场窄带信号的来波方向角估计值。

本发明的优点在于：

1.一种脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，对于接收阵列数据，不同阵元先进行循环自相关运算，再计算阵元间的互相关，得到不同阵元间的协方差矩阵。

2.一种脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，对于运算得到的M×M组互循环相关协方差矩阵，算术融合后得到一个互循环相关协方差矩阵，再继续进行奇异值分解与谱峰搜索。

3、一种脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，在强脉冲噪声环境下，利用信号的循环平稳特性，互相关在自相关的基础上进一步地去掉随机噪声和杂波成分，进一步提升了强脉冲噪声背景下来波方向角估计性能。

附图说明

图1为本发明一种脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法流程图；

图2为本发明构建包括M个阵元接收N个远场窄带信号的示意图；

图3为本发明计算的搜索功率谱密度P(θ)的局部谱图；

图4为本发明在三种不同算法下广义信噪比GSNR与信源个数识别准确率Pd之间的对比图；

图5为本发明在三种不同算法下广义信噪比GSNR与均方根误差RMSE之间的对比图；

图6为本发明在三种不同算法下信源角度差与信源个数识别准确率Pd之间的对比图；

图7为本发明在三种不同算法下信源角度差与均方根误差RMSE之间的对比图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

当脉冲噪声进一步提高，阵列接收数据的广义信噪比GSNR进一步下降时，自相关引入的低频噪声会影响DOA估计的准确性能，本发明基于互相关理论与交织思想，交错阵元，计算接收阵列的多组互循环相关协方差矩阵，通过算术融合，将多组互循环相关协方差矩阵融合成一个互循环相关协方差矩阵，替代传统MUSIC算法中的协方差矩阵，再进行奇异值分解，把互循环相关协方差矩阵分解成信号子空间和噪声子空间，计算接收数据的功率谱，对功率谱进行谱峰搜索，局部峰值即为波达方向角的估计值。

如图1所示，具体步骤如下：

步骤一、构建包括M个阵元和N个远场窄带信号的远场窄带信号源模型；

所述信号源模型如下：

如图2所示，M个阵元组成均匀等距的阵列天线，阵元间距为d，N个远场频率已知的窄带信号入射到该阵列上，入射角度为θ₁,θ₂,...,θ_N，各阵元接收到信号后经各自的传输通道送到处理器，处理器处理来自M个通道的数据。

步骤二、以第一个阵元为参考阵元，计算第m个阵元t时刻的接收数据x_m(t)；

输出的表达式为：

si(t)为第i个远场窄带信号在t时刻入射到阵列的的复包络；θi为第i个远场窄带信号入射到阵列天线上时的入射角度；λ为远场窄带信号的波长；n_m(t)表示服从α稳定分布的第m个阵元在_t时刻的加性噪声；各阵元噪声相互独立，信号与噪声之间相互独立，(1)式写成矩阵形式：X＝AS+N。

远场窄带信号的复包络形式表示为：

式中，u_i(t)是第i个远场窄带信号在t时刻接收信号的幅度，ω_i是接收信号的频率，ω_i＝2πf_i，f_i表示第i个远场窄带信号的入射频率； 是t时刻接收信号的相位，j表示复数；τ表示时延；s_i(t-τ)表示第i个远场窄带信号在t-τ时刻入射到阵列的复包络。

假设：得到：

阵列天线第l个阵元接收信号的数学模型是：

其中，x_l(t)表示第l个阵元在t时刻的接收信号，a_li为第l个阵元对第i个远场窄带信号的导向矢量，n_l(t)表示服从α稳定分布的第l个阵元在t时刻的加性噪声；τ_li表示第i个远场窄带信号到达第l个阵元时相对于参考阵元的时延。

式(4)写成矢量形式，为：

X(t)＝A(θ)S(t)+N(t) (5)

式中，X(t)为阵列的M×1维快拍数据矢量，S(t)为空间信号的N×1维矢量，N(t)为阵列的M×1维噪声数据矢量，A(θ)为空间阵列的M×N维流型矩阵(导向矢量矩阵)，为：

A＝[a(θ₁),a(θ₂),...,a(θ_k),...,a(θ_N)] (6)

导向矢量a(θ_k)为：

a(θ_k)＝[1,exp(-j(2π/λ)dsinθ_k),...,exp(-j(2π/λ)(M-1)dsinθ_k)]^T (7)

α稳定分布是一类适用范围很宽并得到广泛应用的随机信号模型，包括高斯分布(α＝2)和分数低阶α稳定分布(0<α<2)两种情况，一般可用α稳定分布噪声表示脉冲噪声的数学模型，特征函数表示成以下形式：

ψ(t)是脉冲噪声随时间无规则出现的突发性干扰值；表示t和α的函数关系；sgn(t)表示符号函数：通过参数a，α，β和γ可确定一个稳定分布的特征函数，其中：(1)a∈(-∞,+∞)为位置参数，表示概率密度函数的位置，a的变化引起概率密度函数的水平移动；

(2)α∈(0,2]为特征参量，α越大，概率密度分布曲线越“矮胖”，拖尾越小；α越小，概率密度分布曲线越“高瘦”，拖尾越重；

(3)β∈[-1,1]为倾斜度参量，如果β>0，表示概率密度函数曲线向右倾斜，反之向左倾斜；

(4)γ∈[0,+∞)为分散性参数，表示数据集中程度，γ越大，数据以均值为中心的分散程度越大。

步骤三、针对t时刻，利用任意两个阵元的接收数据，计算这两个阵元的接收数据的自相关，得到互循环协方差矩阵R(m,k)；

针对第m个阵元与第k个阵元的接收数据，互循环相关协方差矩阵定义式如下：

R(m,k)＝<E[x_m(t)x_m(t+τ)](E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})])^H>_t (9)

其中，m＝1,2,...,M，k＝1,2,...,M，<·>_t表示时间平均，(·)^H表示共轭转置，x_m(t+τ)表示第m个阵元t+τ时刻的接收信号，τ表示时延；x_k(t)表示第k个阵元t时刻的接收信号，ε为循环频率；E[x_m(t)x_m(t+τ)]是x_m(t)的自相关函数，E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})]是x_k(t)的循环自相关矩阵。R(m,k)是M×M的方阵。

步骤四、将M×M个互循环相关协方差矩阵R(m,k)求和取均值，得到一个互循环相关协方差矩阵R_mutual；

步骤五、对互循环相关协方差矩阵R_mutual进行奇异值分解；

互循环相关协方差矩阵R_mutual大小为M×M，经过奇异值分解之后，U，∑和V均为M×M的方阵；U的元素向量称为左奇异向量，∑是实数对角矩阵，对角线外的元素都为0，对角线上的元素是按从大到小排列的奇异值；V的元素向量被称为右奇异向量；s代表信号子空间，n代表噪声子空间。U_s为信号子空间的左奇异向量方阵；U_n为噪声子空间的左奇异向量方阵；∑_s为信号子空间的实数对角矩阵，趋近于零；∑_n为噪声子空间的实数对角矩阵；V_s为信号子空间的右奇异向量方阵；V_n为噪声子空间的右奇异向量方阵；

步骤六、对奇异值分解后的互循环相关协方差矩阵R_mutual计算功率谱密度P(θ)：

a(θ)＝[1,exp(-j(2π/λ)dsinθ),...,exp(-j(2π/λ)(M-1)dsinθ)]^T；

步骤七、搜索功率谱密度P(θ)的局部谱峰即得远场窄带信号的来波方向角估计值。

如图3所示，纵坐标是功率谱密度P(θ)，横坐标是角度θ，搜索局部谱峰，谱峰个数是远场窄带信号的个数，谱峰对应的角度是估计出的来波方向角。

本发明提出的DOA估计算法，对比循环互相关MUSIC算法和经典MUSIC算法，根据信源个数识别准确率Pd和均方根误差RMSE两个指标，衡量本发明的算法性能。如下：

循环互相关MUSIC算法的协方差矩阵表示为：R_XX(ε,τ)＝<X(t)X^H(t+τ)e^-j2πεt>_t；

经典MUSIC算法的协方差矩阵表示为：R＝E[XX^H]＝AE[SS^H]A^H+σ²I＝AR_sA^H+σ²I；

R_s是信号相关矩阵，σ²I是噪声相关矩阵，σ²是噪声功率，I是M×M阶的单位矩阵。

(1)信源个数识别准确率Pd的定义式如下：

(2)均方根误差RMSE的定义式如下：

其中，y_i表示通过谱峰搜索出的入射方向角，i＝1,2,...,N，θ_i表示远场窄带信号实际入射到阵列的来波方向角。

(3)脉冲噪声下，广义信噪比GSNR等效为高斯噪声下的信噪比SNR，定义为：其中，表示信号的方差，γ表示脉冲噪声的分散性参数，与脉冲噪声的集中程度有关系。

循环互相关MUSIC算法用cyclic曲线表示；经典MUSIC算法用music曲线表示；本发明提出的DOA估计算法性能用mutual曲线表示；当快拍数是200或400时，广义信噪比GSNR与信源个数识别准确率Pd的关系如图4所示，本发明提出的方案在GSNR＝8dB时，准确率高达93％，两种快拍数下，本发明提出的方案检测准确率均高于其他两种算法。当快拍数是200或400时，广义信噪比GSNR与均方根误差RMSE之间的对比如图5所示，两种快拍数下，本发明提出的方案在GSNR<0时，均方根误差RMSE均低于其他两种算法，低GSNR下角度估计相对准确。

在广义信噪比GSNR＝0时，信源角度差与信源个数识别准确率Pd的关系如图6所示。从图6可以看出，本发明提出的方案，随着信源角度差的增大，检测准确率均高于其他两种算法。

信源角度差与均方根误差RMSE之间的关系如图7所示，随着信源角度差的增大，三种算法的RMSE均在减小，本发明提出的方案，整体性能优于其他两种算法。

仿真结果表明，本发明提出的DOA估计算法，与循环互相关MUSIC算法和经典MUSIC算法比较，高脉冲噪声时信源个数和来波方向估计保持较高准确率，是一种高鲁棒性的角度估计算法。

Claims (2)

1.脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，其特征在于，具体步骤如下：

步骤一、构建包括M个阵元和N个远场窄带信号的远场窄带信号源模型；

步骤二、计算第m个阵元t时刻的接收数据x_m(t)；

输出的表达式为：

s_i(t)为第i个远场窄带信号的复包络；θ_i为第i个远场窄带信号入射到阵列天线上时的入射角度；λ为远场窄带信号的波长；n_m(t)表示服从α稳定分布的第m个阵元在t时刻的加性噪声；

步骤三、针对t时刻，利用任意两个阵元的接收数据，计算这两个阵元的互循环协方差矩阵R(m,k)；

针对第m个阵元与第k个阵元的接收数据，互循环协方差矩阵定义式如下：

R(m,k)＝＜E[x_m(t)x_m(t+τ)](E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})])^H＞_t

其中，m＝1,2,...,M，k＝1,2,...,M，＜·＞_t表示时间平均，(·)^H表示共轭转置，x_m(t+τ)表示第m个阵元t+τ时刻的接收信号，τ表示时延；x_k(t)表示第k个阵元t时刻的接收信号，ε为循环频率；E[x_m(t)x_m(t+τ)]是x_m(t)的自相关函数，E[(x_k(t)e^-j2πεt)(x_k(t+τ)e^{-j2πε(t+τ)})]是x_k(t)的循环自相关矩阵。

步骤四、将M×M个互循环相关协方差矩阵R(m,k)求和取均值，得到一个互循环相关协方差矩阵R_mutual；

步骤五、对互循环相关协方差矩阵R_mutual进行奇异值分解；

互循环相关协方差矩阵R_mutual大小为M×M，经过奇异值分解之后，U，∑和V均为M×M的方阵；U的元素向量称为左奇异向量，∑是实数对角矩阵，对角线外的元素都为0，对角线上的元素是按从大到小排列的奇异值；V的元素向量被称为右奇异向量；s代表信号子空间，n代表噪声子空间。U_s为信号子空间的左奇异向量方阵；U_n为噪声子空间的左奇异向量方阵；∑_s为信号子空间的实数对角矩阵；∑_n为噪声子空间的实数对角矩阵；V_s为信号子空间的右奇异向量方阵；V_n为噪声子空间的右奇异向量方阵；

步骤六、对奇异值分解后的互循环相关协方差矩阵R_mutual计算功率谱密度P(θ)：

a(θ)＝[1,exp(-j(2π/λ)dsinθ),...,exp(-j(2π/λ)(M-1)dsinθ)]^T；

步骤七、搜索功率谱密度P(θ)的局部谱峰即得远场窄带信号的来波方向角估计值。

2.如权利要求1所述的脉冲噪声环境下基于互循环相关MUSIC算法信源个数与来波方向角估计方法，其特征在于，所述信号源模型如下：

M个阵元组成均匀等距的阵列天线，阵元间距为d，N个远场窄带信号入射到该阵列上，各阵元接收到信号后经各自的传输通道送到处理器，处理器处理来自M个通道的数据。