CN104036295B - 多光谱遥感影像道路中心线自动检测方法 - Google Patents

Abstract

本发明公开了一种多光谱遥感影像道路中心线自动检测方法。其实现步骤是：1.输入多光谱遥感图像，调整其数据格式；2.对调整后的多光谱遥感图像进行道路光谱标记；3.计算多光谱遥感图像聚类中心及对应类别号；4.比较聚类中心和道路光谱标记的距离，获得道路类别号；5.将道路类别号对应像素点灰度值设为1，获得二值图像；6.滤除二值图像中细小斑点区域，获得去噪图像；7.滤除去噪图像中非道路特征区域，获得特征图像；8.滤除建筑物与道路耦合区域，获得区域图像；9.从区域图像中提取出道路中心线。本发明具有检测性能高、区域适应性强的优点，可用于郊区多光谱遥感影像中道路中心线和城区多光谱遥感影像中道路中心线的检测。

Description

多光谱遥感影像道路中心线自动检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像的目标检测，具体说是一种遥感图像道路目标检测方法。该方法适用于道路目标的自动检测。

背景技术

多光谱遥感影像是利用卫星多光谱扫描系统对地扫描获得的高分辨率图像，其丰富的光谱信息为地物影像计算机识别与检测提供了可能。但由于其光谱维度低，空间分辨率较低，使得土地利用、覆盖信息的提取成为遥感信息处理的难点问题之一。道路作为重要的人造地物是现代交通体系的主体，具有重要的地理、政治、经济意义，道路也是地图和地理信息系统中主要记录和标识对象。由于数字化地理交通信息的要求，迫切需要一种道路目标的自动检测算法。

Shukla V.和Chandrakanth R.在文章“Semi-Automatic Road ExtractionAlgorithm for High Resolution Image Using Path Following Approach”(ICV GIPSpace Applications Centre Ahmedabad，2002)中提出一种基于边缘跟踪的半自动道路检测方法。该方法先给定初始的种子点和初始方向，然后开始边缘跟踪，直到边缘终止处作为新的跟踪种子点，基于边缘跟踪方法得到的路径就是一条候选道路段。采用半自动检测算法得到的道路较为准确，但是需要人机交互提供给计算机先验信息，再由计算机完成剩余操作，其显著缺点是不能够自适应完成道路检测工作且人工干扰大、工作量大。

Thierry Geraud和Jean-Baptiste Mouret在文章“Fast Road NetworkExtraction in Satellite Image Using Mathematical Morphology and Markov RandomFields”(EURASIP Journal on Applied Signal Processing，2004)中提出了一种基于数学形态学和马尔科夫随机场的道路自动提取方法。该方法首先使用区域形态学算子对图像进行预处理，然后进行分水岭变换，保留道路信息，再用马尔科夫随机场模型对道路建模，按照道路节点的连接情况建立区域邻接图，并通过道路连接点数判断是否道路，从而去除非道路段。该方法对山区、乡村道路提取效果较好，但对于城市道路中建筑物与道路的粘连情况难以有效处理。

Zhang Q.和Couloigner I.在文章“A wavelet approach to road extractionfrom high spatial resolution remotely-sensed imagery(Geomatica,2004)中提出了一种多尺度道路提取方法。该方法运用小波变换提取道路中心线，根据道路交叉点和道路中心线在小波变换后不同层级的小波系数呈现出的不同规律来提取交叉点和中心线。该方法不太稳定，易受局部特征影响。

S.Idbraim和D.Mammass在文章“An automatic system for urban roadextraction from satellite and aerial images”(WSEAS Transactions on SignalProcessing archive，2008)中提出了一种高分辨率遥感影像道路提取算法。该方法采用基于自适应方向滤波器和分水岭变换，对图像沿道路方向进行逐块滤波，得到道路主方向后，将分水岭算法应用于梯度图上，纠正上一阶段的错误，计算分水岭盆地的表面曲率系数，判断是否道路。该方法存在的不足是由于道路边缘的建筑物和路面呈现近似的灰度特征，导致道路边界模糊，错误分割。

发明内容

本发明的目的在于针对上述现有道路检测算法的不足，提出了一种多光谱遥感影像道路中心线自动检测方法，以避免大量人工操作，提高道路检测精度和检测速度，实现对遥感影像中道路目标的自动检测。

实现本发明目的的技术方案包括如下步骤：

(1)输入高分辨率16位多光谱遥感图像，并调整这些遥感图像的数据格式构建光谱矢量SV；

(2)利用光谱矢量完善光谱标记库：

