CN101697007B - 一种基于雷达图像的飞鸟目标识别与跟踪方法 - Google Patents

Abstract

本发明涉及一种基于雷达图像的飞鸟目标识别与跟踪方法，该方法通过从探鸟雷达PPI图像中减去背景图像，得到背景差分的雷达图像；然后对背景差分雷达图像进行杂波抑制，然后进行目标信息提取和多目标跟踪，最后通过数据融合，实现目标识别与跟踪。本发明提供的方法不仅填补了我国“雷达探鸟”领域的空白，为民航鸟击防范工作提供了技术保障，而且也适用于低空低速情况下其它雷达目标信号的处理，具有广阔的应用前景；该方法可以用于指导机场工作人员实施驱鸟作业。

Description

一种基于雷达图像的飞鸟目标识别与跟踪方法

技术领域

本发明属于图像处理技术领域，涉及一种基于雷达图像的飞鸟目标识别与跟踪方法，主要应用于机场雷达探鸟系统中雷达图像的后端处理。

背景技术

鸟击是指航空器起降或飞行过程中和鸟类、蝙蝠等飞行物相撞的事件。自从有了航空器就有了鸟击事件，据统计，仅北美地区每年鸟击造成的军用和民用航空损失就超过5亿美元并且危及乘客生命安全，欧洲每万次飞机起降平均有5.7次鸟击。中国民航在2001-2007年间，报告的鸟击事件共1055次，导致事故征候215次，造成中等程度以上损伤217起，在维修中统计的直接损失超过2.7亿人民币，间接损失更加难以计算。观测和统计鸟类迁徙和在机场附近的规律，对于鸟击防范具有重要意义。

传统的机场鸟情观测依靠人工，但在目测困难的黎明、黄昏和夜晚，恰恰是鸟击事件的高发期。据中国民航在2001～2007年间的统计，黎明、黄昏和夜晚等能见度差的条件下发生的鸟击事件占鸟击事件总数的55％，迫切需要相关技术手段的支持来降低鸟击事件的发生。雷达探鸟是鸟情观测的重要技术手段，其优点在于不受能见度和恶劣天气等因素的限制，能够全天候自动运行。机场雷达探鸟系统是鸟击防范的有效技术手段，但是目前对于基于雷达图像的飞鸟目标识别与跟踪方法还没有系统的成熟的方法和实践。

发明内容

本发明针对现有雷达探鸟技术中存在不足，为了降低由于鸟击事件而带来的飞行事故，从而减少航空损失，提出了一种基于雷达图像的飞鸟目标识别与跟踪方法。

所述的基于雷达图像的飞鸟目标识别与跟踪方法，包括背景差分、噪声抑制、目标信息提取、多目标跟踪和数据融合五个步骤，具体如下：

步骤一、背景差分，从探鸟雷达PPI图像中减去背景图像，即得到背景差分的雷达图像；背景图像的生成方法为平均值法或主成分分析法。

步骤二、对经过背景差分的雷达图像进行杂波抑制；采用恒虚警阈值分割和形态学的方法对经过背景差分的雷达图像进行降噪处理，去除其中残留的杂波信息；

步骤三、进行目标信息提取；确定目标在图像中是否为独立的区域，共有几个区域，对不同的区域进行标识。提取飞鸟目标信息包括鸟类数量、鸟类大小、坐标位置、飞行速度。

步骤四、多目标跟踪；多目标跟踪分为目标状态估计和数据关联两个主要方面，基于蒙特卡罗的思想，目标状态估计采用卡尔曼滤波，多目标数据关联采用粒子滤波数据关联方法对量测值进行目标关联。

步骤五、数据融合，实现目标识别与跟踪。将多目标跟踪获得的飞鸟目标平滑轨迹与卫星地图或坐标系相融合，生成含飞鸟目标轨迹的融合图像，便于机场工作人员观测使用。

本发明提供的飞鸟目标识别与跟踪方法不仅填补了我国“雷达探鸟”领域的空白，为民航鸟击防范工作提供了技术保障，而且也适用于低空低速情况下其它雷达目标信号的处理，具有广阔的应用前景。其优点主要体现在以下几点：

(1)本发明提供的基于雷达图像的飞鸟目标识别与跟踪方法，将飞鸟目标运动的平滑轨迹从复杂的雷达PPI图像序列中提取出来，并将其与机场地图或坐标系相融合，方便观测人员使用。

