CN107390187B - 基于混合高斯分布的杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于混合高斯分布的杂波抑制方法，主要解决现有技术对频谱较宽的杂波抑制效果不好，对目标周围杂波抑制不掉的问题。其方案是：1)接收回波数据；2)对回波数据做动目标显示处理；3)对动目标显示处理后的数据做二维傅里叶变换；4)统计二维傅里叶变换后的纯杂波数据获得杂波混合高斯模型；5)用混合高斯对消法对二维傅里叶变换后的含有目标和杂波的数据做杂波抑制处理；6)对杂波抑制处理后的数据做质心凝聚处理得到目标数据矩阵。本发明能抑制掉任何多普勒分布的杂波，具有通用性好和精确性高的优点，可用于雷达回波信号的目标检测。

Description

基于混合高斯分布的杂波抑制方法

技术领域

本发明属于通信技术领域，更进一步涉及一种杂波抑制方法，可用于雷达回波信号的目标检测。

背景技术

传统的目标检测方法，通过使用动目标显示技术和自适应滤波器抑制杂波，然后利用恒虚警技术检测目标，但是动目标显示技术只能抑制零频附近的杂波，并且目标速度较低时，动目标显示技术会同时抑制目标，而自适应滤波器也只能抑制特定频率的杂波，使目标检测性能下降。

胡可欣，胡爱明在其发表的论文“自适应杂波抑制技术在雷达中的应用”(现代电子技术，2006年第8期总第223期24～26页)中提出一种杂波抑制方法，针对雷达杂波抑制中常用的级联动目标显示滤波器，提出采用时变加权原理，通过对动杂波速度估值，实现对动目标杂波抑制滤波器权值的优化，最后给出了这两种杂波抑制级联型动目标显示的工程实现。该方法的不足之处是不能消除多普勒展宽较宽的全部杂波。

温博、郑远、刘兵在其发表的论文“基于零频抑制与杂波图的地面慢速目标检测”(电子设计工程，2014年第22卷第10期85～87页)中提出一种检测地面慢速目标的方法。采用Kalmus滤波器提高零频杂波抑制能力，采用空间邻域插值方式完成目标区域的杂波背景积累，该方法虽说可解决慢速目标检测不连续的问题，但却存在目标速度较大时目标会被滤波器抑制，从而检测不到目标的不足。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于混合高斯分布的杂波抑制方法。以充分利用雷达回波数据中目标和杂波幅值不同的分布特性，在没有目标出现时统计出杂波幅值分布特性，在有目标时通过杂波的分布特性抑制杂波，实现对目标的准确检测。

实现本发明目的的技术方案是：首先将雷达回波数据做动目标显示处理，并分别在距离和多普勒维做傅里叶变换得到目标和杂波的距离、速度和幅值的三维信息，然后通过对多组纯杂波数据进行统计分析，得出杂波幅值分布服从混合高斯分布特性，判断含有目标和杂波的数据，数据中符合此分布特性的点视为杂波，剩余不符合此分布特性的为目标，具体步骤包括如下：

(1)采集雷达接收到的回波数据X，并以矩阵的形式进行存储；

(2)通过动目标显示MTI滤波器滤除回波数据中的零频静物杂波；

(3)对滤波器输出的数据分别进行距离维和多普勒维傅里叶变换，得到含有目标、杂波的距离和多普勒信息的数据矩阵Z；

(4)根据数据矩阵Z中的纯杂波数据的分布特性，获得一个表示杂波的混合高斯分布模型：

(4a)统计数据矩阵Z中的多组纯杂波数据的幅值得到单高斯分布的方差

并计算数据矩阵Z中第一组杂波数据的平均幅值μ₀；

(4b)将四阶混合高斯的方差初始化为

将第一阶高斯分布的均值初始化为μ₀，将其他三阶高斯分布的均值初始化为零，初始化四阶高斯分布的权系数ω₀和学习因子ρ₀；

(4c)每次从数据矩阵Z中读取一组杂波数据，对每组数据中的每个点的值进行判断：

