基于多级框架显著性特征的遥感影像道路识别方法

本发明以遥感影像为研究对象，公开了一种基于多级框架显著性特征的遥感影像道路识别方法，涉及图像处理领域，具体涉及遥感处理识别领域。

近年来，随着传感器、遥感平台、数据通讯等技术的飞速发展，遥感技术进入了一个能够动态、快速、准确、及时、并用多手段提供各种对地观测数据的新阶段，为人们获取环境信息、认识环境提供了一个重要途径。

遥感影像是从远离地面的工作平台上，通过传感器对地球表面的电磁波信息进行探测，然后再经过信息传输、处理与分析，最终得到遥感影像。它具有很高的空间分辨力和光谱分辨力，可以在几百个甚至上千个连续光谱波段获取地物目标影像。与此同时，随着智慧城市的提出，道路作为重要的人造地物是现代交通体系的主体，具有重要的地理、政治、经济意义，道路也是地图和地理信息系统中主要的记录和标识对象。上世纪70年代中期，由于数字化地理交通信息的需要，道路影像自动提取技术随之出现并逐步发展。如今，多光谱高分辨率的遥感卫星、成像雷达、无人驾驶飞机的出现，使得对地观测手段更加完备，地理影像数据日益丰富。另一方面，地图绘制、地理信息系统更新、城市观测和规划等诸多应用需求的出现和不断增长，促使自动道路提取技术不断发展。因此，遥感影像中的道路检测与识别成为遥感领域的研究热点之一。

基于显著性的遥感影像识别方法，依据影像内容的显著性建立适用于遥感影像的显著性模型，对其光谱特征及纹理、结构等特征进行综合分析，能够在复杂的影像环境中快速定位识别区域，在一定程度上消除了影像的冗余信息，突出了影像的主要内容，降低了影像处理分析的复杂度并减小语义鸿沟。能在没有先验信息的条件下检测与周围背景存在差异的目标(如自然背景中的人造目标)。不需要背景或目标的先验光谱信息，实用性较强，为遥感影像目标识别、影像分类与检索奠定基础。为此，本发明将显著性特征引入遥感影像的道路识别中，提出了 一种基于多级框架显著性特征的遥感影像道路识别方法。

本发明与已有的遥感影像道路识别方法不同，针对遥感影像中道路区域显著性特征明显的特点，对传统的遥感影像识别技术进行改进，提出了多级框架显著性特征的道路识别方法。该发明改进传统的显著性特征，对遥感影像进行光谱差异和局部线性轨迹的多级框架显著性特征提取，再对所得到的道路网进行优化：采用后处理技术去除道路网中的孤立斑块和噪声，采用部分区域分割(Region Part Segmentation,RPS)方法消除非道路区域，通过中轴变换(Medial Axis Transform,MAT)方法去除类似道路光谱特征的区域，例如停车场、大型建筑的屋顶。最后采用基于领域覆盖方法的增量学习对道路网进行训练，实现遥感影像的道路识别。主要过程如附图1所示，可分为以下几个步骤：基于多级框架显著性特征的道路网络提取、基于多级框架显著性特征的道路网络优化和基于增量学习的道路种类识别。

步骤1：多级框架显著性特征的提取

步骤1.1：基于主导奇异度量方法的局部线性显著性特征提取

采用梯度矩阵奇异值分解(Singular Value Decomposition,SVD)的方法计算像素点的主导方向算子，然后设置主导方向阈值，去除不具有主导方向的背景区域，实现道路边缘像素的提取；

步骤1.1.1：梯度矩阵计算

对于影像中的每个像素，首先计算该像素点周围N×N大小窗口的梯度矢量，影像f(x,y)中的点(xk,yk)的梯度通过下式计算：

${\mathrm{GF}}_k={[\frac{\mathrm{\delta f}(x_k,y_k)}{\delta \left(x\right)},\frac{\mathrm{\delta f}(x_k,y_k)}{\delta \left(y\right)}]}^T$

其中 <math><mrow> <mi>&delta;f</mi> <mo>/</mo> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>f</mi> <mi>y</mi> </msub> <mo>&CircleTimes;</mo> <mi>&delta;G</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>&delta;f</mi> <mo>/</mo> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>f</mi> <mi>x</mi> </msub> <mo>&CircleTimes;</mo> <mi>&delta;G</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>,</mo> <msub> <mi>f</mi> <mi>y</mi> </msub> <mo>=</mo> <mi>f</mi> <mo>&CircleTimes;</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>,</mo> <msub> <mi>f</mi> <mi>x</mi> </msub> <mo>=</mo> <mi>f</mi> <mo>&CircleTimes;</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>,</mo> </mrow></math> 表示一维卷积；

步骤1.1.2：梯度矩阵的奇异值分解