(2)本发明采用平均值法提取背景图像信息，计算量小，能够满足实时要求。

(3)本发明采用的杂波抑制方法对经过背景差分的雷达PPI图像进行处理，能够去除背景差分图像中残留的边缘杂波信息，初步锁定了运动目标，杂波抑制中设定较低的分割阈值，提高了检测率。

(4)本发明采用的目标信息提取方法对经过杂波抑制的雷达PPI图像进行处理，提取出包含部分杂波的雷达量测值。

(5)本发明采用的多目标跟踪方法通过数据关联和目标状态估计的方法，对真实的飞鸟目标进行跟踪的同时剔除杂波，降低了虚警率。

(6)本发明采用的数据融合方法将飞鸟目标平滑轨迹标定在机场地图或坐标系上，能直观地反映出飞鸟数量和所在位置，指导机场工作人员实施驱鸟作业。

附图说明

图1为本发明的飞鸟目标识别与跟踪方法流程图；

图2为CFAR检测器结构示意图；

图3为二值图像示意图；

图4为二值图像扫描示意图；

图5为二值图像标记示意图；

图6为二值图像目标标识示意图；

图7为雷达多目标跟踪算法流程图；

图8a是一幅原始的探鸟雷达图像；

图8b是背景差分后的雷达图像；

图8c是杂波抑制后的雷达图像；

图8d是含目标运动轨迹的融合图像。

具体实施方式

下面结合附图和实施例对本发明提供的基于雷达图像的飞鸟目标识别与跟踪方法进行详细说明。

本发明提供的飞鸟目标识别与跟踪方法通过如下步骤实现：

步骤一：背景差分。

所述的背景差分是指从探鸟雷达PPI图像中减去背景图像，所述的背景图像通过平均值法或主成分分析法进行构造。

平均值法是最常用、最简单的背景构造方法，这种方法通常适应于场景内的目标滞留时间较短、目标出现并不频繁的情况，针对每一帧图像重新构造背景，计算公式如下：

$B_k=\frac{1}{N}(f_k+f_{k-1}+\·\·\·+f_{k-N+1})$ (1)

$=B_{k-1}+\frac{1}{N}(f_k-f_{k-N})$

式中，N为重构背景所用的图像数，B_k为重构后的图像，b_k-1为针对上一帧构造的背景图像，f_k为第k帧图像。本发明中，每隔一定时间(5～10min)就重新构造一次背景信息。

主成分分析(PCA)是一种线性降维方法，它把高维数据投影到低维子空间，是一个从可能的高维数据集中提取特征的极端有力的工具。PCA是对描述观测数据的坐标系统的一个正交变换，它旨在用原始变量的线性组合获得较少的互不相关的新变量，同时尽可能多地保持输入数据集所包含的信息。雷达图像中的每帧图像都包含背景和运动目标(飞鸟)，因此背景可视为该图像序列最大的主成分。归结起来，可按以下步骤计算探鸟雷达PPI图像序列的主成分：

(1)计算雷达图像序列的样本方差矩阵S。设x₁，...，x_n为n个由雷达图像数据展开的一维观测样本向量，组成观测样本矩阵X，进而计算该矩阵的方差矩阵S；

(2)计算方差矩阵S的特征值及其对应的特征向量；

(3)特征值按降序排列；

(4)选择最大的特征值对应的特征向量为最大的主成分，将该特征向量还原为二维雷达PPI图像，即探鸟雷达背景图像。

从探鸟雷达PPI图像中减去背景图像，即得到背景差分的雷达图像。

步骤二：对经过背景差分的雷达图像进行杂波抑制。

对于经过背景差分的雷达图像，去除了主要的背景信息，但仍含有大量杂波，特别是边缘杂波，需要通过恒虚警(CFAR)阈值分割和形态学的方法进行杂波抑制。

CFAR阈值分割能够根据背景杂波功率的变化，自动的获取检测门限，以保持虚警不变的特性，是一种提供检测阈值的雷达信号处理方法。所述的CFAR阈值分割通过CFAR检测器实现，所述CFAR检测器的结构如图2所示，该检测器包括参考单元1、检测单元2、保护单元3，以及比较器4和乘法器5，参考单元1、检测单元2和保护单元3共占用N+M+1个单元，其中，前N/2和后N/2个单元为参考单元1，中间的1个单元为检测单元2，检测单元2两侧与参考单元1各相隔M/2个保护单元3，信号以串行的方式进入检测器各单元。CFAR检测器根据N个参考单元信号得到一个背景强度相对估计值Z，其估计方法与所采用的CFAR检测方式有关，可以采用均值类、统计排序类、削减平均类等方法进行背景强度估计。在乘法器5中，将估计值Z乘上一个门限加权系数T，得到判决门限TZ。门限加权系数T通常由下式确定：