若光谱标记库未建立，则在调整后的遥感图像上人工标记N个典型光谱矢量V_i，i＝1…N,N为人工标记道路目标总数,V_i分别对应所有人工标记道路目标的典型光谱矢量，计算N个典型光谱矢量V_i的均值V作为道路目标光谱标记；

若光谱标记库已建立且标记库中包含道路目标的光谱标记，则取出光谱标记库中道路目标光谱标记V；

(3)设定聚类数目K，K≥3，采用K-Means聚类算法获得聚类后的聚类中心Q_i,以及Q_i对应的类别号R_i，i＝1…K；

(4)计算聚类中心Q_i与道路光谱标记V的欧氏距离Δ_i，i＝1…K，其中Δ_i的最小值为Δ_min，最小值Δ_min对应的类别号为R_min，类别号R_min为包含道路目标类别；

(5)将类别号R_min所对应光谱矢量的灰度值设定为1，剩余其它类别号R_i所对应光谱矢量的灰度值设定为0，i＝1…K，得到包含道路目标的二值图像BW；

(6)滤除包含道路目标的二值图像BW中连通域面积S小于50的细小斑点区域，得到去噪图像BW₂；

(7)滤除去噪图像BW₂中连通域长宽比C小于1.3的连通区域，并滤除该连通域矩形度P大于0.3的连通区域，得到满足道路形态学约束的特征图像BW₃；

(8)采用top-hat变换滤除特征图像BW₃中指定结构元素strel的区域，得到去除部分建筑物的耦合图像BW₄，该指定结构元素strel是指类型为圆盘型，大小为9的结构元素；

(9)采用多方向线形滤波方法完全滤除耦合图像BW₄中与道路粘连的建筑物区域，得到完全去除建筑物的区域图像BW₅；

(10)采用形态学细化算法对区域图像BW₅进行细化操作，得到道路中心线road。

本发明与现有方法相比具有如下优点：

第一，本发明较选用道路种子点和初始方向的半自动检测方法，减少了人工操作。

第二，本发明利用道路目标同其它地物具有不同光谱特征的性质，采用K-means算法获得道路检测结果，克服了现有道路检测方法对遥感光谱信息利用率低的问题。

第三，本发明利用一系列形态学处理获得最终道路中心线检测结果，由于形态学操作简单，处理速度快，使得本发明具有简单快速的特点。

第四，本发明利用道路目标在图像中呈现细长、线状特性并具有明显方向性的特点，使用top-hat变换和多方向线形滤波去除与道路目标耦合的建筑物区域，能很好的实现城市中道路目标的检测。

附图说明

图1为本发明的整体流程图；

图2为用本发明对郊区多光谱遥感图像的检测实验效果图；

图3为用本发明对城区多光谱遥感图像的检测实验效果图；

图4为用本发明对西安电子科技大学新校区的检测实验效果图。

具体实施方式

下面结合附图，对本发明的实现和效果作进一步的详细描述。

参照图1，本发明的实现步骤如下：

步骤1，输入多光谱遥感图像g，该多光谱遥感图像的行数为h，列数为w，光谱层数为p,其中h>0,w>0,p≥4。

步骤2，对多光谱遥感图像g的每个光谱层分别作累计直方图调整，得到调整后的多光谱遥感图像g′。

(2a)将多光谱遥感图像g中第k光谱层的子图像保存为g_k，k＝1…p，计算该子图像g_k的最小灰度值min_k和最大灰度值max_k，并在min_k～max_k范围内做子图像g_k的累计直方图；

(2b)将累计直方图2％位置处的灰度值保存为Zmin，将累计直方图98％位置处的灰度值保存为Zmax，将子图像g_k中所有小于Zmin的灰度值用Zmin代替，将子图像g_k中所有大于Zmax的灰度值用Zmax代替；

(2c)将子图像g_k在Zmin～Zmax范围内的灰度映射到0～65535灰度值，最后将各层子图像g_k按光谱层顺序拼合得到调整后的多光谱遥感图像g′。

步骤3，添加道路目标光谱标记V。

(3a)用调整后的多光谱遥感图像g′中的每个像素点在不同光谱层的灰度值构成矢量，即光谱矢量SV_q，q＝1…M，M表示多光谱遥感图像g′中光谱矢量的总数；

(3b)采用真彩色模式显示多光谱图像g′，在图像g′上人工标记N个最具代表性的典型光谱矢量V_i，对N个人工标记的典型光谱矢量V_i求取均值得到道路光谱标记V。

步骤4，采用K-Means聚类算法获得聚类后的聚类中心Q_i，以及Q_i对应的类别号R_i。