若符合在上一次更新得到的混合高斯分布中的任意一阶高斯分布，则按如下公式增加此阶高斯分布的权重：

y_i,t＝(1-α)·y_i,t-1+α

否则，按如下公式减小此阶高斯分布的权重：

y_i,t＝(1-α)·y_i,t-1

其中，y_i,t-1为读取第t组数据时的第i阶高斯分布，y_i,t为读取第t组数据更新之后的第i阶高斯分布，α为权值因子；

(4d)按如下公式更新此阶高斯分布的均值和方差：

μ_i,t＝(1-ρ)·μ_i,t-1+ρX_t

其中，μ_i,t-1为读取第t组数据时的第i阶高斯分布的均值，μ_i,t为读取第t组数据更新之后的第i阶高斯分布的均值，X_t为读取的第t组数据的每一点的幅值，ρ为学习因子，

为读取第t组数据时的第i阶高斯分布的方差，

为读取第t组数据更新之后的第i阶高斯分布的方差，T表示转置；

(4e)按如下公式更新杂波的混合高斯分布模型：

其中，P(X_t)为读取第t组数据时的混合高斯分布模型，X_t为读取的第t组数据的每一点的幅值，ω_i,t为第t组数据第i阶高斯分布所占权值，y(X_t,μ_i,t,σ_i,t)为第i阶高斯分布，μ_i,t为第t组数据第i阶高斯分布的均值，σ_i,t为第t组数据第i阶高斯分布的协方差，μ_t为第t组数据更新之后的混合高斯分布模型的均值，σ_t为第t组数据更新之后的混合高斯分布模型的协方差；

(4f)重复步骤(4c)～(4e)，更新得到稳定的杂波的混合高斯分布模型；

(5)根据数据矩阵Z获得只含有目标的数据矩阵G：

(5a)依次取数据矩阵Z中含目标和杂波的每组数据中每一点k(i,j)的幅值z_i,j，按如下公式进行混合高斯判断处理得到杂波抑制后输出的数据矩阵A：

其中，A_i,j为数据矩阵A中第i行第j列的点；

(5b)对数据矩阵A做质心凝聚处理获得的目标数据矩阵G。

本发明与现有技术相比具有以下优点：

第一，本发明采用混合高斯分布杂波抑制的方法，可以有效地统计出所有杂波幅值分布特性，从而抑制掉任何距离和多普勒分布的杂波，克服了现有技术只能抑制特定范围多普勒杂波的缺点，使得本发明具有通用性高的优点。

第二，本发明采用逐点判断是否为杂波的方法，克服了现有技术不能对消掉目标周围杂波的缺点，使得本发明具有精确性高，杂波抑制效果好的优点。

附图说明

图1为本发明的实现流程图；

图2为本发明的实验结果图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1.接收回波数据。

连续采集雷达阵元接收到的多组回波数据，并将多组回波数据保存到不同的txt文件中；

用fopen函数依次读取多个txt文件中的回波数据，并保存到回波数据矩阵X中。

步骤2.动目标显示处理。

通过动目标显示MTI滤波器滤除回波数据中零频静物杂波，得到输出数据Y；

Y＝W*X，

其中，W表示MTI滤波器的响应函数，X表示输入的回波数据矩阵，*表示卷积操作。

步骤3.二维傅里叶变换。

按照如下公式对滤波器输出的数据Y进行二维傅里叶变换，得到含有目标、杂波的距离和多普勒信息的数据矩阵Z；

其中，Y(m,n)是步骤2中MTI滤波器的输出数据，m,n分别表示数据中的快时间和慢时间信息，Z(d,r)是傅里叶变换之后的数据矩阵，d表示目标与杂波的距离信息，r表示目标与杂波的多普勒信息，D表示目标与杂波的最大距离，R表示目标与杂波的最大多普勒。

步骤4.获得杂波混合高斯模型。

(4a)统计数据矩阵Z中的多组纯杂波数据的幅值，得到单高斯分布的方差