$T=P_{f0}^{-1/N}-1$

其中P_f0表示虚警率。在比较器4中，判决门限TZ与检测单元信号进行比较判决，如果检测单元信号强度大于判决门限TZ，则判定为目标，反之为杂波。

经过CFAR阈值分割的雷达图像称为二值图像，其中目标由亮色表示，初步确定了目标区域。此二值图像由形态学作进一步处理，去除像素数较少的杂波区域。

形态学处理是指以腐蚀和膨胀这两种方法为基础，分别运用这两种方法及其组合可以进行图像形状和结构的分析和处理，包括图像分割、特征抽取、边界检测、图像滤波和恢复等。

1)图像腐蚀；

腐蚀的作用是消除目标边界点，使边界向内部收缩的过程，可以把小于结构元素的目标去除。这样选取不同大小的结构元素，就可以去除不同大小的目标。如两个物体间有细小的连通，通过腐蚀可将两个物体分开。腐蚀的数学表达式是

$S=X\⊗B=\{x,y|B_{\mathrm{xy}}\⊆X\}---\left(2\right)$

式中，S表示腐蚀后的二值图像集合，B表示用来进行腐蚀的结构元素，结构元素内的每个元素取值为0或1，它可以组成任何一种形状的图形，在B图形中有一个中心点；X表示原图像经过二值化后的像素集合。此公式的含义是用B来腐蚀X得到的集合S，S是由B完全包括在X中时B的当前位置的集合。

对于经过CFAR阈值分割的二值雷达图像，通常是拖动结构元素B在雷达图像域X中移动，横向移动间隔取1个像素，纵向移动间隔取1个扫描行。在每一个位置上，当结构元素B的中心点平移到X图像上的某一点(x，y)，如果结构元素内的每一个像素都与以(x，y)为中心的相同邻域中对应像素完全相同，那么就保留(x，y)像素点，对于原图中不满足条件的像素点则全部删除，从而达到使物体边界向内收缩的效果。腐蚀实际上是把图像的外围去掉，同时保留图像内部的部分。

2)图像膨胀；

膨胀的作用与腐蚀的作用正好相反，它是对二值化目标边界点进行扩充，将与目标区域接触的所有背景点合并到该目标中，使边界向外部扩张的过程。如果两个目标之间的距离比较近，则膨胀运算可能会把两个目标连通到一起，膨胀对填补图像分割后目标中的空洞很有用。膨胀的数学表达式是

$S=X\⊕B=\{x,y|B_{\mathrm{xy}}\∩X\≠\mathrm{\Φ}\}---\left(3\right)$

此公式的含义是用B来膨胀X得到的集合S，S是由B映像的位移与X至少有一个像素相同时B的中心点位置的集合。

对于经过CFAR阈值分割的二值雷达图像，通常是拖动结构元素B在雷达图像域X中移动，横向移动间隔取1个像素，纵向移动间隔取1个扫描行。在每一个位置上，当结构元素B的像素与目标区域至少有一个像素相交，那么就保留(x，y)像素点，从而达到使目标边界向外扩张的效果。膨胀实际上是把图像的外围进行扩充，同时保留图像内部的部分。

在实际的雷达图像处理过程中，膨胀和腐蚀经常结合使用。一幅图像往往经过一系列的膨胀与腐蚀处理，膨胀和腐蚀的次数可以任意选择，采用相同或不同的结构元素。

步骤三：进行目标信息提取。

经过背景差分和杂波抑制后获得的二值图像，需要进行目标信息提取，首先必须确定目标在图像中是否为独立的区域，共有几个区域，这就需要对不同的区域进行标识。提取飞鸟目标信息包括但不限于鸟类数量、鸟类大小、坐标位置、飞行速度。

二值图像的部分区域示例如图3所示，图中A代表目标区域，O代表背景，规定用四连通准则加标记。由于扫描有一定的次序，对任一点来说，当前点的左前一点和上一点必然是已经扫描过了的点，在扫描过程中遇到目标区域上点P，则其上点及左点必然是已经标记过了的点，对P点加标记的方法是由左点及上点来确定的，主要有下面几种不同的情况：

