CN102945374A - 一种高分辨率遥感图像中民航飞机自动检测方法 - Google Patents

Abstract

本发明是基于圆周频率滤波与梯度方向直方图相结合的一种高分辨率遥感图像中民航飞机自动检测方法，该方法共有五大步骤。步骤一：读取高分辨率遥感图像数据，并进行多尺度圆周频率滤波；步骤二：提取训练正、负样本的梯度方向直方图特征；步骤三：训练基于梯度方向直方图特征的Adaboost分类器；步骤四：利用Adaboost分类器进行精检测；步骤五：通过区域生长聚类，定位飞机目标。本发明克服了高分辨率遥感图像背景复杂、伪目标繁多的问题，可以在虚警率较低的前提下保证较高的检测率，在高分辨率遥感图像目标检测技术领域里具有广阔的应用前景。

Description

一种高分辨率遥感图像中民航飞机自动检测方法

技术领域：

本发明涉及基于圆周频率滤波与梯度方向直方图（Histograms of Oriented Gradients,HOG）特征的一种高分辨率遥感图像中民航飞机自动检测方法，属于高分辨率遥感图像目标检测技术领域。

背景技术：

一直以来，对图像中感兴趣目标的检测与识别都是一项的重要研究课题。近年来，随着遥感技术的不断发展和图像分辨率的不断提高，从高分辨率遥感图像中自动地检测与识别各种目标不仅具有重要的军事应用价值，也越来越成为民航遥感领域关注的热点。其中高分辨率遥感图像中的民航飞机自动检测，作为目标检测的重要组成部分，由于其在城市规划、机场管理等方面的重要指导作用，也一直在研究和探索中。但通常情况下，由于高分辨率遥感图像的复杂性以及飞机机身的大小、对比度等因素的影响，想要在虚警较少的前提下保证较高的检测率并不容易。

目前存在的飞机检测方法有很多，大多数方法都是根据飞机机身的某些单一的特征进行检测的。比如利用飞机的角点或直线特征，还有仿射尺度不变特征或仿射不变矩描述子特征。除此之外，也有利用模板匹配技术，完成飞机目标的检测。但是这些特征所携带的信息比较单一，仅用其中一、两种往往会产生大量的虚警，难以应对高分辨率遥感图像复杂的背景；而把多种特征结合起来又毫无疑问会大大增加算法的复杂性和实现难度。本发明针对上述情况，提出一种基于圆周频率滤波与梯度方向直方图特征相结合的方法，首先通过多尺度圆周频率滤波，进行粗检测，得到民航飞机目标的候选点集；再利用训练得到的Adaboost分类器，进行精检测，取得了良好的检测效果。

发明内容：

1、目的：本发明的目的是提供基于圆周频率滤波与梯度方向直方图特征相结合的一种高分辨率遥感图像中民航飞机自动检测方法。该方法利用圆周频率滤波与基于梯度方向直方图特征的Adaboost分类器，来进行民航飞机目标检测。

2、技术方案：本发明是通过以下技术方案实现的：

本发明一种高分辨率遥感图像中民航飞机自动检测方法，它包括如下步骤：

步骤一：读入图像，并对输入图像进行圆周频率滤波，本发明读入的遥感图像分辨率约为1米。

所谓圆周频率滤波就是计算图像中每一个像素点的圆周频率响应，并把响应值大于某一阈值的像素点作为飞机目标的候选点集保留下来；小于该阈值的点则认为是非飞机区域而排除掉。其计算过程如下：

假设图像上有一坐标为(i,j)的点，f_k(k=0,1,...,N-1)表示以(i,j)为中心以r为半径的圆周上逆时针方向像素的灰度值所组成的一维数列，其中r为这个圆周的半径，N为这个圆周序列的点数。求这个一维数列的傅立叶变换的幅值f(i,j)，如下所示：

$f(i,j)={\left(\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_k\cos \frac{\mathrm{c\πk}}{N}\right)}^2+{\left(\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_k\sin \frac{\mathrm{c\πk}}{N}\right)}^2$

