基于杂波轮廓图的动目标显示滤波器滤波方法与流程

本发明属于信号处理技术领域，特别涉及一种基于杂波轮廓图的动目标显示滤波器滤波方法，是一种海杂波抑制技术，适用于目标检测与跟踪，便于工程实现。

背景技术：

杂波抑制技术是雷达在杂波背景下检测目标的一种重要手段，目前用于杂波抑制的措施有很多，其中，从频率域利用杂波和目标回波的多普勒信息的不同来抑制杂波，是雷达信号处理中最常用的方法，称之为动目标显示mti或动目标检测mtd技术。为解决盲速问题，常用的方法是利用参差重复频率，此时只能选择mti技术；而工程实现中，杂波抑制常采用传统自适应mti或者级联mti技术。

传统自适应mti技术采用运动杂波谱中心补偿抑制法，先通过杂波谱中心估计，在对谱中心补偿，然后进行杂波抑制；工程实现中，测频不仅会占用fpga的大量逻辑资源和存储资源，在频率估计不准确时还会导致设计的滤波器凹口对不准杂波中心频率而无法完全对消杂波。

级联mti技术多采用两级级联mti滤波器，第一级mti滤波器对消地杂波，第二级mti滤波器对消剩余的动杂波；级联mti技术不仅需要测频，还要对第一级的滤波结果进行相位校正，不仅占用大量的资源，而且计算量大。

无论是传统的mti技术还是优化设计的mti滤波器，在工程实现中都需要对雷达回波数据作测频处理；基于fpga的杂波抑制的工程实现中，测频不仅会占用fpga的大量逻辑资源和存储资源，在频率估计不准确时还会导致设计的滤波器凹口对不准杂波中心频率而无法完全对消杂波。

技术实现要素：

针对上述现有技术存在的不足，本发明的目的在于提出一种基于杂波轮廓图的动目标显示滤波器滤波方法，该种基于杂波轮廓图的动目标显示滤波器滤波方法不通过测频根据建立的杂波轮廓图选择滤波器组的输出结果，以减少工程实现中的运算量，节约资源。

本发明的技术思路:将雷达作用域(处理范围)划分成多个杂波单元，采用“区域扩展”技术建立杂波轮廓图；根据雷达的工作频段和最大海情计算海杂波最大多普勒频移，在频域上按最大多普勒频移分成若干多普勒通道，每个通道有不同的中心频率；计算一组mti滤波器系数，mti滤波器的个数与凹口中心与划分的多普勒通道一致；根据杂波轮廓图选择滤波器组的输出结果。

为达到上述技术目的，本发明采用如下技术方案予以实现。

一种基于杂波轮廓图的动目标显示滤波器滤波方法，包括以下步骤：

步骤1，确定雷达，并对雷达的检测范围进行划分，得到m”×n”个目标分辨单元和m'×n'个杂波单元，将m”×n”个目标分辨单元记为m”×n”维目标分辨单元矩阵，将m'×n'个杂波单元记为m'×n'维杂波单元矩阵；m”>m'，n”>n'，m”、n”、m'、n'分别为大于0的正整数；

初始化：令na表示第na个天线扫描周期，na∈{0,1,…,y}，na的初始值为0；y表示设定的天线扫描周期个数；

设定第na个天线扫描周期时的完整杂波图，且第na个天线扫描周期时的完整杂波图距离向上杂波单元个数为k1，第na个天线扫描周期时的完整杂波图方位向上杂波单元个数为k2，其中，表示向下取整；

令i∈{1+(1-1)m',1+(2-1)m',1+(3-1)m',……,1+(k1-1)m'}，j∈{1+(1-1)n',1+(2-1)n',1+(3-1)n',……,1+(k2-1)n'}，m'表示m'×n'维杂波单元矩阵距离向最大编号，n'表示m'×n'维杂波单元矩阵方位向最大编号；

步骤2，计算得到第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为i、方位向编号为j处的杂波单元平均值