(4a)设置K-Means聚类算法的参数：选定距离测度为欧氏距离，选定收敛系数ζ为0.01，选定最大迭代次数T为50；

(4b)选择多光谱遥感图像g′的第一个光谱矢量SV₁作为当前光谱矢量d，随机选定K个初始聚类中心作为当前聚类中心Q_i；

(4c)计算当前光谱矢量d到所有当前聚类中心Q_i的距离，并将当前光谱矢量d归到距离最短的中心所在的类，得到当前光谱矢量d对应的类别号；

(4d)若当前光谱矢量d不是多光谱遥感图像g′的最后一个光谱矢量SV_M，则设置当前光谱矢量d的下一个光谱矢量为当前光谱矢量d，并重复步骤(4c)；

(4e)根据欧氏距离求每类的当前聚类中心Q_i；

(4f)重复(4b)到(4e)，直到所有中心的偏移都小于收敛系数ζ或达到最大迭代次数T时停止迭代，得到最终的聚类中心Q_i和Q_i对应的类别号R_i，i＝1…K，K≥3，K表示聚类数目。

步骤5，计算聚类中心Q_i与道路光谱标记V的欧氏距离Δ_i，i＝1…K，K≥3，K表示聚类数目，其中Δ_i的最小值为Δ_min，最小值Δ_min对应的类别号为R_min。

步骤6，将类别号R_min所对应光谱矢量的灰度值设定为1，剩余其它类别号R_i所对应光谱矢量的灰度值设定为0，得到包含道路目标的二值图像BW。

(6a)选择多光谱遥感图像g′的第一个光谱矢量SV₁作为当前光谱矢量d；

(6b)若当前光谱矢量d对应的类别号为R_min，则将当前光谱矢量d的灰度值设定为1，若当前光谱矢量d对应的类别号为其它类别号R_i，则将当前光谱矢量d的灰度值设定为0；

(6c)若当前光谱矢量d不是多光谱遥感图像g′的最后一个光谱矢量SV_M，则设置当前光谱矢量d的下一个光谱矢量为当前光谱矢量d；

(6d)重复步骤(6b)-(6c)，直到多光谱遥感图像g′的最后一个光谱矢量SV_M的灰度值为0或1，即得到二值图像BW。

步骤7，滤除二值图像BW中连通域面积S小于50的细小斑点区域，得到去噪图像BW₂。

采用K-Means算法获得的二值图像BW中常常包含细小斑点区域，这些斑点区域易对后续的道路检测造成干扰。为了得到更好的道路检测效果，本发明采用连通域方法滤除斑点区域，其步骤如下：

(7a)统计二值图像BW中所有初始连通域CR_j的面积S_j，j＝1…n，n为二值图像中初始连通域总数；

(7b)若初始连通域CR_j的面积S_j小于50，则将该初始连通域CR_j中所有像素点的灰度值设定为0，即得到去噪图像BW₂。

步骤8，滤除去噪图像BW₂中连通域长宽比C小于1.3的连通区域，并滤除该连通域矩形度P大于0.3的连通区域，得到满足道路形态学约束的特征图像BW₃。

(8a)统计去噪图像BW₂中所有去噪连通域CS_k的长宽比C_k，k＝1…m，m为去噪图像中去噪连通域总数，其中C_k＝l_k1/l_k2，l_k1、l_k2分别表示去噪连通域CS_k最小外接矩形的长和宽；

(8b)若去噪连通域CS_k的长宽比C_k小于1.3，则将该去噪连通域CS_k中所有像素点灰度值设定为0；

(8c)统计去噪图像BW₂中所有去噪连通域CS_k的矩形度P_k，k＝1…m，m为去噪图像中去噪连通域总数，其中P_k＝S_k1/S_k2，S_k1为去噪连通域CS_k的面积，S_k2为去噪连通域CS_k最小外接矩形面积；

(8d)若去噪连通域CS_k的矩形度P_k大于0.3，则将该去噪连通域CS_k中所有像素点灰度值设定为0，即得到特征图像BW₃。

步骤9，采用top-hat变换滤除特征图像BW₃中指定结构元素strel的区域，得到去除部分建筑物的耦合图像BW₄，该指定结构元素strel是指类型为圆盘型，大小为9的结构元素。

步骤10，采用多方向线形滤波方法完全滤除耦合图像BW₄中与道路粘连的建筑物区域，得到完全去除建筑物的区域图像BW₅。

(10a)构造多方向线形结构元素

