基于超像素的协同分割变化检测方法

孙 扬， 朱 凌

(北京建筑大学 测绘与城市空间信息学院， 北京 100044)

在同一区域不同时间获取的遥感图像中，对地表变化进行定量分析和确定其变化的类型以及其空间分布情况，是遥感变化检测的过程. 在地表覆盖产品的生产、土地利用检测、地理国情监测和全球的资源环境监测等方面，都需要借助变化检测技术来实现，因此吸引了大量学者对变化检测技术进行研究.

传统的变化检测方法分为基于像素和基于对象两种[1]，然而这两种方法存在一定的局限性. 基于像素的变化检测对噪声敏感，易产生椒盐噪声[2]，只考虑单个像素光谱特征，不考虑像素之间的属性关系；基于对象的变化检测方法不具有普遍性，很难找到最佳分割尺度，而且形成多时相对象边界不一致[3]. 协同分割的概念最早是在计算机视觉领域出现的，其目的是为了快速准确地从一组图像中提取出共同包含的感兴趣目标，不需要人机交互即可完成图像分割. 这种分割技术十分适合从大量图像数据组中分割共同目标区域，广泛应用于各个科学领域. 为了解决传统变化检测的缺陷，我国学者袁敏等[3]将多视图图像协同分割思想引入到变化检测当中，提出了一种面向多时相高分辨率遥感图像变化检测的协同分割方法. 这种协同分割变化检测方法是在变化强度图的引导下，结合各时相的图像自身特征进行分割，通过能量函数的构建和优化，直接生成边界准确、空间对应的多时相变化对象. 谢振雷等[4]在此基础上改变了能量函数的优化方法，利用最大流/最小割的Dinic算法求得能量函数的最小割.

这种基于协同分割的变化检测方法相较于传统的检测方法有以下几点优势：

1)能够克服基于像素变化检测方法难以避免的“椒盐现象”；

2)与基于对象的检测方法相比，能够生成边界准确一致、空间对应的多时相变化对象.

但是协同分割变化检测的最小割/最大流算法，是将每一个像素视为一个节点来构建网络流图，导致计算次数与图中节点和边的数量直接相关. 如运行1 000×1 000像素以上的图像，节点和边的数量已是百万级，在此基础上再进行迭代计算会导致迭代次数过高，大大降低了该方法的计算效率.

超像素的概念最早是由REN等[5]于2003年提出的. 超像素即为具有相似的纹理、亮度和颜色的匀质相邻像素构成的不规则像素块. 超像素分割利用相邻像素之间的相似性将像素进行分组，利用少量超像素代替大量像素来表达图片特征，能够大幅度降低图像后期处理的复杂度，这种技术常用作分割算法的预处理. 因此，本文将超像素分割引入协同变化检测中，减少像素节点的数量，以此来改善其运算时间长，不能有效处理大范围影像的问题.

1 实验方法

图1为实验流程图，实验方法主要分为3步：数据预处理、超像素分割、协同分割变化检测.

图1 基于超像素的协同分割变化检测方法流程图Fig.1 Flow chart of base on superpixel cosegmentation change detection method

1.1 数据预处理

由于实验区的原始影像均为未做处理的原始数据，为了保证分割结果的质量，需要对两种影像分别进行预处理. 本研究中所有的预处理均在ENVI软件中完成.

1)正射校正：利用RPC参数文件以及DEM文件对高分一号卫星影像进行内定向、外定向和正射校正[6]，TM影像不需要正射校正.

2)几何校正：正射校正后，两时相影像的几何位置会有偏差，需要选择4个以上的控制点来进行配准，保证两期影像的误差小于一个像素即可.

3)大气校正：用高分一号辐射定标系数对图像进行辐射定标，定标计算公式为：

L=GainDN+Bias

(1)

式中：Gain为定标斜率，DN为卫星载荷观测值，Bias为定标斜率，L为定标后亮度值. 定标后进行FLAASH大气校正，消除大气、光照以及季节气候等因素对地物反射率的影响.

4)相对大气校正：由于两时相影像拍摄的相机不同，使得两者之间光谱值差异较大，需要利用线性矫正法对图像进行相对大气矫正[7].

1.2 超像素分割