(a)当左前一点和上一点皆为背景O，则P点加新标记；

(b)当左前一点和上一点有一个为O，另一个为已加标记，则点P和已知标记的邻点加上相同标记；

(c)当左前一点和上一点两个邻点皆为已加标记，则P点标记与左点标记相同。

根据上面的三原则，在第一次扫描后所有目标区域上皆已加标记，如图4，标记依次为1，2，3，4，...这时图像中的同一目标区域可能有几种不同的标记，因此需要第二次扫描，来把同一目标上的标记统一起来，只要是4连通的都属于同一目标，其标记都应该一致，如图5所示标记为同一个目标区域。同理，对于8连通加标记也是如此，这时任一点P加标记需有当前点P的上点、左上点、右上点和左前点共4点来决定；若4个点皆为背景点，则该点加新标记，若其中有一邻点已加标记，则该点也加以相同的标记；若4邻点有两个以上已加标记，则可以判断与左前点加相同标记，或与右上点加相同的标记。第二遍扫描把同一目标区域的标记统一化，这样就得到每个目标的单一的、完整的标记了。

目标信息提取方法以存在着一定面积白区域的二值图像为处理对象，在此二值图像中，相互连接的白像素的集合称为一个白区域。目标信息提取时首先采用8连通判别方法进行目标区域标识，如图6中所示的A、B、C三个不连通的目标区域，对这三个目标的标识过程如下：

1)从左到右，从上到下逐个像素扫描；

2)若该点的左上、正上、右上及左前点共4个点的像素值都不为目标，则把数标加1，且此数组值为1；

3)采用(行坐标，列坐标)方式标记目标，若遇到(1，1)像素为目标A，依次判断该像素点的右上点(0，2)，正上(0，1)，左上(0，0)及左前(1，0)是否为该目标，优先级依次降低的顺序为右上点(0，2)，正上点(0，1)，左上点(0，0)及左前点(1，0)。

4)若右上点为目标，则当前点跟右上点，并标记和右上点相同的值。如：当前点(2，2)，则其右上点(1，3)为目标，所以当前点(2，2)标记与右上点(1，3)相同的值。

5)若右上点不是目标，则判断正上点。如：当前点(5，4)，则其右上点(4，5)不是目标，则判断正上点(4，4)为目标，所以当前点(5，4)标记与正上点(4，4)相同的值。

6)同理，若当前点的右上点、正上点都不为目标，则相同的方法依次判断左上点，若左上点也不为目标，则再判断左前点。

7)若如：当前点(1，8)的右上点(0，9)，正上点(0，8)，左上点(0，7)及左前点(1，7)都不为目标，则当前点的值在原来的标记上加1，以此标记作为与原目标的区别。

8)其中有一特别调整：从图6可见(10，2)是一个新加标记的点，当前点(10，3)的右上点(9，4)及左前点(10，2)为不同标记，正上点和左上点不为目标，则当前点(10，3)标记同右上点(9，4)置相同的值。此时，从头到尾扫描图像，把所有同(10，2)标记相同的像素值都标记成与右上点(9，4)同样的值。有多少个像素点转换，则统计右上点像素值(标记值)的数组就加多少，而把统计左前点的像素值(标记值)的数组置O。对图像中每个像素进行标记操作之后，将目标区域的像素值改为标号，计算各个标号的总合，即求得不同目标区域的像素数目n。PPI图像中心定义为坐标原点，x轴水平向右，y轴垂直向上。中心坐标(x₀，y₀)由(2)式得到，其中S为单一目标连通区域，量程系数C反映了选择不同量程时每个像素代表的实际距离。

$x_0=C\underset{(x,y)\∈S}{\mathrm{\Σ}}x/n,$ $y_0=C\underset{(x,y)\∈S}{\mathrm{\Σ}}y/n---\left(2\right)$

步骤四：多目标跟踪。经过目标信息提取的探鸟雷达PPI图像，初步提取了雷达量测信息，包括每个区域的中心坐标和像素数目。通过多目标跟踪，将跟踪到的雷达量测最终确定为飞鸟目标。

多目标跟踪大致可以分为目标状态估计和数据关联两个主要方面：前一方面提供跟踪需要的状态估计(预测)值，主要问题是数据精度；后一方面提供量测与目标的对应关系，即量测与航迹的关联以及杂波的剔除，主要问题是数据关联的正确性。雷达多目标跟踪方法流程如图7所示。该方法基于蒙特卡罗的思想，判断结果表现为一组离散的样本，目标状态估计采用卡尔曼滤波。每一时刻共获得T个测量值，包括n个目标测量值以及一定数量的杂波，具体步骤为：