上式中f(i,j)为圆周上的数列f_k经过离散傅立叶变换后的幅值。c是常量，代表求离散傅立叶变换时正弦和余弦函数的周期数。由于飞机主要由机头、机翼（一对）和尾翼组成，当半径r的取值使得圆周扫过这4部分时，它们分别对应于f_k的4个峰（当飞机比背景亮）或4个谷（当飞机比背景暗）。因此飞机区域的圆周序列f_k通常有4个峰和4个谷，即近似4周期的正弦信号。而当傅里叶变换中正弦或余弦函数的周期与f_k的周期相同时，该傅里叶变换的幅值f(i,j)的绝对值最大，即在飞机中心的圆周频率响应的幅值大，因此令常数c=8，代表傅里叶变换的正弦与余弦函数有四个周期。

本发明选择圆周半径r为4像素，采样点数N为40个，常数c为8。这样，长度约为30到50像素的飞机区域会得到相对较强的响应，平坦的背景区域的响应相对较弱。此外为了检测更大的飞机，本发明还对输入图像进行两次降采样，采样率为1.30、1.80，使得大飞机也有较强的滤波响应，这在后文将称为多尺度圆周频率滤波。还有一些复杂的背景如房屋、道路边缘以及与飞机形状近似的建筑也会有较强的响应，这些点将在后面分类器精检测的过程中被排除掉。

步骤二：提取训练样本的梯度方向直方图特征。因为本发明欲利用分类器对圆周频率滤波结果进行精检测，因此这一步是为后面训练分类器服务。训练样本为在Google Earth软件上截取的大小为40*40的图片（关于训练的具体细节会在步骤三详述），用来提取特征。本发明采用的是梯度方向直方图特征，特征提取的流程如下：

首先将样本图像计算图像梯度，之后再将梯度图像分成小的连通区域，此连通区域被称为胞元(cell)。然后统计出胞元(cell)中所有像素点的梯度方向直方图。最后把这些直方图组合起来便可以构成特征描述符。为了使特征描述符对光照变化和阴影干扰具有更强的不变性，算法将这些局部直方图在图像的更大的范围内(被称为块区或block)进行对比度归一化(contrast-normalized)；归一化所采用的方法是：先计算各直方图在其所在块区(block)中的能量密度，然后根据这个能量密度对区间中的各个胞元做归一化。

步骤三：训练Adaboost分类器。步骤二结束后，用提取的梯度方向直方图特征训练Adaboost分类器。由于梯度方向直方图特征不具有方向不变性，因此需要训练多个分类器以检测不同方向的飞机。本发明训练了0°到315°等间隔的8个方向的分类器。训练样本为在Google Earth软件上截取的大小为40*40的图片。其中负样本，即非飞机图片5000幅；8个方向的正样本，即飞机图片，每个方向各有3000幅。

训练Adaboost分类器的流程如下：假定训练集为(x₁,y₁),...,(x_m,y_m)，训练样本的权重分布为D_t(i)(1≤i≤m)；其中x_i为样本特征，y_i为样本标签。假定起初时，所有的权重D_t(i)都相同，并且Y={-1,+1}。AdaBoost算法是在一系列的迭代（假设迭代次数t=T）中重复利用给定的弱分类器，并根据样本分布D_t(i)做出较弱的判定(weak hypothesis)：h_t:X →{-1,+1}。在每一次判定后算法会由较弱的判定结果求出弱分类器的分类误差ε_t，并根据此误差ε_t调整权重分布D_t(i)，增加错误样本的权重，使得弱分类器进一步专注于训练集中的困难样本。这样，当T次迭代结束后，就得到了最终的分类器，其分类假设为：

其中h_t(x)为原弱分类器分类假设，α_t是由每次分类误差ε_t求得的权重系数。