本文所采用的超像素分割方法为Achanta等[8]于2010年提出的Slic(Simple Linear Iterative Cluster)分割方法. 该算法是将彩色图像转换到Lab颜色空间，并与xy坐标合并形成一个五维的特征向量，用简单的线性迭代对图像进行局部聚类，生成均匀的超像素[9]. 这种方法相较于其他的超像素算法的优越性在于：算法速度快，显著地减少了距离计算的数量，并且与超像素k的数量无关；可以通过修改参数来控制超像素的尺寸及其紧凑度. 算法具体步骤如下：

2)迭代聚类. 算法为了提高计算速度，以每一个种子点为中心建立一个2S×2S的搜索区域，计算每一个搜索区域的所有像素与种子点的距离度量D，取最短距离作为该像素的聚类中心. 完成全部种子点的搜索后，重新计算每个超像素块的质心，将并其作为下一次迭代的新种子点，直至迭代的次数大于某一个阈值即可. 距离度量D由颜色距离dc和空间距离ds构成，计算公式如下：

(2)

(3)

(4)

式(2)和式(3)中的L,a,b,x,y为用Lab颜色空间以及像素的xy坐标构建出的五维的特征向量，i代表每一个像素，k代表每一个种子点. 公式(4)中的m为控制超像素紧凑度的常数，m越大，得到的超像素就更加紧凑，m的取值范围为[1,40]，本研究中取值为10；S为相邻种子点之间的步长.

1.3 协同分割变化检测

协同分割变化检测部分的中心思想是将两幅或多幅图像之间的信息连通起来，利用这些信息构建能量函数. 然后将图像构建成网络流图，把能量函数的值视为网络流图中各类边的权重，利用最小割/最大流算法求得网络流图的最大流，根据Ford and Fulkerson定理[10]获得最小割集. 最后将最小割集中的边割开，即可获得最终的变化检测结果. 协同分割变化检测部分大致分为两部分，一是构建能量函数，二是用最小割/最大流方法优化能量函数.

1)构建能量函数. 本文的能量函数均遵循袁敏等[3]和谢振雷等[4]的能量函数构建方式，以此求得每一个超像素中所有像素的平均变化特征项和图像特征项，函数构建方式：

Energy=λE1+E2

(5)

式中：E1为变化特征项，E2为图像特征项，λ为变化特征项的权重，用来调整图像和变化特征参与分割的程度.

变化特征项E1是以变化强度图作为引导得到的，变化强度图中值越高的位置表示该区域的变化越明显[11].

图像特征项E2则是将图像中每个超像素块的光谱特征以及纹理特征加权组合得到的. 光谱特征是通过计算每一个超像素的所有光谱波段的平均值获得的. 纹理特征的提取方法采用灰度共生矩阵(GLCM)的纹理分析方法[12]，本文选取基于灰度共生矩阵的相关性(Correlation)、对比度(Contrast)、熵(Entropy)、角二阶矩(Angular Second Moment)4种特征来计算每一个超像素的纹理特征. 在袁敏等[3]的文献中，每一个时相的图像均会产生各自的图像特征. 将两个图像特征分别与变化强度图进行分割，则会产生对应于两时相图像的两个变化检测结果. 在本文中，通过赋予两时相图像特征一定的权重，来构成综合图像特征，利用综合图像特征进行协同分割只会产生一个变化检测结果. 这种方式的优点在于，构成的综合图像特征可以同时考虑两个时相的图像特征，而不是仅仅依赖于一个图像的特征进行协同分割.

E2=λt1Et1+(1-λt1)Et2

(6)

式中：Et1和Et2表示t1和t2时相图像的图像特征，λt1表示t1时相图像的图像特征所占的权重.

2.优化能量函数. 能量函数构建完成后，需要将能量函数最优化，从而获得最终的图像分割结果. 然而能量函数的优化是通过求得网络流图的最小割来实现的. 最小割/最大流算法是一种能量函数优化的方法，算法中将每一个超像素作为一个普通节点，然后设置两个特殊节点S(目标)、T(背景)，节点之间以边相连，将两个特殊节点分别与每个普通节点相连，普通节点与其相邻的普通节点相连(本文采用四邻域). 然后将变化特征项当作两个特殊节点与普通节点间边的权重，称为t连接；将图像特征项当作两个普通节点间边的权重，称为n连接. 由这种方式将超像素图像映射为网络流图，求取网络流图中的最大流.

根据Ford and Fulkerson定理[10]，在网络流图中最大流的值等于最小割的容量，最小割的表示公式：

(7)