得到梯度矩阵以后，计算梯度矩阵SVD的协方差公式如下：

$\mathrm{Cov}\left(G\right)={\Sigma }_i{\mathrm{GF}}_{i·}{\mathrm{GF}}_i^T=\left[\begin{array}{cc}{s}_1& s_2\\ s_3& s_4\end{array}\right]$

其中s1,s2,s3,s4为梯度矩阵的奇异值分量，s1为主方向奇异值；

步骤1.1.3：主导奇异算子的计算

主导奇异算子(Dominant Singular Operator,DSO)用来计算主方向奇异值和奇异值总和的比率，当该比率接近1时说明该像素点的周围区域为细长形状；因此，DSO定义如下：

$\mathrm{DSO}=\frac{s_1}{\underset{i=2}{\overset{4}{\Sigma }}{s}_i},s_1>s_2,s_3,s_4>0$

由上式可以看出，当所有的梯度分量具有相同方向时，只有一个奇异值是非零的，即DSO的值为1；如果四个奇异值相等，则DSO的值为0.25，因此DSO值的范围是[0.25,1]；本发明定义一个阈值，当DSO的值小于该阈值时，说明对应的影像块是嘈杂的，没有主导方向，没有体现出道路局部线性的特征，因此被视为背景区域即非道路区域；相反则被看作是道路区域；

步骤1.2：基于概率支持向量机的光谱差异显著性特征提取

首先将遥感影像进行分块处理，然后提取图像块的光谱特征，并用特征向量表示，最后通过概率支持向量机计算像素点的概率标签；

采用Tile分块方式对遥感影像进行预处理，将遥感影像库的目标影像从左至右、从上至下分成规则、不重叠且大小相等的子图像块；

提取图像块的光谱特征，选取训练样本，训练集分成非道路子影像和道路子影像块，对于每一个子影像样本，在红、绿、蓝三个波段内提取它的直方图作为特征矢量，包含平均值、能量、标准差和熵，每个样本用特征矢量表示，并计算分块以后的遥感数据，再通过SVM分类器训练这些样本的特征矢量；

通过概率支持向量机获取识别样本的概率P：

$P(y=1|f_s)=\frac{1}{1+\exp (A{f}_s+B)}$

其中y是二分类标签，fs是从SVM决策函数得到的标签，参数A和B通过最大似然估计训练集获得，对于每一个影像块，如果P值大于所规定的阈值α，那么认定该样本为道路区域，反之则认定为非道路区域；

步骤1.3：基于赢者取全神经网络的多级框架显著性特征融合

首先将道路边缘显著图记为Sl，道路区域显著图记为Sg，对二者进行归一化合并即可得到总的显著图S，计算公式如下:

$S=\frac{1}{2}[N\left(S_l\right)+N\left(S_g\right)]$

其中N(·)为归一化因子，归一化算子N(.)是通过以下三个步骤进行计算：1)将每一张特征图的灰度范围归一化至某个特定范围[O,M]，以消除不同特征模式下的幅度差异；2)出每一张特征图的最大值M以及所有其他局部最大值的平均值 3)将特征图都乘以通过使用N(.)算子进行归一化，使得不同特征图上的局部较大值得到增强，而那些局部最大值不明显的特征图得到抑制，从而使得显著的区域得到突出，而均匀的非显著区域被忽略；

然后通过WTA神经网络计算每一个像素点的概率标签，神经元的概率标签在每一次迭代后采用“赢者取全”的方式进行更新，若该神经元的第k类为优胜者，则对第k类的概率增加一个常数；反之则减少一个常数；

WTA方法反复迭代地增加道路区域像素和删除非道路区域像素，因此本质上是一个校正机制的融合过程，本发明方法需要道路边缘和区域特征作为输入，在融合之后给予像素一个最终的标签道路或非道路，这个标签作为WTA神经元的输出，像素点属于两个类的概率：道路k＝1，非道路k＝2，约束为P(1)+P(2)＝1；比较两个类的概率，像素点被确定为道路或非道路，最终实现道路边缘和道路区域融合；

步骤2：多级框架显著性特征的优化

步骤2.1：道路网的后处理

本发明采取连通区域标记的方法保留不相交的网络线段，去除非道路区域；计算每个连通区域的面积，删除小于面积阈值的区域，同时计算每一个联通区域的偏心距，删除小于偏心距阈值的区域，偏心距的计算公式如下：

${\mu }_{20}=M_{20}-\bar{x}{M}_{10}$

${\mu }_{02}=M_{02}-\bar{y}{M}_{10}$

E＝μ20/μ02