步骤四：用训练好的Adaboost分类器进行精检测。分类器训练完成后，把候选点集中所有点的梯度方向直方图特征输入Adaboost分类器，落在分类超平面正方向一侧的点则认为是最终检测到的民航飞机目标，反之则认为是伪目标而排除掉。此外为了检测到更大的飞机，此步对原始输入图像也进行了两次降采样，采样率为1.30、1.80。同时为了使该发明对飞机的方向变化更加鲁棒，本发明对原图及降采样后的图像顺时针旋转22.5°，这样实际上就使本发明可以检测的方向范围由原来的8方向细化为16方向。

步骤五：为了定位检测到的民航飞机，本发明对步骤四中的检测结果进行了区域生长聚类，并把每一聚类区域的中心作为飞机的中心。至此，完成了民航飞机检测的工作。

其中，步骤二中所述的提取训练样本的梯度方向直方图特征，其具体实现过程为：

（1）对40*40的输入图像计算图像的梯度，得到40*40的梯度图像。

（2）选择梯度图像左上角16*16的部分作为第一个块区(block)，该块区(block)包含大小为8*8的胞元(cell)共2*2个。对于每个胞元(cell)，都统计出一个9维的梯度方向直方图，其中的方向间隔采用无符号梯度，采样间隔为180°/9=20°。这样，把2*2个胞元(cell)的梯度方向直方图组合起来，再用该块区的能量密度进行归一化，就得到该块区的2*2*9=36维的特征描述符，记为F₀。

（3）让块区(block)以8为步长在图像中遍历，每一步都会得到36维的特征描述符F_i，其中0≤i≤15为移动的步数。最后把F_i(0≤i≤15)组合起来，即求得到该样本图像的16*36=576维的特征描述符。

其中，步骤三中所述的训练Adaboost分类器，其具体实现过程为：

给定样本：(x₁,y₁),...,(x_m,y_m)，(x_i∈X,y_i∈Y={-1,+1})，

（1）初始化D_t(i)=1/m (1≤i≤m)，m为训练样本总数。

（2）t=1，开始迭代（共迭代T次）。

（3）利用D_t(i)，(1≤i≤m)训练弱分类器，并利用之进行较弱判定(weak hypothesis)：

h_t:X →{-1,+1}

计算其分类误差为：

${\mathrm{\ϵ}}_t={\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]=\underset{i:h_t\left(x_i\right)\≠y_i}{\mathrm{\Σ}}{D}_t\left(i\right).$ 其中 ${\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]$ 表示分类错误的概率。

（4）更新权重：

${\mathrm{\α}}_t=\frac{1}{2}\ln \left(\frac{1-{\mathrm{\ϵ}}_t}{{\mathrm{\ϵ}}_t}\right),$ $\begin{array}{c}{D}_{t+1}\left(i\right)=\frac{D_t\left(i\right)}{Z_t}\×\left\{\begin{array}{cc}{e}^{-\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)=y_i\\ e^{\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)\≠y_i\end{array}\right.\\ =\frac{D_t\left(i\right)\exp (-{\mathrm{\α}}_t{y}_i{h}_t\left(x_i\right))}{Z_t}\end{array}$

其中Z_t是一个归一化系数，以使得D_t+1满足成为概率分布的条件,α_t为中间值。

（5）经过T次迭代，得到最终分类假设：

$H\left(x\right)=\mathrm{sign}\left(\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\α}}_t{h}_t\left(x\right)\right)$

本发明采用的弱分类器为决策树，层数为3层，迭代的次数为T=30次。

3、优点及功效：本发明的优点是：针对高分辨率遥感图像背景复杂的特点，利用民航飞机的圆周频率滤波响应与梯度方向直方图特征，可以在虚警较少的前提下保证较高的检测率，总体效果优于其它方法。

通过圆周频率滤波可以滤掉高分辨率遥感图像中的平滑背景区域，而梯度方向直方图则可以排除房屋、道路等伪目标。本发明的检测率与虚警率计算公式如下：

在Google Earth上截取的20张5000*5000的检测图中，有10张检测率在90%以上，而虚警率低于0.1%。另外10张检测率在80%以上，虚警率低于0.5%，效果较好。