cut为图中最小割集，cut是权值最小的边的合集，w(p,q)表示权值最小的超像素p、q之间的边. 若V1和V2分别为属于目标和背景的普通节点的子集，当超像素节点p、q之间的边在网络流图中权值最小时，超像素p、q分别属于两个不同的子集p∈V1和q∈V2，这条边即可划分进最小割集cut. 找到图中所有最小割集的边将其断开后，每一个超像素块都会被划分至V1或V2中，将V1子集中的超像素赋值为变化区域，V2子集赋值为非变化区域，即可完成协同分割变化检测.

本文利用基于最短路径的Dinic算法[13]求得图的最小割集. 将从节点S(目标)出发到节点T(背景)为止的路径距离固定为k，寻找到路径长度为k的路径后，进行流量更新，至少使路径中的一条边达到饱和. 将这条边划分至最小割集中，然后再次从源点出发，寻找新的路径长途为k的路径. 当没有路径长度为k的路径后，再重新进行分层，并寻找长度为k+1的路径,直到没有任何路径从原点到达汇点，算法结束. 然后将最小割集中的边断开，获得最终的协同分割变化检测结果.

2 试验结果和分析

2.1 研究数据介绍

本文的试验区共有2块，第一块选取中国江西省南昌市南昌县的高分一号卫星影像(图2). 本研究选取的实验区大小为1 350×1 350像素，两时相影像获取日期分别为2015-04-14和2017-04-29，空间分辨率为16 m. 影像共有4个波段，分别为蓝、绿、红3个可见光波段以及1个近红外波段. 从原影像中可以看出，该地区大部分的变化为人造覆盖的变化，新增了道路房屋和公共设施等.

图2 南昌县的高分一号卫星影像Fig.2 Remote sensing image of GF-1 satellite in Nanchang County

第二块实验区选取美国密西西比州杰克逊市的Landsat TM影像(图3)，实验区大小为1 200×1 200像素，面积约为1 296 km2，影像轨道号为P022R038，两时相影像获取日期分别为2000-09-17和2010-09-05，影像空间分辨率为30 m. TM影像共有7个波段，分别为蓝、绿、红、近红外、中红外、热红外以及追加的中红外波段. 实验区内发生明显变化的地类大部分为人类活动引起的变化，有小范围的人造覆盖变化.

图3 美国密西西比州杰克逊市Landsat TM影像Fig.3 Remote sensing image of Landsat in Jackson, Mississippi, USA

2.2 试验结果

文中的算法是在MATLAB软件中完成编写并得到最终变化检测结果. 图4为南昌实验区的变化检测结果，其中白色部分为未变化区域，黑色为变化区域. 为测试不同步长对试验结果的影响，将超像素分割步骤中种子点之间的步长设置为9和7，图4(a)(b)(c)为不加入超像素的结果. 为了更有效率的测试小步长结果的效果，在南昌试验区选取了一块450×450像素的试验区，将步长设置为3，其变化检测结果如图4(e)(f). 试验中加入和未加入超像素的协同分割变化检测结果，均使用的是同一种协同分割变化检测算法，即文中1.3小节中介绍的方法.

图5为杰克逊市试验区的变化检测结果，其中(a)、(b)是加入超像素分割后的变化检测结果，(c)是未加入超像素的变化检测结果，该试验区超像素分割步长分别为9和7.

图4 南昌实验区检测结果Fig.4 Results of Nanchang experimental area

图5 美国实验区检测结果Fig.5 Results of American experimental area

为了检验分割结果的准确性，本文在实验区内选择一定量的变化和非变化样本点，使用混淆矩阵在eCognition软件中进行精度检验. 混淆矩阵也称误差矩阵，是表示精度评价的一种标准格式，用n行n列的矩阵形式来表示[14]. 矩阵中的评价指标有总体精度、制图精度和用户精度，这些精度指标从不同的侧面反映了图像分类的精度. 在混淆矩阵中增加了可以用于一致性检验的Kappa系数来评定分割结果. 本文在实验区内选择一定量的变化和非变化样本，通过计算误差矩阵以及Kappa系数来对协同分割结果的准确性进行评价.

2.3 结果分析

首先，超像素分割的加入使得算法在运行时间上有极大幅度的提升. 由表1可看出，计算时间和超像素步长相关，最长不超过4 h，不需要将实验区分块处理，既保持了图斑的完整性又提高了算法运行速度，从而提高了算法的运行效率.