(a)当前时刻测量值与当前时刻的目标状态预估值进行多目标数据关联；

所述的当前时刻的目标状态预估值由卡尔曼状态预估生成，目标状态预估值对实时获得的当前时刻T个测量值进行目标关联，包括目标的起始、延续和终结过程，排除杂波干扰。

卡尔曼状态预估是指卡尔曼滤波器根据上一时刻的量测更新值预估出系统下一时刻的状态，用方程表示如下：

$m_k^{-}=A_{k-1}{m}_{k-1}$ (3)

$P_k^{-}=A_{k-1}{P}_{k-1}{A}_{k-1}^T+Q_{k-1}$

式(3)中，m_k ^-和P_k ^-是k时刻获得量测之前预估的均值和方差，A_k为状态转移矩阵，Q_k为噪声矩阵。

所述的多目标数据关联采用粒子滤波数据关联方法对量测值进行目标关联。从某种程度上说，目标状态估计是整个多目标跟踪算法的关键。因为量测与航迹的关联也需要知道量测到达时刻每条航迹的状态预测值，否则量测与航迹的正确关联几乎是不可能的，精确的状态预测值，是得到正确关联的前提。基于粒子滤波(particle filtering，PF)的数据关联方法，关联结果表现为一组离散的样本，给定重要性分布π(λ_k|λ_1:k-1 ⁽ⁱ⁾，y_1:k)，一组粒子{w_k-1 ⁽ⁱ⁾，λ_k-1 ⁽ⁱ⁾，m_k-1 ⁽ⁱ⁾，P_k-1 ⁽ⁱ⁾:i＝1，...，N}，量测值y_k，一组粒子{w_k ⁽ⁱ⁾，λ_k ⁽ⁱ⁾，m_k ⁽ⁱ⁾，P_k ⁽ⁱ⁾:i＝1，...，N}处理步骤如下：

(1)基于先前生成的潜在关联变量λ_k-1 ⁽ⁱ⁾，对于每个粒子i＝1，...，N的均值m_k-1 ⁽ⁱ⁾和协方差P_k-1 ⁽ⁱ⁾进行卡尔曼预估。

(2)由相应的重要性分布给出每个粒子i＝1，...，N的新潜在关联变量λ_k ⁽ⁱ⁾

λ_k ⁽ⁱ⁾～π(λ_k|λ_1:k-1 ⁽ⁱ⁾，y_1:k)           (5)

(3)计算新的(非归一化)权重：

$w_k^{*\left(i\right)}\∝w_{k-1}^{*\left(i\right)}\frac{p\left(y_k\right|{\mathrm{\λ}}_{1:k}^{\left(i\right)},y_{1:k-1})p\left({\mathrm{\λ}}_k^{\left(i\right)}\right|{\mathrm{\λ}}_{k-1}^{\left(i\right)})}{\mathrm{\π}\left({\mathrm{\λ}}_k\right|{\mathrm{\λ}}_{1:k-1}^{\left(i\right)},y_{1:k})}---\left(6\right)$

其中似然项为卡尔曼滤波器边缘测量相似度

$p\left(y_k\right|{\mathrm{\λ}}_{1:k}^{\left(i\right)},y_{1:k-1})$ (7)

$=N\left(y_k\right|H_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right)m_k^{-\left(i\right)},H_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right)P_k^{-\left(i\right)}{H}_k^T\left({\mathrm{\λ}}_k^{\left(i\right)}\right)+R_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right))$

其中卡尔曼滤波器的模型参数基于潜在关联变量λ_k ⁽ⁱ⁾确定。

(4)权重归一化

$w_k^{\left(i\right)}=\frac{w_k^{*\left(i\right)}}{{\mathrm{\Σ}}_{j=1}^N{w}_k^{*\left(i\right)}}---\left(8\right)$

(5)基于潜在关联变量λ_k ⁽ⁱ⁾给出每个粒子的卡尔曼滤波更新。

(6)估计有效的粒子数目

$n_{\mathrm{eff}}\≈\frac{1}{{\mathrm{\Σ}}_{i=1}^N{\left(w_k^{\left(i\right)}\right)}^2}---\left(9\right)$

如果粒子的有效数目太低(比如n_eff＜N/10)，则进行重采样。确定了一组粒子之后，滤波分布近似为