附图说明：

图1为本发明所述方法检测流程框图

具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步描述：

本发明在Visual Studio 2010编程环境下实现。计算机读取高分辨率遥感图像数据后，首先进行多尺度圆周频率滤波，然后训练基于梯度方向直方图特征的Adaboost分类器，并利用之进行精检测，最后，对分类结果进行区域生长聚类，得到最终检测结果。

本发明的流程见图1所示，计算机配置采用：Intel Core i7-2600K处理器，主频3.4GHz，内存8GB，操作系统为Windows 7Ultimate，编程环境为Visual Studio 2010。该民航飞机自动检测方法包括以下步骤：

步骤一、读入图像，并对其进行多尺度圆周频率滤波。计算机在Visual Studio 2010环境下读取高分辨率遥感图像数据，其分辨率约为1米。之后对输入图像进行降采样和多尺度圆周频率滤波，得到飞机目标的候选点集。

步骤二、提取训练样本的梯度方向直方图特征，用以训练分类器。

步骤三、训练Adaboost分类器。通过30次的迭代，训练出以3层决策树为弱分类器的Adaboost分类器，其流程如下：

（1）(x₁,y₁),...,(x_m,y_m),(x_i∈X,y_i∈Y={-1,+1})为训练样本，m为训练样本总数，初始化权重D_t(i)=1/m(1≤i≤m)。

（2）t=1，开始迭代（共迭代T次）。

（3）利用D_t(i)，(1≤i≤m)训练弱分类器，并利用之进行较弱判定(weak hypothesis)：

h_t:X →{-1,+1}

计算其分类误差为：

${\mathrm{\ϵ}}_t={\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]=\underset{i:h_t\left(x_i\right)\≠y_i}{\mathrm{\Σ}}{D}_t\left(i\right).$ 其中 ${\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]$ 表示分类错误的概率。

（4）更新权重：

${\mathrm{\α}}_t=\frac{1}{2}\ln \left(\frac{1-{\mathrm{\ϵ}}_t}{{\mathrm{\ϵ}}_t}\right),$ $\begin{array}{c}{D}_{t+1}\left(i\right)=\frac{D_t\left(i\right)}{Z_t}\×\left\{\begin{array}{cc}{e}^{-\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)=y_i\\ e^{\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)\≠y_i\end{array}\right.\\ =\frac{D_t\left(i\right)\exp (-{\mathrm{\α}}_t{y}_i{h}_t\left(x_i\right))}{Z_t}\end{array}$

其中Z_t是一个归一化系数，以使得D_t+1满足成为概率分布的条件,α_t为中间值。

（5）经过T次迭代，得到最终分类假设：

$H\left(x\right)=\mathrm{sign}\left(\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\α}}_t{h}_t\left(x\right)\right)$

步骤四、用训练好的Adaboost分类器对候选点集进行精检测。

步骤五、通过区域生长聚类，定位飞机目标，得到最终检测结果。

有益效果：

实验结果：为了验证本发明的有效性，我们使用该方法进行实验，取得了较好的检测效果。实验所用数据为来自于在Google Earth的软件截图，截图大小为5000*5000，分辨率为1米。在20张实验图像中，其有10张检测率在90%以上，而虚警率低于0.1%；另外10张检测率在80%以上，虚警率低于0.5%，效果较好。

从实验结果来看，我们发明的方法很好的解决了高分辨率遥感图像背景复杂的问题，特别是对于房屋、道路等伪目标，也可有效排除，它具有广阔的应用前景和实用价值。

Claims (3)

1.一种高分辨率遥感图像中民航飞机自动检测方法，其特征在于：它包括如下步骤：

步骤一：读入图像，并对输入图像进行圆周频率滤波，读入的遥感图像分辨率为1米；

所述圆周频率滤波就是计算图像中每一个像素点的圆周频率响应，并把响应值大于某一阈值的像素点作为飞机目标的候选点集保留下来；小于该阈值的点则认为是非飞机区域而排除掉；其计算过程如下：