步骤3，分别令i分别取1+(1-1)×m'至1+(k2-1)n'，重复执行步骤2，进而得到第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为1+(1-1)×m'、方位向编号为j处的杂波单元平均值至第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为1+(k1-1)m'、方位向编号为j处的杂波单元平均值记为第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为j处的k1个杂波单元平均值

步骤4，分别令j取1+(1-1)n'至1+(k2-1)n'，将i的值置为1，重复执行步骤2和步骤3，进而分别得到第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为1+(1-1)n'处的k1个杂波单元平均值至第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为1+(k2-1)n'处的k1个杂波单元平均值并记为第na个天线扫描周期时的完整杂波图

步骤5，根据第na个天线扫描周期时的完整杂波图，得到第na个天线扫描周期时的杂波轮廓图；其中q为大于0的正整数；

步骤6，令na的值加1，返回步骤2，直到得到第y个天线扫描周期时的杂波轮廓图，并将此时得到的第y个天线扫描周期时的杂波轮廓图，作为记为最终杂波轮廓图；

步骤7，计算海杂波最大多普勒频移；

步骤8，根据海杂波最大多普勒频移，确定多普勒区间，进而得到d个多普勒通道以及每个多普勒通道的中心频率；d为大于0的正整数；

步骤9，根据每个多普勒通道的中心频率，计算得到p个k阶mti滤波器的权系数；其中，p表示mti滤波器个数，k表示每个mti滤波器的最大阶数，p、k分别为大于0的正整数；

步骤10，确定雷达输入脉冲，并根据p个k阶mti滤波器的权系数和最终杂波轮廓图，得到最终的k阶mti滤波器输出结果，所述最终的k阶mti滤波器输出结果为基于杂波轮廓图的动目标显示滤波器滤波结果。

本发明与现有mti方案相比具有以下优点：

1)避免频率估计不准的问题

传统mti方案根据估计的海杂波多普勒中心频率设置滤波器凹口中心，但实际工程中由于中心频率估计不准造成杂波抑制效果不好；本发明无需进行频率估计，避免了频率估计不准的问题，而且在动杂波频率变化较大的情况下，有很好的杂波抑制效果。

2)计算量少，节约资源

级联mti滤波器虽然可以滤除地杂波和海杂波，但需要相位校正；本发明只进行一级mti滤波，不需要进行相位校正，计算量少，节约逻辑资源和存储资源。

3)易于工程实现

本发明方法简单，容易用fpga完成工程实现。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1为本发明的一种基于杂波轮廓图的动目标显示滤波器滤波方法流程图；

图2(a)为二维杂波单元的方位-距离划分示意图；

图2(b)为图2(a)中阴影部分的放大图；

图3为一个驻留时间内脉压处理后的幅度图；

图4为传统自适应mti滤波器的滤波结果幅度图；

图5为本发明方法的滤波结果幅度图；

图6为非相干积累结果的幅度图。

具体实施方式

参照图1，为本发明的一种基于杂波轮廓图的动目标显示滤波器滤波方法流程图；所述基于杂波轮廓图的动目标显示滤波器滤波方法，包括以下步骤：

步骤1，确定雷达，本实施例中采用脉冲多普勒雷达；然后对雷达的检测范围进行划分，得到m”×n”个目标分辨单元和m'×n'个杂波单元，每个目标分辨单元分别由目标距离尺寸δr和目标方位尺寸δθ来表示，δr＝cτ/2，c表示光速，τ表示雷达发射信号的脉冲宽度，其中一个脉冲重复周期tr对应雷达方位上的天线转动角度为方位尺寸δθ；每个杂波单元分别用杂波距离尺寸δrc和杂波方位尺寸δa表示，δrc＝m'·δr，δa＝n'δθ；m”>m'，n”>n'，m”、n”、m'、n'分别为大于0的正整数。

为了减少存储容量，本发明中设定杂波单元总个数大于目标分辨单元总个数，且本实施例中m'∈{32,33,……,256}，n'∈{16,17,……,128}；将m”×n”个目标分辨单元记为m”×n”维目标分辨单元矩阵，将m'×n'个杂波单元记为m'×n'维杂波单元矩阵，且每一个杂波单元的平均值为m”×n”个目标分辨单元矩阵的二维平均，则m'×n'维杂波单元矩阵中距离向编号为i、方位向编号为j处的杂波单元平均值为