$p(x_k,{\mathrm{\λ}}_k|y_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)}\mathrm{\δ}({\mathrm{\λ}}_k-{\mathrm{\λ}}_k^{\left(i\right)})N\left(x_k\right|m_k^{\left(i\right)},P_k^{\left(i\right)}).---\left(10\right)$

(b)对每个关联后的目标进行卡尔曼状态更新，分别获得每个目标的状态更新值。更新(update)部分根据当前时刻的量测值估计出系统的当前状态，用方程表示如下：

$v_k=y_k-H_k{m}_k^{-}$

$S_k=H_k{P}_k^{-}{H}_k^T+R_k$

$K_k=P_k^{-}{H}_k^T{S}_k^{-1}---\left(4\right)$

$m_k=m_k^{-}+K_k{v}_k$

$P_k=P_k^{-}-K_k{S}_k{K}_k^T$

式(4)中，y_k为k时刻获得的量测，H_k为k时刻的量测矩阵，m_k ^-和P_k ^-是k时刻获得量测之前预估的均值和方差；m_k和P_k是k时刻获得量测之后估计的均值和方差；v_k是k时刻的测量修正；S_k是k时刻的测量预估方差；K_k是滤波增益，其给出了k时刻预估值应该修正的程度。

(c)对全部卡尔曼滤波结果进行卡尔曼平滑处理，得到每个目标的平滑轨迹。

对卡尔曼滤波器获得的滤波结果进行平滑处理，其均值m_k ^s和方差P_k ^s由下式计算：

$m_{k+1}^{-}=A_k{m}_k$

$P_{k+1}^{-}=A_k{P}_k{A}_k^T+Q_k$

$C_k=P_k{A}_k^T{\left[P_{k+1}^{-}\right]}^{-1}---\left(11\right)$

$m_k^s=m_k+C_k[m_{k+1}^s-m_{k+1}^{-}]$

$P_k^s=P_k+C_k[P_{k+1}^s-P_{k+1}^{-}]C_k^T$

式中

·m_k ^s和P_k ^s是k时刻状态平均值和方差的平滑估计；

·m_k和P_k是k时刻状态平均值和方差的滤波估计；

·m_k+1 ^-和P_k+1 ^-是k+1时刻预估的状态平均值和方差，同卡尔曼滤波中的情况相同；

·C_k是k时刻的平滑增益，给出了该时刻平滑估计需要修正的程度。

步骤五：数据融合，实现目标识别与跟踪。将多目标跟踪获得的飞鸟目标平滑轨迹与卫星地图或坐标系相融合，生成含飞鸟目标轨迹的融合图像，便于机场工作人员观测使用。

下面结合某探鸟雷达PPI图像处理的全过程，对本发明提供的飞鸟目标识别与跟踪方法作详细介绍。

步骤一、背景差分。

一幅原始的探鸟雷达PPI图像如图8a所示，背景差分将背景图像从探鸟雷达图像中去除，生成背景图像的方法为平均值法或者主成分分析法。考虑到背景图像随时间存在细微变化，背景图像必须经常更新。背景差分后的雷达图像如图8b所示。

步骤二、杂波抑制对经过背景差分的雷达图像进行降噪处理，去除其中残留的杂波信息。杂波抑制包括恒虚警阈值分割和形态学。其中，恒虚警阈值分割能够根据不同的雷达图像自适应地选择阈值。经过阈值分割的图像通过形态学模块进一步降噪，去除那些像素过少的目标。杂波抑制后的雷达图像如图8c所示。

步骤三、目标信息提取从经过杂波抑制的雷达图像中提取雷达量测信息，提取的信息包括目标数量、大小和坐标位置。提取的雷达量测信息见表1。

表1雷达量测信息提取

步骤四、目标跟踪基于目标信息提取的雷达量测信息对飞鸟目标进行跟踪，输出飞鸟目标平滑轨迹，同时进行记录。通过对雷达图像序列进行多目标跟踪处理的结果可知，表1中的量测1、2、3形成了平滑轨迹，量测4、5只出现在一幅雷达图像中，未能形成轨迹，因此量测值1、2、3为目标，4、5为杂波。

步骤五、数据融合将多目标跟踪方法获得的飞鸟目标轨迹信息与卫星地图或坐标系相融合，生成含飞鸟目标平滑轨迹的融合图像，便于机场工作人员观测，如图8d所示。