设图像上有一坐标为(i,j)的点，f_k(k=0,1,...,N-1)表示以(i,j)为中心以r为半径的圆周上逆时针方向像素的灰度值所组成的一维数列，其中r为这个圆周的半径，N为这个圆周序列的点数，求这个一维数列的傅立叶变换的幅值f(i,j)，如下所示：

$f(i,j)={\left(\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_k\cos \frac{\mathrm{c\πk}}{N}\right)}^2+{\left(\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_k\sin \frac{\mathrm{c\πk}}{N}\right)}^2$

上式中f(i,j)为圆周上的数列f_k经过离散傅立叶变换后的幅值；c是常量，代表求离散傅立叶变换时正弦和余弦函数的周期数；由于飞机主要由机头、一对机翼和尾翼组成，当半径r的取值使得圆周扫过这4部分时，它们分别对应于f_k当飞机比背景亮的4个峰或当飞机比背景暗的4个谷，因此飞机区域的圆周序列f_k通常有4个峰和4个谷，即近似4周期的正弦信号；而当傅里叶变换中正弦或余弦函数的周期与f_k的周期相同时，该傅里叶变换的幅值f(i,j)的绝对值最大，即在飞机中心的圆周频率响应的幅值大，因此令常数c=8，代表傅里叶变换的正弦与余弦函数有四个周期；

选择圆周半径r为4像素，采样点数N为40个，常数c为8，这样，长度为30到50像素的飞机区域会得到相对较强的响应，平坦的背景区域的响应相对较弱；此外为了检测更大的飞机，还对输入图像进行两次降采样，采样率为1.30、1.80，使得大飞机也有较强的滤波响应，称为多尺度圆周频率滤波，还有一些复杂的背景如房屋、道路边缘以及与飞机形状近似的建筑也会有较强的响应，这些点将在后面分类器精检测的过程中被排除掉；

步骤二：提取训练样本的梯度方向直方图特征；因为欲利用分类器对圆周频率滤波结果进行精检测，因此这一步是为后面训练分类器服务；训练样本为在Google Earth软件上截取的大小为40*40的图片，用来提取特征；这里采用的是梯度方向直方图特征，特征提取的流程如下：

首先将样本图像计算图像梯度，之后再将梯度图像分成小的连通区域，此连通区域被称为胞元cell，然后统计出胞元cell中所有像素点的梯度方向直方图，最后把这些直方图组合起来便构成特征描述符；为了使特征描述符对光照变化和阴影干扰具有更强的不变性，算法将这些局部直方图在图像的更大的范围内进行对比度归一化即contrast-normalized；归一化所采用的方法是：先计算各直方图在其所在块区中的能量密度，然后根据这个能量密度对区间中的各个胞元做归一化；

步骤三：训练Adaboost分类器；步骤二结束后，用提取的梯度方向直方图特征训练Adaboost分类器，由于梯度方向直方图特征不具有方向不变性，因此需要训练复数个Adaboost分类器以检测不同方向的飞机，这里训练了0°到315°等间隔的8个方向的分类器，训练样本为在Google Earth软件上截取的大小为40*40的图片；其中负样本，即非飞机图片5000幅；8个方向的正样本，即飞机图片，每个方向各有3000幅；

训练Adaboost分类器的流程如下：假定训练集为(x₁,y₁),...,(x_m,y_m)，训练样本的权重分布为D_t(i)(1≤i≤m)；其中x_i为样本特征，y_i为样本标签；假定起初时，所有的权重D_t(i)都相同，并且Y={-1,+1}；AdaBoost算法是在一系列的迭代中重复利用给定的弱分类器，并根据样本分布D_t(i)做出较弱的判定h_t:X →{-1,+1}在每一次判定后算法会由较弱的判定结果求出弱分类器的分类误差ε_t，并根据此误差ε_t调整权重分布D_t(i)，增加错误样本的权重，使得弱分类器进一步专注于训练集中的困难样本；这样，当T次迭代结束后，就得到了最终的分类器，其分类假设为：