其中，i表示多方向线形结构元素的序列号，i＝-30…30，α_i表示多方向线形结构元素的方向角，α_i＝i×3，L表示多方向线形结构元素的长度，L∈[70,100]，x_i表示多方向线形结构元素在图像中的横坐标，y_i表示多方向线形结构元素在图像中的纵坐标；

(10b)对耦合图像BW₄作多方向线形结构元素为的开运算，得到区域图像BW₅：

其中，∪表示求并集运算，ο表示求开运算。

步骤11，对区域图像BW₅进行形态学细化操作，得到道路中心线road。

(11a)将区域图像BW₅中道路目标的边界像素点的灰度值设定为0；

(11b)重复步骤(11a)，直到道路目标的像素宽度为1，即得到道路中心线road。

本发明的效果可以通过以下仿真实验做进一步的说明。

1.仿真条件

本发明是在中央处理器为Intel(R)Core i5-34703.20GHZ、内存8G、WINDOWS7操作系统上，运用MATLAB软件进行的实验仿真。

本发明仿真实验使用的多光谱遥感图像来自于SPOT-6卫星，有蓝(0.455-0.525μm)；绿(0.530-0.590μm)；红(0.625-0.695μm)；近红外(0.760-0.890μm)四个波段，空间分辨率为2m。

2.仿真内容

仿真1，用本发明对西安郊区多光谱遥感图像进行道路中心线检测，仿真结果如图2，其中图2a为西安郊区多光谱遥感图像的原图，图2b为西安郊区多光谱遥感图像的道路中心线检测图；

仿真2，用本发明对西安城区多光谱遥感图像进行道路中心线检测，仿真结果如图3，其中图3a为西安城区多光谱遥感图像的原图，图3b为西安城区多光谱遥感图像的道路中心线检测图；

仿真3，用本发明对西安电子科技大学新校区多光谱遥感图像进行道路中心线检测，仿真结果如图4，其中图4a为西安电子科技大学新校区多光谱遥感图像的原图，图4b为西安电子科技大学新校区多光谱遥感图像的道路中心线检测图；

从图2，图3和图4可见，本发明提出的道路中心线检测方法具有检测效果好、适应性强、检测速度快的优点。

综上，本发明可以有效检测郊区多光谱遥感影像中道路中心线和城区多光谱遥感影像中道路中心线。

Claims (6)

1.一种多光谱遥感影像道路中心线自动检测方法，包括如下步骤：

(1)输入高分辨率16位多光谱遥感图像，并调整这些遥感图像的数据格式构建光谱矢量SV：

(1a)将多光谱遥感图像g的行数保存为h、列数保存为w、光谱层数保存为p；

(1b)对多光谱遥感图像g各光谱层的图像分别做累计直方图调整，得到调整后的多光谱遥感图像g′：

(1b1)将多光谱遥感图像g中第k光谱层的子图像保存为g_k，k＝1,…,p，计算该子图像g_k的最小灰度值min_k和最大灰度值max_k，并在min_k～max_k范围内做子图像g_k的累计直方图；

(1b2)将累计直方图2％位置处的灰度值保存为Zmin，将累计直方图98％位置处的灰度值保存为Zmax，将子图像g_k中所有小于Zmin的灰度值用Zmin代替，将子图像g_k中所有大于Zmax的灰度值用Zmax代替；

(1b3)将子图像g_k在Zmin～Zmax范围内的灰度映射到0～65535灰度值，最后将各层子图像g_k按光谱层顺序拼合得到调整后的多光谱遥感图像g′；

(1c)用调整后的多光谱遥感图像g′中的每个像素点在不同光谱层的灰度值构成矢量，即光谱矢量SV；

(2)利用光谱矢量完善光谱标记库：

若光谱标记库未建立，则在调整后的遥感图像上人工标记N个典型光谱矢量V_i，i＝1,…,N,N为人工标记道路目标总数,V_i分别对应所有人工标记道路目标的典型光谱矢量，计算N个典型光谱矢量V_i的均值V作为道路目标光谱标记；

若光谱标记库已建立且标记库中包含道路目标的光谱标记，则取出光谱标记库中道路目标光谱标记V；

(3)设定聚类数目K，K≥3，采用K-Means聚类算法获得聚类后的聚类中心Q_i,以及Q_i对应的类别号R_i，i＝1,…,K；

(4)计算聚类中心Q_i与道路目标光谱标记V的欧氏距离Δ_i，i＝1,…,K,其中Δ_i的最小值为Δ_min，最小值Δ_min对应的类别号为R_min，类别号R_min为包含道路目标类别；