Claims (7)

1.一种基于雷达图像的飞鸟目标识别与跟踪方法，其特征在于如下步骤：

步骤一、背景差分，从探鸟雷达PPI图像中减去背景图像，即得到背景差分的雷达图像；

步骤二、对经过背景差分的雷达图像进行杂波抑制处理，采用恒虚警阈值分割和形态学的方法对经过背景差分的雷达图像进行降噪处理，去除其中残留的杂波信息；

步骤三、进行目标信息提取；目标信息提取过程为：

对经过背景差分和杂波抑制后获得的二值图像进行不同区域的标识，采用从上到下，从左到右的扫描次序，目标区域中任意一点P的标记过程主要有下面几种不同的情况：

(a)当左前一点和上一点皆为背景O，则P点加新标记；

(b)当左前一点和上一点有一个为O，另一个为已加标记，则点P和已知标记的邻点加上相同标记；

(c)当左前一点和上一点两个邻点皆为已加标记，则P点标记与左点标记相同；

如果同一目标区域有几种不同的标记，则进行第二次扫描，将同一目标区域的标记统一化，得到每个目标的单一的、完整的标记；

步骤四、多目标跟踪；多目标跟踪步骤为：

(a)当前时刻测量值与当前时刻的目标状态预估值进行多目标数据关联；所述的多目标数据关联采用粒子滤波数据关联方法对量测值进行目标关联，关联结果表现为一组离散的样本，给定重要性分布

一组粒子量测值y_k，一组粒子 $\{w_k^{\left(i\right)},{\mathrm{\λ}}_k^{\left(i\right)},m_k^{\left(i\right)},P_k^{\left(i\right)}:i=1,...,N\}$ 处理步骤如下：

(1)基于先前生成的潜在关联变量对于每个粒子i＝1，...，N的均值

和协方差

进行卡尔曼预估；

(2)由相应的重要性分布给出每个粒子i＝1，...，N的新潜在关联变量

${\mathrm{\λ}}_k^{\left(i\right)}~\mathrm{\π}\left({\mathrm{\λ}}_k\right|{\mathrm{\λ}}_{1:k-1}^{\left(i\right)},y_{1:k})---\left(5\right)$

(3)得到新的权重：

$w_k^{*\left(i\right)}\∝w_{k-1}^{*\left(i\right)}\frac{p\left(y_k\right|{\mathrm{\λ}}_{1:k}^{\left(i\right)},y_{1:k-1})p\left({\mathrm{\λ}}_k^{\left(i\right)}\right|{\mathrm{\λ}}_{k-1}^{\left(i\right)})}{\mathrm{\π}\left({\mathrm{\λ}}_k\right|{\mathrm{\λ}}_{1:k-1}^{\left(i\right)},y_{1:k})}---\left(6\right)$

其中似然项为卡尔曼滤波器边缘测量相似度

$p\left(y_k\right|{\mathrm{\λ}}_{1:k}^{\left(i\right)},y_{1:k-1})$ (7)

$=N\left(y_k\right|H_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right)m_k^{-\left(i\right)},H_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right)P_k^{-\left(i\right)}{H}_k^T\left({\mathrm{\λ}}_k^{\left(i\right)}\right)+R_k\left({\mathrm{\λ}}_k^{\left(i\right)}\right))$

其中卡尔曼滤波器的模型参数基于潜在关联变量

确定；

(4)权重归一化

$w_k^{\left(i\right)}=\frac{w_k^{*\left(i\right)}}{{\mathrm{\Σ}}_{j=1}^N{w}_k^{*\left(i\right)}}---\left(8\right)$

(5)基于潜在关联变量给出每个粒子的卡尔曼滤波更新；

(6)估计有效的粒子数目

$n_{\mathrm{eff}}\≈\frac{1}{{\mathrm{\Σ}}_{i=1}^N{\left(w_k^{\left(i\right)}\right)}^2}---\left(9\right)$

如果粒子的有效数目太低，则进行重采样；确定了一组粒子之后，滤波分布近似为