其中h_t(x)为原弱分类器分类假设，α_t是由每次分类误差ε_t求得的权重系数；

步骤四：用训练好的Adaboost分类器进行精检测；分类器训练完成后，把候选点集中所有点的梯度方向直方图特征输入Adaboost分类器，落在分类超平面正方向一侧的点则认为是最终检测到的民航飞机目标，反之则认为是伪目标而排除掉；此外为了检测到更大的飞机，此步对原始输入图像也进行了两次降采样，采样率为1.30、1.80；同时为了对飞机的方向变化更加鲁棒，对原图及降采样后的图像顺时针旋转22.5°，这样就使检测的方向范围由原来的8方向细化为16方向；

步骤五：为了定位检测到的民航飞机，对步骤四中的检测结果进行区域生长聚类，并把每一聚类区域的中心作为飞机的中心，至此，完成了民航飞机检测的工作。

2.根据权利要求1所述的一种高分辨率遥感图像中民航飞机自动检测方法，其特征在于：步骤二中所述的提取训练样本的梯度方向直方图特征，其具体实现过程为：

（1）对40*40的输入图像计算图像的梯度，得到40*40的梯度图像；

（2）选择梯度图像左上角16*16的部分作为第一个块区，该块区包含大小为8*8的胞元cell共2*2个；对于每个胞元cell，都统计出一个9维的梯度方向直方图，其中的方向间隔采用无符号梯度，采样间隔为180°/9=20°；把2*2个胞元cell的梯度方向直方图组合起来，再用该块区的能量密度进行归一化，就得到该块区的2*2*9=36维的特征描述符，记为F₀；

（3）让块区以8为步长在图像中遍历，每一步都会得到36维的特征描述符F_i，其中0≤i≤15为移动的步数；最后把F_i(0≤i≤15)组合起来，即求得到该样本图像的16*36=576维的特征描述符。

3.根据权利要求1所述的一种高分辨率遥感图像中民航飞机自动检测方法，其特征在于：步骤三中所述的训练Adaboost分类器，其具体实现过程为：

给定样本：(x₁,y₁),...,(x_m,y_m)，(x_i∈X,y_i∈Y={-1,+1})，

（1）初始化D_t(i)=1/m(1≤i≤m)，m为训练样本总数；

（2）t=1，开始迭代，共迭代T次；

（3）利用D_t(i)，(1≤i≤m)训练弱分类器，并利用之进行较弱判定即weak hypothesis：

h_t:X →{-1,+1}

计算其分类误差为：

${\mathrm{\ϵ}}_t={\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]=\underset{i:h_t\left(x_i\right)\≠y_i}{\mathrm{\Σ}}{D}_t\left(i\right);$ 其中 ${\mathrm{Pr}}_{i~D_t}[h_t\left(x_i\right)\≠y_i]$ 表示分类错误的概率；

（4）更新权重：

${\mathrm{\α}}_t=\frac{1}{2}\ln \left(\frac{1-{\mathrm{\ϵ}}_t}{{\mathrm{\ϵ}}_t}\right),$ $\begin{array}{c}{D}_{t+1}\left(i\right)=\frac{D_t\left(i\right)}{Z_t}\×\left\{\begin{array}{cc}{e}^{-\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)=y_i\\ e^{\mathrm{\αt}}& \mathrm{if}{h}_t\left(x_i\right)\≠y_i\end{array}\right.\\ =\frac{D_t\left(i\right)\exp (-{\mathrm{\α}}_t{y}_i{h}_t\left(x_i\right))}{Z_t}\end{array}$

其中Z_t是一个归一化系数，以使得D_t+1满足成为概率分布的条件,α_t为中间值；

（5）经过T次迭代，得到最终分类假设：

$H\left(x\right)=\mathrm{sign}\left(\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\α}}_t{h}_t\left(x\right)\right)$

这里采用的弱分类器为决策树，层数为3层，迭代的次数为T=30次。

Images