(5)将类别号R_min所对应光谱矢量的灰度值设定为1，剩余其它类别号R_i所对应 光谱矢量的灰度值设定为0，i＝1,…,K，得到包含道路目标的二值图像BW；

(6)滤除包含道路目标的二值图像BW中连通域面积S小于50的细小斑点区域，得到去噪图像BW₂；

(7)滤除去噪图像BW₂中连通域长宽比C小于1.3的连通区域，并滤除该连通域矩形度P大于0.3的连通区域，得到满足道路形态学约束的特征图像BW₃；

(8)采用top-hat变换滤除特征图像BW₃中指定结构元素strel的区域，得到去除部分建筑物的耦合图像BW₄，该指定结构元素strel是指类型为圆盘型，大小为9的结构元素；

(9)采用多方向线形滤波方法完全滤除耦合图像BW₄中与道路粘连的建筑物区域，得到完全去除建筑物的区域图像BW₅；

(10)采用形态学细化算法对区域图像BW₅进行细化操作，得到道路中心线road。

2.根据权利要求1所述的多光谱遥感影像道路中心线自动检测方法，其中所述步骤(2)中人工添加光谱标记，是采用真彩色模式显示多光谱图像，对于需要添加的道路目标，在图像上人工标记N个最具代表性的典型光谱矢量，并对这N个最具代表性的典型光谱矢量V_i求取均值得到道路目标光谱标记V。

3.根据权利要求1所述的多光谱遥感影像道路中心线自动检测方法，其中所述步骤(3)采用K-Means聚类算法获得聚类后的聚类中心Q_i，以及Q_i对应的类别号R_i，按如下步骤进行：

(3a)初始化：选定距离测度为欧氏距离，选定收敛系数ζ为0.01，选定最大迭代次数T为50，随机选定初始聚类中心；

(3b)对调整后的多光谱遥感图像g′的每个光谱矢量，求其到所有聚类中心的距离，并将该光谱矢量归到距离最短的中心所在的类；

(3c)根据欧氏距离求每类的中心；

(3d)重复(3b)、(3c)，直到所有中心的偏移都小于收敛系数ζ或达到最大迭代次数T时停止迭代，得到最终的聚类中心Q_i和Q_i对应的类别号R_i，i＝1,…,K。

4.根据权利要求1所述的多光谱遥感影像道路中心线自动检测方法，其中所述步 骤(6)中滤除二值图像BW中连通域面积S小于50的细小斑点区域，按如下步骤进行：

(6a)统计二值图像BW中所有初始连通域CR_j的面积S_j，j＝1,…,n,n为二值图像中连通域总数；

(6b)若初始连通域CR_j的面积S_j小于50，则将该连通区域中所有像素点的灰度值设定为0。

5.根据权利要求1所述的多光谱遥感影像道路中心线自动检测方法，其中所述步骤(7)中滤除去噪图像BW₂中连通域长宽比C小于1.3的连通区域，并滤除该连通域矩形度P大于0.3的连通区域，按如下步骤进行：

(7a)统计去噪图像BW₂中所有去噪连通域CS_k的长宽比C_k，k＝1,…,m，m为BW₂中去噪连通域总数，其中C_k＝l_k1/l_k2，l_k1、l_k2分别表示去噪连通域CS_k最小外接矩形的长和宽；

(7b)若去噪连通域CS_k的长宽比C_k小于1.3，则将该连通域区域中所有像素点灰度值设定为0；

(7c)统计去噪图像BW₂中所有去噪连通域CS_k的矩形度P_k，k＝1,…,m，m为BW₂中去噪连通域总数，其中P_k＝S_k1/S_k2，S_k1为去噪连通域CS_k的面积，S_k2为去噪连通域CS_k最小外接矩形面积；

(7d)若去噪连通域CS_k的矩形度P_k大于0.3，则将该连通域区域中所有像素点灰度值设定为0。

6.根据权利要求1所述的多光谱遥感影像道路中心线自动检测方法，其中所述步骤(9)中采用多方向线形滤波方法完全滤除耦合图像BW₄中与道路粘连的建筑物区域，按如下步骤进行：

(9a)构造多方向线形结构元素

其中，i表示多方向线形结构元素的序列号，i＝-30,…,30, α_i表示多方向线形结构元素的方向角，α_i＝i×3，L表示多方向线形结构元素的长度， L∈[70,100]，x_i表示多方向线形结构元素在图像中的横坐标，y_i表示多方向线形结构元素在图像中的纵坐标；

(9b)对耦合图像BW₄作多方向线形结构元素为的开运算，得到区域图像BW₅：

其中，∪表示求并集运算，表示求开运算。