m'表示m'×n'维杂波单元矩阵距离向最大编号，n'表示m'×n'维杂波单元矩阵方位向最大编号，n∈{1,2,……,n”}，m∈{1,2……,m”}，xi+m,j+n表示m'×n'维杂波单元矩阵中距离向编号为i+m、方位向编号为j+n处的杂波数据，i+m∈{1,2,……,m'}，j+n∈{1,2……,n'}，。

确定噪声单元电平估计值：雷达最大探测距离为lkm，将距离雷达中心9l/10km处视作为远区；选取远区m×n个目标分辨单元，并记为m×n维目标分辨单元矩阵，m×n维目标分辨单元矩阵中每一个目标分辨单元的噪声电平值为m×n维目标分辨单元矩阵的二维平均，将m×n维目标分辨单元矩阵中距离向编号为i'、方位向编号为j'处的噪声电平值记为即：

其中，n'∈{1,2,……,n}，m'∈{1,2……,m}，i'∈{1,2,……,m}，j'∈{1,2……,n}，本实施例中m'∈{64，65……512},n'∈{32，33……256}；m>m'，n>n'；xi'+m',j'+n'表示m×n维目标分辨单元矩阵中距离向编号为i'+m'、方位向编号为j'+n'处的目标数据。

初始化：令na表示第na个天线扫描周期，na∈{0,1,…,y}，na的初始值为0；y表示由雷达开机时的不稳定状态达到稳定状态过程中对应的天线扫描周期个数，本实施例中y取值为8。

设定第na个天线扫描周期时的完整杂波图，且第na个天线扫描周期时的完整杂波图距离向上杂波单元个数为k1，第na个天线扫描周期时的完整杂波图方位向上杂波单元个数为k2，其中，表示向下取整。

令i∈{1+(1-1)m',1+(2-1)m',1+(3-1)m',……,1+(k1-1)m'}，j∈{1+(1-1)n',1+(2-1)n',1+(3-1)n',……,1+(k2-1)n'}，m'表示m'×n'维杂波单元矩阵距离向最大编号，n'表示m'×n'维杂波单元矩阵方位向最大编号。

步骤2，杂波单元值更新：计算得到第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为i、方位向编号为j处的杂波单元平均值其表达式为：

其中，表示第na-1个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为i、方位向编号为j处的杂波单元平均值，xi+m,j+n(na)表示第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为i+m、方位向编号为j+n处的杂波数据，k表示设定的常数，k∈[0,1]，本实施例中k取值为7/8；n∈{1,2,……,n”}，m∈{1,2……,m”}，m”表示m”×n”维目标分辨单元矩阵距离向最大编号，n”表示m”×n”维目标分辨单元矩阵方位向最大编号。

步骤3，分别令i分别取1+(1-1)×m'至1+(k2-1)n'，重复执行步骤2，进而得到第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为1+(1-1)×m'、方位向编号为j处的杂波单元平均值至第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为1+(k1-1)m'、方位向编号为j处的杂波单元平均值记为第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为j处的k1个杂波单元平均值

步骤4，分别令j取1+(1-1)n'至1+(k2-1)n'，将i的值置为1，重复执行步骤2和步骤3，进而分别得到第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为1+(1-1)n'处的k1个杂波单元平均值至第na个天线扫描周期时m'×n'维杂波单元矩阵中方位向编号为1+(k2-1)n'处的k1个杂波单元平均值并记为第na个天线扫描周期时的完整杂波图所述第na个天线扫描周期时的完整杂波图距离向上杂波单元个数为k1，第na个天线扫描周期时的完整杂波图方位向上杂波单元个数为k2。

步骤5，确定第na个天线扫描周期时的杂波轮廓图。

5.1初始化：令k1表示第na个天线扫描周期时的完整杂波图距离向上杂波单元编号，k2表示第na个天线扫描周期时的完整杂波图方位向上杂波单元编号，其中，表示向下取整；的初始值分别为1。