其中 $\bar{x}=\frac{M_{10}}{M_{00}};\bar{y}=\frac{M_{01}}{M_{00}};M_{\mathrm{pq}}=\underset{x}{\Sigma }\underset{y}{\Sigma }{x}^i{y}^{i}I(x,y);$

步骤2.2：基于部分区域分割方法的道路网优化

采取部分区域分割方法消除道路边缘的非道路区域，步骤如下：

步骤2.2.1：对道路网络的内外轮廓线进行平滑处理；

步骤2.2.2：计算平滑轮廓的曲率；将一条曲线表示为参数形式，其中t表示路径长度，x和y是轮廓的坐标：

r(t)＝(x(t),y(t))

曲率被定义为斜率的变化率，通过以下公式计算：

$K_s\left(t\right)=\frac{x_s^{\prime }{y}_s^{\prime \prime }-y_s^{\prime }{x}_s^{\prime \prime }}{{(x_s^{\prime 2}+y_s^{\prime 2})}^{3/2}}$

步骤2.2.3：确定局部极值，即局部曲率的导数为0的点；

$K_s^{\prime }\left(t\right)=\frac{d{K}_s\left(t\right)}{\mathrm{dt}}=0$

步骤2.2.4：通过跟踪区域的外/内轮廓，将局部极值点作为凸/凹主导点(Convex/Concave Dominant Points,CDPs)，其中外轮廓的凸主导点设为CDPcx，内轮廓的凹主导点设为CDPce；

步骤2.2.5：CDPcx沿其法线方向进行内向运动，CDPce沿其法线方向进行外向运动，当遇到同一个轮廓线内另外一个移动的CDP时停止；

1)追踪所有冻结的CDPs的起点，并且连接相应的CDPs，该部分区域将被 去除；

步骤2.3：基于中轴线变换理论的道路网优化

本发明提取道路的中轴线，然后将道路网中小于平均宽度的部分保留，删除宽度过大的非道路区域，具体的步骤如下：

步骤2.3.1：中轴线骨架的确定：首先确定图像中的边缘轮廓，平滑边缘像素点；然后跟踪一条边缘轮廓的每一像素点，到它所对应的另一条边缘的关联点；再连接关联点，将连接线的中点作为中轴点；最后连接整条边缘轮廓的所有中轴点，得到图像的中轴骨架；

步骤2.3.2：非道路区域的消除：得到中轴线骨架以后，计算每一部分的属性，包含该段的宽度、长度以及和标准宽度的标准差；然后再从道路网的源节点开始遍历，将不满足假设的部分去除，即去除平均宽度大于标准宽度或标准差过大的区域；

步骤3：基于增量学习的道路种类识别

步骤3.1：道路网络特征的选取

根据高速公路、城市道路、山区道路的特点，提取道路网的特征因子：

面积反映道路所占区域的大小；细长度为区域外界矩形的长宽比，能够反映道路的宽度；交叉点个数统计道路网中的交叉点个数；曲率表示道路网的完全程度；纹理特征为利用灰度共生矩阵，计算能量、灰度相关、局部平稳、熵和惯性矩5个纹理特征系数；

然后将这些特征因子向量化得到特征向量，最后采用增量学习方法对特征向量进行训练和识别；

步骤3.2：基于增量学习的特征识别

本发明采取基于覆盖的增量学习方法对图像特征进行训练和识别，基于增量学习方法的训练与识别方法包括：(1)对初始样本的特征进行训练，生成基础分 类器，构建初始分类模型；(2)对当前一批新增数据进行识别，根据识别结果对已有分类器进行调整，应用调整后的分类器再对下一批新增数据进行识别；(3)如此反复，每次依据当前识别结果调整已有分类器，直至所有新增数据识别完为止；

步骤3.2.1：初始分类模型的建立

本发明采用领域覆盖方法训练初始样本，为基于增量学习的道路影像识别构建初始分类模型，方法如下：

1)领域覆盖方法的初始化

首先对道路网络样本进行初始化：

(1)设道路网络样本集K＝{(x1,y1),(x2,y2),...,(xp,yp)},其中K中有p个样本，xi为道路网络的特征向量，每个输入样本xi(i＝1,2,...,p)有n维属性，yi为xi对应道路网络的种类；

(2)设p个样本对应的输出中不同的yi有q个，令I(t)(t＝1,2,...,q)分别代表令样本的输出为yt的所有样本标号的集合；将对应的输入集合记为P(t),P(t)＝{xi|i∈I(t)}；

2)领域覆盖方法构建初始分类模型

在对道路网络样本初始化以后，领域覆盖方法通过以下步骤构建初始分类模型：

(1)设道路网络训练样本集为X＝{(xt,yt),t＝1,2,...,p}，样本共分为s类；

(2)求样本集X中样本的最大模r，如(19)所示；然后将X中的所有样本点依次按公式(20)变换，并投射到中心在原点、半径为R(R＝r+1)的超球面上，