计算数据矩阵Z中第一组杂波数据的均值μ₀；

(4b)取混合高斯模型的阶数为四阶，并进行如下初始化：

将四阶混合高斯模型的第一阶高斯分布的方差

第二阶高斯分布的方差

第三阶高斯分布的方差

第四阶高斯分布的方差

都初始化为

将四阶混合高斯模型的第一阶高斯分布的均值μ₁初始化为μ₀，将四阶混合高斯模型的第二阶高斯分布的均值μ₂、第三阶高斯分布的均值μ₃、第四阶高斯分布的均值μ₄都初始化为零；

(4c)初始化第i阶高斯分布函数：

根据上面初始化的第i阶高斯分布的方差

均值μ_i获得初始化的第i阶高斯分布函数：

其中，exp(·)是指数函数，σ_i是第i阶高斯分布的均方差，

i＝1,2,3,4，x表示服从此阶高斯分布的随机变量；

(4d)更新第i阶高斯分布函数：

从数据矩阵Z中读取第t组纯杂波数据，按照如下公式判断出每组数据中的每个点的幅值z_k是否属于第i阶高斯分布：

若i＝1,2,3,4，则读取的数据中的点的幅值z_k属于第i阶高斯分布，并进行如下操作：

增加第i阶高斯分布的权重：y_i,t＝(1-α_i)·y_i,t-1+α_i

减小第j阶高斯分布的权重：y_j,t＝(1-α_j)·y_j,t-1，

更新第i阶高斯分布的均值：μ_i,t＝(1-ρ_i)·μ_i,t-1+ρ_iX_t，

更新第i阶高斯分布的和方差：

其中，y_i,t-1为读取第t组数据时的第i阶高斯分布，y_i,t为读取第t组数据更新之后的第i阶高斯分布，α_i为第i阶高斯分布的权值因子，y_j,t-1为读取第t组数据时的第j阶高斯分布，y_j,t读取第t组数据更新之后的第j阶高斯分布，α_j为第j阶高斯分布的权值因子，j≠i并且j＝1,2,3,4，μ_i,t-1为读取第t组数据时的第i阶高斯分布的均值，μ_i,t为读取第t组数据更新之后的第i阶高斯分布的均值，X_t为读取的第t组数据的每一点的幅值，ρ_i为第i阶高斯分布的学习因子，

为读取第t组数据时的第i阶高斯分布的方差，

为读取第t组数据更新之后的第i阶高斯分布的方差，T表示转置；

若i＝5，则读取的数据中的点的幅值z_k不属于第i阶高斯分布，将第四阶高斯分布的参数更新为：μ₄＝z_k，

(4e)根据(4d)更新的第i阶高斯分布，更新杂波的混合高斯分布模型：

其中，P(X_t,μ_t,σ_t)为第t组数据更新之后的混合高斯分布模型；X_t为读取的第t组数据的每一点的幅值；ω_i,t为读取第t组数据时第i阶高斯分布所占权值；y(X_t,μ_i,t,σ_i,t)为读取第t组数据更新之后的第i阶高斯分布；μ_i,t为读取第t组数据更新之后的第i阶高斯分布的均值；σ_i,t为读取第t组数据更新之后的第i阶高斯分布的协方差；μ_t为第_t组数据更新之后的混合高斯分布模型的均值；σ_t为第_t组数据更新之后的混合高斯分布模型的协方差；

(4f)设数据矩阵Z中纯杂波数据的总组数为T_c，重复步骤(4d)～(4e)，直到步骤(4d)中从数据矩阵Z中读取第t组纯杂波数据中的t＝T_c，获得稳定的混合高斯分布P(x_p,μ,σ)：

其中，y(x_i,μ_i,σ_i)是最后得到的第i阶高斯分布，x_i是符合第i阶高斯分布的随机变量，μ_i是第i阶高斯分布的均值，σ_i是第i阶高斯分布的均方差，ω_i是第i阶高斯分布所占权值，x_p表示服从此混合高斯分布的随机变量，μ是此混合高斯分布的均值，σ是此混合高斯分布的均方差。