$p(x_k,{\mathrm{\λ}}_k|y_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_k^{\left(i\right)}\mathrm{\δ}({\mathrm{\λ}}_k-{\mathrm{\λ}}_k^{\left(i\right)})N\left(x_k\right|m_k^{\left(i\right)},P_k^{\left(i\right)})---\left(10\right)$

(b)对每个关联后的目标进行卡尔曼状态更新，分别获得每个目标的状态更新值；

根据当前时刻的量测值估计出系统的当前状态，用方程表示如下：

$v_k=y_k-H_k{m}_k^{-}$

$S_k=H_k{P}_k^{-}{H}_k^T+R_k$

$K_k=P_k^{-}{H}_k^T{S}_k^{-1}---\left(4\right)$

$m_k=m_k^{-}+K_k{v}_k$

$P_k=P_k^{-}-K_k{S}_k{K}_k^T$

式(4)中，y_k为k时刻获得的量测，H_k为k时刻的量测矩阵，和

是k时刻获得量测之前预估的均值和方差；m_k和P_k是k时刻获得量测之后估计的均值和方差；v_k是k时刻的测量修正；S_k是k时刻的测量预估方差；K_k是滤波增益，其给出了k时刻预估值应该修正的程度；

(c)对全部卡尔曼滤波结果进行卡尔曼平滑处理，得到每个目标的平滑轨迹；

对卡尔曼滤波器获得的滤波结果进行平滑处理，其均值和方差

由下式得到：

$m_{k+1}^{-}=A_k{m}_k$

$P_{k+1}^{-}=A_k{P}_k{A}_k^T+Q_k$

$C_k=P_k{A}_k^T{\left[P_{k+1}^{-}\right]}^{-1}---\left(11\right)$

$m_k^s=m_k+C_k[m_{k+1}^s-m_{k+1}^{-}]$

$P_k^s=P_k+C_k[P_{k+1}^s-P_{k+1}^{-}]C_k^T$

式中和

是k时刻状态平均值和方差的平滑估计；m_k和P_k是k时刻状态平均值和方差的滤波估计；

和

是k+1时刻预估的状态平均值和方差，同卡尔曼滤波中的情况相同；C_k是k时刻的平滑增益，给出了该时刻平滑估计需要修正的程度；

步骤五、数据融合，实现目标识别与跟踪。

2.根据权利要求1所述的飞鸟目标识别与跟踪方法，其特征在于：所述的步骤一中背景图像的生成方法为平均值法或主成分分析法。

3.根据权利要求2所述的飞鸟目标识别与跟踪方法，其特征在于：所述的平均值法生成的背景图像为：

$B_k=\frac{1}{N}(f_k+f_{k-1}+\·\·\·+f_{k-N+1})$ (1)

$=B_{k-1}+\frac{1}{N}(f_k-f_{k-N})$

式中，N为重构背景所用的图像数，B_k为重构后的图像，B_k-1为针对上一帧构造的背景图像，f_k为第k帧图像。

4.根据权利要求2所述的飞鸟目标识别与跟踪方法，其特征在于：所述的主成分分析法具体步骤为：

(1)确定雷达图像序列的样本方差矩阵S：设x₁，...，x_n为n个由雷达图像数据展开的一维观测样本向量，组成观测样本矩阵X，进而计算该矩阵的方差矩阵S；

(2)确定方差矩阵S的特征值及其对应的特征向量；

(3)特征值按降序排列；

(4)选择最大的特征值对应的特征向量为最大的主成分，将该特征向量还原为二维雷达PPI图像，即探鸟雷达背景图像。

5.根据权利要求1所述的飞鸟目标识别与跟踪方法，其特征在于：所述步骤二中的恒虚警阈值分割通过恒虚警检测器实现，所述恒虚警检测器包括参考单元、检测单元、保护单元，以及比较器和乘法器，参考单元、检测单元和保护单元共占用N+M+1个单元，其中，前N/2和后N/2个单元为参考单元，中间的1个单元为检测单元，检测单元两侧与参考单元各相隔M/2个保护单元，单元信号以串行的方式进入检测器各单元；恒虚警检测器根据N个参考单元信号得到一个背景强度相对估计值Z，在乘法器中，将估计值Z乘上一个门限加权系数T，得到判决门限TZ，判决门限TZ与检测单元信号在比较器中进行比较判决，如果检测单元信号强度大于判决门限TZ，则判定为目标，反之为杂波。

6.根据权利要求5所述的飞鸟目标识别与跟踪方法，其特征在于：门限加权系数T由下式确定：

$T=P_{f0}^{-1/N}-1$

其中P_f0表示虚警率。

7.根据权利要求1所述的飞鸟目标识别与跟踪方法，其特征在于将由多目标跟踪算法提取的飞鸟目标平滑轨迹与卫星地图或坐标系相融合，生成便于观测的融合图像。

Images