其中T(xk)为样本点投射后的新坐标；

r＝max{|xk||xk∈X}

$T\left(x^k\right)=(x^k,\sqrt{R^2-{\left|x^k\right|}^2})$

(3)求投射后样本集X所对应的I(t),P(t)(t＝1,2,...,s)；

(4)随机选取P(t)中未被覆盖的一个样本点xk,令

$d^1\left(k\right)=\underset{m\notin I\left(t\right)}{\max }\{⟨x^k,x^m⟩\}$

$d^2\left(k\right)=\underset{m\in I\left(t\right)}{\max }\{⟨x^k,x^m⟩|⟨x^k,x^m⟩>d^1\left(k\right)\}$

$d\left(k\right)=\frac{d^1\left(k\right)+d^2\left(k\right)}{2}$

其中〈xk,xm〉表示xk和xm进行内积运算，d1(k)表示与xk异类点内积的最大值，d2(k)表示与xk同类点内积的最小值；

(5)作以x^k为球形领域中心、半径为θ＝d(k)的覆盖 其中表示第i类样本的第j个覆盖；

(6)将P(t)中已被覆盖的点做上标记，判断训练样本集合X中是否所有的点均被标记；若是，本发明方法结束；若否，继续判断P(t)内样本点是否全部被覆盖并标记，i＝i+1、j＝1、t＝t+1,返回(4)，若否，令j＝j+1，返回(4)；

本发明方法最终求出一组覆盖 $C=\{C_1^1,C_2^1,...,C_{n_1}^1,C_1^2,C_2^2,...,C_{n_2}^2,...,C_1^s,...,C_{n_s}^s\},$ C即为训练初始道路网络特征所得到的覆盖模型；

步骤3.2.2：新增样本的识别

基于覆盖的增量方法步骤如下：

(1)对初始的道路网络样本进行学习，构建初始分类模型，得到一组覆盖C；

(2)用已有模型对当前一批新增的道路网络数据Rn进行识别；

(3)根据识别结果，对已有的分类模型进行调整，得到新的模型，方法如下：首先从Rn中到识别错误的样本，并把初始模型中错误样本对应的覆盖Cm删除；然后将识别错误的样本和拒识样本重新进行训练，得到新的覆盖Ca；最后将新得到的覆盖Ca加入到原有的分类模型当中，至此，得到新的分类模型Cn，且Cn＝C-Cm+Ca；

(4)重复(2)(3)直至所有新增样本识别完成。

本发明与现有技术相比，具有以下明显的优势和有益效果：

本发明以遥感影像为研究对象，充分利用遥感影像中道路的光谱、物理以及形态特征，对传统的显著性特征进行改进，提出遥感影像的光谱差异和局部线性轨迹两种显著性特征，同时，对所提取的道路网络进行近一步分析，得出该路段所属类型。首先提取遥感影像的显著性特征，并将光谱差异和局部线性轨迹两种特征进行融合，建立多级框架显著图；然后采取RPS方法、MAT方法等对道路网络进行优化，近一步消除非道路区域；最后提取道路网络的特征，并采用增量学习方法对道路进行分类，实现多级框架显著性特征的遥感影像道路识别。

图1基于多级框架显著性特征的遥感影像道路识别方法流程；

图2道路边缘提取流程；

图3道路区域提取流程；

图4多级框架显著性特征融合流程；

图5中轴变换方法示意图；

图6基于增量学习的识别流程；

图7道路网络样本集与覆盖在二维空间的分布示意图；

图8多级框架显著性特征提取；

图9多级框架显著性特征优化。

根据上述描述，以下是一个具体的实施流程，但本专利所保护的范围并不限于该实施流程。本发明选取512×512大小的多光谱卫星影像作为实验样本，实验环境为2.10GHz CPU，5.00G内存的PC机，实现环境为Windows 7+Visual Studio 2010。

首先进行多级框架显著性特征的提取，采用DSM方法和P-SVM分别提取遥感影像的光谱差异和局部线性显著性特征，然后再采用“赢者取全”神经网络进行多级框架显著性特征融合，实现多级框架显著性特征的提取。实验结果如附图8所示，其中(a)为512×512大小的原始输入图像，(b)为DSM道路边缘提取结果，(c)P-SVM道路区域提取结果，(d)神经网络特征融合结果。

然后进行多级框架显著性特征的优化，分别采用图像后处理、部分区域分割方法和中轴线变换方法对道路网络进行优化，去掉道路网络中的非道路区域，进一步提高道路网的拓扑性和完整性。实验结果如附图9所示，其中(a)为图像后处理结果，(b)RPS方法处理结果，(c)MAT方法处理结果。