5.2比较第na个天线扫描周期时的完整杂波图中距离向编号为方位向编号为处的杂波单元平均值与对应m×n维目标分辨单元矩阵中距离向编号为方位向编号为处的噪声电平值的大小，若则第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为方位向编号为处的杂波单元内有杂波，并标记为1；否则无杂波，并标记为0；其中q为大于0的正整数，本实施例中q取值为3。

5.3令不变，且令分别取1至k2，重复执行5.2，进而分别将第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为方位向编号为1处的杂波单元至第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为方位向编号为k2处的杂波单元标记为1或0，记为第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为处的杂波轮廓图。

5.4令分别取1至k1，将的值置为1，重复执行5.2和5.3，进而分别得到第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为1处的杂波轮廓图至第na个天线扫描周期时m'×n'维杂波单元矩阵中距离向编号为k1处的杂波轮廓图，记为第na个天线扫描周期时的杂波轮廓图。

步骤6，令na的值加1，返回步骤2，直到得到第y个天线扫描周期时的杂波轮廓图，所述第y个天线扫描周期时的杂波轮廓图稳定，且第y个天线扫描周期时的杂波轮廓图中包含的所有元素取值全部为0或者全部为1；并将此时得到的第y个天线扫描周期时的杂波轮廓图，作为最终杂波轮廓图。

步骤7，计算海杂波最大多普勒频移|fdmax|，|fdmax|＝2×|vmax|/λ，|vmax|表示海杂波最大风速，λ为雷达发射电磁波波长。

具体的，海表面状态有轻风、微风、和风、大风、强风、巨风、狂风，将海表面状态为强风时视作高海情，此时风速最高达到10m/s，本实施例中海杂波最大风速取值为10m/s。

步骤8，划分多普勒频率通道：根据海杂波最大多普勒频移，确定多普勒区间为(-|fdmax|,|fdmax|)，并对所述多普勒区间(-|fdmax|,|fdmax|)行划分，得到d个多普勒通道，多普勒通道间隔为2×|fdmax|÷d，第1个多普勒通道的中心频率为f1为第d个多普勒通道的中心频率为fd为第d'个多普勒通道的中心频率为其中2≤d'≤d-1，d为大于0的正整数。

步骤9，设计mti滤波器，即计算mti滤波器系数w，mti滤波器个数p与多普勒通道个数d取值相等。

mti滤波器是带阻滤波器，也是fir滤波器；本发明中将每个滤波器设置为双凹口滤波器，其中一个凹口设置在零频处，另一个凹口位置与mti滤波器编号有关，p个mti滤波器凹口中心与d个多普勒通道的中心频率一一对应，p表示k阶mti滤波器个数。

采用零点分配法设计mti滤波器，零点分配法是在凹口处设置频率响应零点的一种方法，w＝[w1,w2,…,wx,…,wp]，对于k阶mti滤波器，第x个mti滤波器的权系数为wx，wx＝(wx,1,wx,2,……,wx,b,……,wx,k)，1≤x≤p，wx,b表示第x个mti滤波器的第b阶系数，b∈[1,2,……,k]，k表示每个mti滤波器的最大阶数；工程中k的取值范围为k∈[2,3……6]，x的初始值为1。

9.1建立第x个mti滤波器的权系数wx的齐次线性方程为：

其中，为第b阶系数的脉冲重复周期的n次幂，wx,b表示第x个k阶mti滤波器的第b+1阶系数，b∈[1,2,……,k]，k表示每个k阶mti滤波器的阶数；fx表示第x个k阶mti滤波器对应的多普勒通道中心频率，1≤x≤p，p表示k阶mti滤波器的总个数，tb表示第x个k阶mti滤波器的第b+1阶系数对应的脉冲重复周期；w0为设定常数，通常设定为1。

9.2根据第x个mti滤波器的权系数wx的齐次线性方程，计算得到第x个k阶mti滤波器的权系数wx。

9.3令x分别取1至p，重复执行9.1和9.2，进而分别计算得到第1个k阶mti滤波器的权系数w1至第p个k阶mti滤波器的权系数wp，记为p个k阶mti滤波器的权系数w，