步骤5.杂波抑制处理。

现有的杂波抑制处理方法有：杂波图法、动目标显示技术、动目标检测技术和脉冲多普勒技术等多种方法，但这些方法对于频谱较宽的杂波和目标周围的杂波抑制效果不好，本发明采用的混合高斯对消法可以达到很好的抑制效果，具体步骤如下：

依次取数据矩阵Z中含目标和杂波的每组数据中每一点k(i,j)的幅值z_i,j，进行混合高斯处理，得到杂波抑制后输出的数据矩阵A：

其中，M是数据矩阵A中总的行数，Q是数据矩阵A中总的列数，数据矩阵A中第i行第j列的点数据为：

步骤6.质心凝聚处理。

(6a)依次取数据矩阵A中相邻的距离区域

其中A_i,j表示数据矩阵A中第i行第j列的点数据，A_i+1,j表示数据矩阵A第i+1行第j列的点数据，i＝1,2,3,…M，j＝1,2,3,…,Q；

(6b)依次取数据矩阵A中相邻的多普勒区域[A_i,j … A_i,j+4]，其中，A_j+4表示数据矩阵A第i行第j+4列的点数据；

(6c)根据以上取得的距离区域和多普勒区域，得到距离多普勒矩阵区域

(6d)对距离多普勒矩阵区域B进行质心凝聚处理，得到目标数据矩阵G：

其中，目标数据矩阵G中第i行第j列的点数据G_i,j表示如下：

下面结合实验对本发明的效果做进一步的描述。

1.实验条件：

本发明实验的运行平台的配置如下：

CPU：Intel(R)Core(TM)i7-4790CPU@3.60GHz，内存8.00GB；

操作系统：Windows 7旗舰版64位操作系统；

实验软件：MATLAB R(2014b)。

本发明的实验参数设置如表1所示：

表1 信号参数及实验参数一览表

参数	取值
		目标距离R<sub>0</sub>	100m
目标速度v<sub>0</sub>	0.95m/s
		第i阶高斯分布的权值因子α<sub>i</sub>	5
第i阶高斯分布的学习因子ρ<sub>i</sub>	3
		第一阶高斯分布的权系数初值ω<sub>1</sub>	8
第二阶高斯分布的权系数初值ω<sub>2</sub>	2
		第三阶高斯分布的权系数初值ω<sub>3</sub>	0.5
第四阶高斯分布的权系数初值ω<sub>4</sub>	0.4

2.实验内容：

对接收到的回波数据进行动目标显示和二维傅里叶变换处理，得到混合高斯对消之前的数据矩阵，并用MATLAB画出此数据矩阵的距离速度二维示意图，并对其杂波抑制。结果如图2，其中：

图2(a)是混合高斯对消前的距离速度图，可以看出有较多的杂波；

图2(b)为对图2(a)混合高斯对消后的结果，其中x轴表示速度维，y轴表示距离维，z轴表示幅度维；

图2(c)是对图2(b)进行质心凝聚处理之后的距离速度图，剩余的一点为目标点。

从图2(b)中可以得出实验获得的目标的速度为0.9455m/s，距离为95.63m；实测距离和速度均在实际距离的误差范围内，说明杂波对消之后所得的目标点位置正确。

从图2(a)和图2(c)中可以看出本目标点附近的杂波被完全抑制，表明用本发明可以得很好的对消效果。

Claims (3)

1.一种基于混合高斯分布的杂波抑制方法，包括：

(1)采集雷达接收到的回波数据X，并以矩阵的形式进行存储；

(2)通过动目标显示MTI滤波器滤除回波数据中的零频静物杂波；

(3)对滤波器输出的数据分别进行距离维和多普勒维傅里叶变换，得到含有目标、杂波的距离和多普勒信息的数据矩阵Z；

(4)根据数据矩阵Z中的纯杂波数据的分布特性，获得一个表示杂波的混合高斯分布模型：