其次，分别对比两种分辨率影像不同步长的结果可以看出，两种分辨率影像的结果精度均是随着超像素步长的增加而降低. 从理论上分析，因为超像素简化了图像中的信息，所以使用超像素分割的结果精度应当低于未使用超像素分割的精度，但高分影像分割步长为3时，精度却高于未使用超像素的检测结果. 这是由于超像素的参与并且超像素的步长较小，导致原本图中误检的零散像素与周围地物合并，将零散像素聚类到了邻近的变化图斑中，从而提高了精度. 若步长设置过小，精度提高但运行效率并未提升，故而不可取.

然后，对比两种分辨率影像之间的结果可看出，高分影像的检测结果较好，基本保证了变化图斑的准确性，而TM影像检测结果准确度很低. 当固定了超像素分割步长后，高分影像对超像素分割的适应性较强，低中分辨率的影像在一个像素中包含的地物信息在分割前本就多于高分影像，在此基础上再用超像素简化图像信息，会造成地物信息的流失，使地物的边界不清晰.

表1 加入超像素前后结果对比

图6 2015年高分一号影像不同分割步长结果 对比(截取自试验区)Fig.6 Comparison of different segmentation step results of GF-1 images in 2015(Intercepted from the test area)

例如，图6、图7为两块试验区经过超像素分割后截取其中某一区域的结果，分割步长分别设置为5，7，9. 从图中可看出，当超像素步长是5或7时，高分影像分割得到的超像素块紧凑且均匀，在简化了图像信息的同时保留住影像中的地物基本特征. 而低中分辨率的影像则对超像素分割适应性较低，当步长为7或9时生成的图像十分模糊，几乎不能完整呈现地物原本的形状和纹理. 当步长为5时，这种现象有所缓解. 所以对于中低分辨率影像来说，步长设定在5较为适合. 对于高分影像来说步长设定在7较为适合，既保障了结果的精度，又提高了算法处理图像的速度.

图7 2010年TM影像不同分割步长结果对比 (截取自试验区)Fig.7 Comparison of different segmentation step results of TM images in 2010(Intercepted from the test area)

最后，无论步长为多少，变化检测出的变化地物边界与实际变化地物边界都会有不一致的情况，实际地物的边界是较为规整的，而检测得到的变化地物边界却呈现出凹凸不平的状态. 如图6的(b)(c)(d)，左下角的线状地物产生了变化，但却并未在变化检测结果中有所体现. 由于超像素分割模糊了地物原本的边界使其不再平整，甚至模糊掉了左下角的线状地物，导致精度降低.

综上所述，图像的分辨率和超像素的分割步长会对检测结果产生影响. 而引起变化检测结果精度降低的因素主要为：检测结果边界与不规则的变化地物边界不一致，变化范围较小的地物被忽略，如道路等现状要素或小范围的房屋建筑.

3 结论

本文针对协同分割变化检测方法效率不高，不能够一次处理大范围的影像这一缺陷，提出了基于超像素分割的协同分割变化检测方法. 超像素的引入使协同分割变化检测的运行速度得到了极大的提升，从超过1 d提升至最少不超过4 h. 与此同时变化检测的结果精度能够保持在0.8左右，保持了检测结果的准确性. 但这种变化检测方法的适用性会随着图像分辨率的下降而降低，所以更适用于高分辨率影像.

这种方法的优势有以下几点：

1)引入超像素分割减少了算法中像素节点的数量，在保证结果准确性的同时，大幅度提高了算法运行的效率；

2)结果不必对试验区分块并行处理，保证了图斑的完整性；

3)加入了超像素后，每一个超像素块中可以挖掘的信息大大增加，例如归一化植被指数、归一化水指数等特征指数、颜色特征、纹理特征以及形状特征等等，这些特征可以加入到能量函数中参与分割，进一步提高检测结果的精度.

同时这种方法也存在缺点，这种方法提取出的变化图斑的边界较粗，不能够很好地贴合原始影像中的变化地物的边界，导致精度下降. 因此可以继续进行的后续研究：

1)两时相影像经过分割后，得到的超像素块大小形状各异，可以将分割后的两期影像进行叠加分析，提取两期影像中不同的部分作为新的超像素加入其中，可以在一定程度上改善变化图斑边界与原始变化地物边界不贴合的问题.

2)由于超像素概念的引入，变化特征项以及图像特征项中可以挖掘超像素块中的信息，继续完善能量函数，以此来改善变化检测的精度.