w＝[w1,w2,…,wx,…,wp]。

步骤10，将雷达输入脉冲记为x(t)，并使用p个k阶mti滤波器的权系数w对雷达输入脉冲x(t)进行滤波，得到p个k阶mti滤波器的滤波输出结果y(t)，其中第x个k阶mti滤波器的滤波输出结果为yx(t)，其表达式为：

yx(t)＝wx,1x(t)+wx,2x(t-t1)+…wx,bx(t-t1-…-tb-1)+…wx,kx(t-t1-…-tb-1-…-tk-1)

其中，t表示时间变量，wx,b表示第x个k阶mti滤波器的第b+1阶系数，b∈[1,2,……,k]，k表示每个k阶mti滤波器的阶数，x(t)表示雷达输入脉冲，x(t-t1)表示雷达在时间上延迟t1后得到的输入脉冲，x(t-t1-…-tb-1)表示雷达在时间上延迟t1至tb-1后得到的输入脉冲，x(t-t1-…-tb-1-…-tk-1)表示雷达在时间上延迟t1至tk-1后得到的输入脉冲，t1表示第x个k阶mti滤波器的第2阶系数对应的脉冲重复周期，tb-1表示第x个k阶mti滤波器的第b阶系数对应的脉冲重复周期，tk-1表示第x个k阶mti滤波器的第k阶系数对应的脉冲重复周期，t0＝0；1≤x≤p，p表示k阶mti滤波器的总个数。

当最终杂波轮廓图中所有元素取值全部为1时，最终的k阶mti滤波器输出结果为

当最终杂波轮廓图中所有元素取值全部为0时，最终的k阶mti滤波器输出结果为

其中，min(·)为求最小值函数，max(·)为求最大值函数。

所述最终的k阶mti滤波器输出结果为基于杂波轮廓图的动目标显示滤波器滤波结果。

通过以下计算和实测数据处理结果对本发明效果作进一步验证说明。

(一)运算量计算

运算量计算强杂波环境下每部脉冲多普勒雷达都有杂波图和噪声电平估值，杂波图数据与噪声电平值的每次比较中只有一次乘法，且为实数相乘，在fpga中可通过移位和加减运算实现。计算量可忽略不计。而测频运算中有复数相乘，无法通过简单的移位和加减元算实现。天线扫描一圈，有m”×n”个目标分辨单元。

假设采用比较简单的dft测频，dft点数为8，则复数乘为：假定每个脉冲有50000个距离单元，即n”＝50000，每圈有4000个脉冲，即m”＝4000，则复乘次数为：运算量十分大，所以本发明大大减少了运算量。

(二)实测数据处理和分析

实测数据处理1雷达架设在某山山腰，周围有树木、海、民房等反射物。选取某一次采集数据中的某个波位数据，图2(a)为二维杂波单元的方位-距离划分示意图，每个杂波单元分别用杂波距离尺寸δrc和杂波方位尺寸δa表示，半径rm表示杂波区的距离范围，图2(b)为图2(a)中阴影部分的放大图，一个脉冲重复周期tr对应在方位上的天线转动角度为δθ，且δrc＝m·δr和δa＝nδθ；图2(b)为图2(a)中阴影部分的放大图，图3为一个驻留时间内脉压处理后的幅度图，图4为传统自适应mti滤波器的滤波结果幅度图，图5为本发明方法的滤波结果幅度图，图6为非相干积累结果的幅度图，虚线为对自适应mti滤波器的滤波方法进行非相干积累处理，实线为对本发明方法的滤波结果进行非相干处理。

(三)实测数据处理分析

对于实测数据处理1，由图4可以看出本发明中使用mti滤波器组方法可以对消掉杂波并凸显出目标，由图5可以看出自适应mti滤波器并没有将杂波完全抑制且目标被抑制，由图6可以看出非相干积累后本发明对杂波的抑制效果更好。

因此本发明方法相比于传统的mti滤波，不仅便于工程实现，效果也更突出。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。