(4a)统计数据矩阵Z中的多组纯杂波数据的幅值得到单高斯分布的方差

并计算数据矩阵Z中第一组杂波数据的平均幅值μ₀；

(4b)将四阶混合高斯的方差初始化为

将第一阶高斯分布的均值初始化为μ₀，将其他三阶高斯分布的均值初始化为零；

(4c)每次从数据矩阵Z中读取一组杂波数据，对每组数据中的每个点的值进行判断：

若符合在上一次更新得到的混合高斯分布中的任意一阶高斯分布，则按如下公式增加此阶高斯分布的权重：

y_i,t＝(1-α)·y_i,t-1+α

否则，按如下公式减小其他阶高斯分布的权重：

y_i,t＝(1-α)·y_i,t-1

其中，y_i,t-1为读取第t组数据时的第i阶高斯分布，y_i,t为读取第t组数据更新之后的第i阶高斯分布，α为权值因子；

(4d)按如下公式更新此阶高斯分布的均值和方差：

μ_i,t＝(1-ρ)·μ_i,t-1+ρX_t

其中，μ_i,t-1为读取第t组数据时的第i阶高斯分布的均值，μ_i,t为读取第t组数据更新之后的第i阶高斯分布的均值，X_t为读取的第t组数据的每一点的幅值，ρ为学习因子，

为读取第t组数据时的第i阶高斯分布的方差，

为读取第t组数据更新之后的第i阶高斯分布的方差，T表示转置；

(4e)按如下公式更新杂波的混合高斯分布模型：

其中，P(X_t,μ_t,σ_t)为读取第t组数据时的混合高斯分布模型，X_t为读取的第t组数据的每一点的幅值，ω_i,t为第t组数据第i阶高斯分布所占权值，y(X_t,μ_i,t,σ_i,t)为第i阶高斯分布，μ_i,t为第t组数据第i阶高斯分布的均值，σ_i,t为第t组数据第i阶高斯分布的协方差，μ_t为第t组数据更新之后的混合高斯分布模型的均值，σ_t为第t组数据更新之后的混合高斯分布模型的协方差；

(4f)重复步骤(4c)～(4e)，更新得到稳定的杂波的混合高斯分布模型；

(5)根据数据矩阵Z获得只含有目标的数据矩阵G：

(5a)依次取数据矩阵Z中含目标和杂波的每组数据中每一点k(i,j)的幅值z_i,j，按如下公式进行混合高斯判断处理得到杂波抑制后输出的数据矩阵A：

其中，A_i,j为数据矩阵A中第i行第j列的点，μ为高斯分布模型的均值；

(5b)对数据矩阵A做质心凝聚处理获得的目标数据矩阵G。

2.根据权利要求1所述的方法，其特征在于：步骤(2)所述的通过动目标显示MTI滤波器滤除回波数据中的零频静物杂波，按照下式进行：

Y＝W*X

其中，W表示MTI滤波器的响应函数，X表示输入的回波数据，*表示卷积操作，Y表示MTI滤波器的输出数据。

3.根据权利要求1所述的方法，其特征在于：步骤(5b)中对数据矩阵A做质心凝聚处理，按如下步骤进行：

(5b1)依次取数据矩阵A中相邻的距离区域

其中A_i,j表示数据矩阵A中第i行第j列的的点，A_i+1,j表示数据矩阵A第i+1行第j列的点，i＝1,2,3,…M，M是数据矩阵A中总的列数，j＝1,2,3,…,Q，Q是数据矩阵A中总的行数；

(5b2)依次取数据矩阵A中相邻的多普勒区域[A_i,j…A_i,j+4]，其中，A_j+4表示数据矩阵A第i行第j+4列的点；

(5b3)根据以上取得的距离区域和多普勒区域，得到距离多普勒矩阵区域

(5b4)按如下公式对距离多普勒矩阵区域B进行处理，得到目标数据矩阵G，其中目标数据矩阵G中第i行第j列的点数据G_i,j表示如下：

Images