基于多尺度分割的高分辨率遥感影像镶嵌线自动提取

温银堂， 王铁柱， 王书涛， 王贵川， 刘诗瑜， 崔 凯

(燕山大学电气工程学院河北省测试计量技术及仪器重点实验室，秦皇岛 066004)

0 引言

随着我国航天事业的蓬勃发展，高分辨遥感影像凭借覆盖范围广、包含信息丰富、周期性观测等优势逐步走入大众的视线，在城市规划、抢险救灾、生态环境监测等方面发挥了重要作用[1]。由于卫星相机视角限制，单景的遥感影像覆盖范围有限，难以展现大范围地区的地物特征。为了获取更大范围、更加完整的高质量遥感影像数据，通常需要将多景影像进行镶嵌处理。遥感影像的镶嵌，就是将从传感器中获取的两幅或两幅以上的遥感影像数据根据一系列的遥感技术手段拼接在一起，形成一整幅大范围影像的过程[2]。镶嵌线的提取是其中非常关键的一步，为避免因不同影像间的灰度差异产生明显的接缝线，镶嵌线不能直接横穿建筑物、湖泊、农田等地物区域。在遥感影像制作过程中，选择合适镶嵌线是消除接缝痕迹的重要手段[3]。

目前常用的商业软件对于镶嵌线的提取结果并不理想，不能有效避让建筑物等明显的地物目标[4]，无法满足较高要求的镶嵌任务，常常还需要人工进行多次编辑修改，耗费大量人力物力，使用算法提取最优镶嵌线仍然是研究难点。现阶段镶嵌线提取方法主要可以分为3种，分别是基于影像重叠区域差异的方法[5]、基于辅助数据的方法[6]和基于形态学的方法[7]。基于重叠区影像差异的方法通常对重叠区域像素对比计算差值，利用一定搜索策略选择最优镶嵌线[8]，方法简单容易实行，但计算量通常较大，重叠区域差异主要是针对单个像素或局部纹理，镶嵌线难免穿过房屋建筑等区域。基于辅助数据的提取方法需要道路矢量、激光雷达(light detection and ranging,LiDAR)点云等数据配合，能够精准有效避开目标地物[9]，镶嵌线质量较高，但辅助数据的获取有一定的难度。基于形态学的方法使用影像分割手段将重叠区域划分成不同的区域，然后利用寻优算法选择镶嵌线通过区域，使接缝线沿着地物本身的边界实现自然过渡，如韩天庆[10]使用标记的分水岭算法和Dijkstra算法寻找最佳接缝线，但分割过程可能产生过分割现象; 岳桂杰等[11]使用Canny边缘检测和A*算法，避开边缘图像的边缘信息寻找最短路径，边缘特征不能很好地反映地物完整性。本文使用超像素合并的方法对遥感影像进行多尺度分割，以保留影像地物纹理信息，使用A*算法在分割线上进行路径寻优，有效减少镶嵌线穿过明显地物，确保地物完整性，对镶嵌线的提取具有重要意义。

1 研究方法

首先对影像进行正射校正等必要的预处理。校正后影像四周会出现无效的黑色区域，通过采样判断有效区域的范围并生成每张影像的有效区域多边形，并以多边形交点作为镶嵌线的起始点。利用简单线性迭代聚类(simple linear iterative cluster，SLIC)算法对于重叠区域影像进行预分割，使颜色、亮度位置等特征相似的像素聚集成紧密的像素块，然后通过区域连接图对相邻相似区域进行合并生成不同尺度的数据集，减少过分割。完成分割后对影像进行掩模，对于分割区域内部设置为障碍，分割线作为可行区域，利用A*算法完成重叠区域从起点到终点的搜索，生成镶嵌线，具体过程如图1所示。

图1 基于多尺度分割的镶嵌线自动提取技术Fig.1 Automatic extraction of mosaic line based on multi-scale segmentation

1.1 SLIC超像素分割

使用SLIC算法对重叠区域影像进行预分割，对明显地物区域进行聚类生成紧密的超像素。超像素是由空间位置相邻，颜色、亮度和纹理等特征相似的像素点组成的图像子区域。SLIC分割能有效避免椒盐噪声带来的影响，可以更好地保留图像地物的纹理信息[12]。生成具有较好边界附着性且具有一定视觉意义的像素块,分割结果的边界紧凑度较高，且是当前效率最高的超像素分割算法[13]。首先要将图像从RGB颜色空间转到Lab颜色空间，图像每个像素对应的(l,a,b)颜色值和所处的位置坐标(x,y)共同构成一个5维的特征向量[14]，根据特征向量将相似的像素聚类生成为超像素块，具体算法步骤[15]：

1)设置生成超像素的个数K，生成K个种子点。将尺寸为M×N的影像分割为K个像素块，每个像素块的大小就是M×N/K，邻近的两个种子的距离为：

(1)

2)调整生成的种子位置。为了防止初始种子点位落在影像物体的轮廓边缘或噪声点上，计算种子点3×3窗口邻域所有像素的梯度值，把初始种子点位置调整到梯度值最小的像素位置，图像梯度计算如下：

G(x,y)=‖I(x+1,y)-I(x-1,y)‖2+
‖I(x,y+1)-I(x,y-1)‖2

(2)

式中：I(x,y)为(x,y)点的Lab向量； ‖‖为求欧几里得范数。

3)计算距离迭代聚类。为了增加算法的运算效率，SLIC只对种子2S×2S的邻域范围内的像素进行搜索，计算像素点与种子点之间的距离，将与种子点距离最近的若干个像素归于一类，距离由颜色距离和空间距离共同决定，计算公式如下：

(3)

(4)

(5)

式中： [lk,ak,bk,xk,yk]为种子点的特征向量； [li,ai,bi,xi,yi]为像素块内待判断的像素点的特征向量；dlab为两个像素点在Lab颜色空间内欧式距离；dxy为两个像素在空间位置上的欧式距离；D为待测像素点和种子点的总距离，可以表示两个像素的相似度；S为种子点的间距大小；m为紧凑度因子，调整颜色信息和距离信息的相对重要程度，值越大则越强调空间的邻近性，在这里取经验值10。SLIC算法种子点搜索范围限制在2S×2S之内，可以减少不必要的计算，有效地提高收敛速度。

所有像素归类完毕后计算超像素里所有像素点的平均向量值，把区域内最接近平均向量值的点调整为聚类中心点，然后重复迭代2)— 3)步骤，直到收敛或达到迭代次数停止计算。多数情况下，10次迭代对于绝大多数影像都能取得比较好的分割效果[16]。对于迭代聚类后的影像进行遍历筛选，其中会有未聚类的孤立点和面积过小的像素点，直接加入到就近的超像素中，以增强连通性。

1.2 基于区域合并的多尺度分割

经过SLIC分割后，图像中生成紧密的过分割超像素块，保留了图像中的纹理信息，但生成过程中只考虑到图像局部特征，对图像分割过于破碎，因此需要对相似区域进行合并消除过分割，只保留相似度差值较大的区域。本文方法通过不断重复提高阈值，使用层次区域合并的方法生成不同的分割尺度集，以超像素里所有像素点的平均向量值为特征，Lab空间和位置的加权和表示超像素间相异度，采用区域邻接连接图(region adjacency graph，RAG)记录超像素间的邻接关系和相似度，逐步合并相似区域。区域连接图的定义为R=(V，E)，其中V为像素区域，E为与之相邻的区域和相似度。生成RAG如图2所示。

(a) 超像素邻接关系(b) 初始区域邻接图(c) 合并后的邻接图图2 区域邻接示意图Fig.2 Schematic of region adjacency graph

多尺度分割主要通过“增大阈值”和“搜索合并小于阈值邻接区域”两部分完成，图2(a)是经过SLIC分割后的超像素的分布关系图，首先对邻接图进行初始化，可以由图2(b)表示，连线表示像素间的相交关系，连线的长度表示超像素间的相似度，然后设定较小的初始阈值，采用深度优先的搜索策略寻找小于合并阈值的相邻超像素对，找到满足条件的相邻像素对(S1，S2)后初始化合并后的区域S12属性信息和连接信息，在RAG中删除原始区域S1，S2的信息，将S12插入到RAG中，并更新所有与S12相邻区域的邻接关系，S1和S2合并后的邻接关系如图2(c)所示，然后继续搜索合并直至所有区域相似度大于阈值。随着阈值的增加，像素块根据相异性不断被合并，当图像由过分割进入欠分割状态时结束，利用二元分割树记录合并过程，可以快速回溯合并中的任意阈值，生成多尺度分割数据集，合并过程如图3所示，选择合适的阈值分割结果。

图3 尺度集层次结构图Fig.3 Illustration of scale-sets hierarchy

分割结果可以通过两个相对的指标来评估： 区域内部同质性和区域间异质性，最优分割需要平衡这两个指标。区域内部同质性表现的是一个区域内像素的一致性，通过光谱特征的局部方差(local variance，LV)来评估，整个分割结果的整体同质性可以通过各个区域加权计算：

(6)

式中：M和N分别为图像的宽度和高度；n为当前阈值的区域数量；si和σi分别为区域i的面积和光谱的标准差。分割结果的整体同质性通过区域的面积来衡量，合并后的较大区域比未合并的小区域权重大，光谱标准差σi的值越小内部同质性越好。

区域间异质性是相邻区域间的相互差异，可以通过空间自相关指标评价，这里使用莫兰指数(Moran’s index，MI)来计算衡量区域间异质性：

(7)

利用区域内部同质性和区域间异质性来估计最优尺度参数，随着阈值的增加区域合并范围变大，区域内同质性逐渐减小，区域间异质性增大，为选择图像由过分割到欠分割的邻接阈值，将LV和MI结合的目标函数GS的最小值作为最优尺度，公式为：

GS(k)=MInorm(k)+LVnorm(k)

(8)

MInorm和LVnorm分别是MI和LV的归一化后的值，即

(9)

(10)

1.3 A*算法路径搜索

A*搜索算法是启发式搜索算法中的一种，可以高效地求解在复杂结构中寻找路径的问题，在电脑游戏寻径计算、最短路径计算、机器人路径规划及最佳航线规划等领域具有广泛的应用[17-20]。A*算法结合了Dijkstra算法和最佳优先算法优势，减少了计算同时找到一条最佳的路径。由配准影像有效区域轮廓交点作为路径寻优的起始点，多尺度分割的区域内部作为障碍物，提取所带纹路信息的分割线作为通路，A*算法在图像分割线中搜索步骤如下：

1)将点A作为起点放入open列表中，然后检测和A相邻的点，如果是障碍物则忽略，如果是可走路径则放入open列表中，同时将A点从open列表移到close列表中，并设置A为父节点。

2)计算open列表中的点的移动属性F，根据属性值F从列表中选择新的扩展点，F可以由以下公式计算：

F=G+H

(11)

式中：G为由起点A点到指定点位移动代价；H为由指定点位不考虑障碍因素移动到终点B的估算代价，这里使用Manhattan距离，即当前点位不考虑斜角，只横向或纵向移动到终点经过的点位数。

3)选取F值为最小的点T作为扩展点并检测其相邻点位，其相邻点位是障碍物或已经在close列表中的直接忽略，作为通路的相邻点没有在open列表中则直接加入open列表中。如果已经在open列表中则按照当前路径重新计算G值并与原G值比较，如果当前路径G值大于原列表G值，则原列表中点位属性不变，如果当前路径G值小于原列表G值，更新列表中的点位属性，并将T作为父节点移入close列表中。

重复步骤2)直到终点B加入open列表中，完成搜索后由B回溯父节点构成镶嵌线。

2 实验及分析

2.1 实验数据

高分二号(GF-2)卫星于2014年8月19日在太原卫星发射中心用长征四号乙运载火箭发射成功，是目前我国分辨率最高的民用光学对地观测卫星，其幅宽相比于高分一号和资源三号卫星是最小的[21]，所以更需要对影像进行拼接。高分二号卫星参数见表1。

表1 高分二号卫星参数Tab.1 Gaofen-2 satellite parameters

本次实验研究影像使用3组正射校正影像，其中2组为天津市地区(图4(a)和(d)，图4(b)和(e))，1组为青岛市地区(图4(c)和(f))，可以看出3组数据中包含农田、建筑、河流、道路等常见典型地物，由于拍摄时间不同，每组影像间的色调有所不同。

(a) 区域1左下(b) 区域2左上(c) 区域3左上

(d) 区域1右上(e) 区域2右下(f) 区域3右下图4 正射校正影像数据Fig.4 Orthorectified image

2.2 实验结果

可以看到GF-2影像经过校正后影像四周出现无效区域，通过采样可以判断有效区域的范围，通过生成影像有效区域多边形并求交集，可以获取重叠区域影像，并以多边形的交点A和B作为镶嵌线的起始点(图5)。

图5 重叠影像获取Fig.5 Overlapping image acquisition

对重叠区域影像，SLIC生成超像素，经过不断地提高阈值和区域合并，得到不同阈值合并后的一系列影像，计算影像的LV和MI同时进行归一化，图6是经过归一化的LV和MI的变化趋势图，

图6 不同分割尺度下LV，MI和GS的变化图Fig.6 Variation diagram of LV, MI and GS under different segmentation scales

可以看出随着合并阈值的增加局部方差逐渐增大，MI逐渐减小，当GS最小时即图像由过分割状态进入欠分割状态，因此可以选择尺度为25的分割结果作为镶嵌线提取的基础。图7展现了在建筑、河流、农田等不同地物的分割效果。分水岭算法是一种基于拓扑理论的数学形态学的分割方法，也是一种优秀且得到广泛应用的分割技术，但可以从图6(b)看出分水岭算法对小面积的建筑区域轮廓提取时不够精细，同时根据图7(f)可以看出对于大区域的河流进行了过分割，不利于后续对镶嵌线的提取。图7中的(c)和(g)是SLIC预分割生成影像，可以看出超像素间非常紧密，内部特性单一，虽然对图像分割过于破碎，但对地物的轮廓非常敏感，很好地保留了地物轮廓信息。图7中的(d)和(h)是区域合并后的影像，从结果可以看出无论对于小面积的区域的农田，还是对于大区域的河流农田，基于区域合并的多尺度分割均保持良好的分割效果。对多尺度分割后的遥感影像进行镶嵌线的提取，设置分割区域内部区域为障碍物不能通过，因此镶嵌线只能在分割线上进行提取，镶嵌线提取结果如图8所示。其中红色的线是使用ENVI软件基于Voronoi图的镶嵌方案提取的镶嵌线，基本只是对于重叠影像的面积进行平分，没有考虑影像特征，接缝线明显； 黄色线是基于重叠影像差异的方法使用Dijkstra算法进行寻优提取得到的，提取的镶嵌线为影像间灰度差异最小，但由于在差值影像上建筑、河流等区域像素是不连续的，提取过程不能考虑图像地物信息，会导致穿越大片建筑、内部均匀的农田和河流等区域； 绿色为使用本文方法提取的镶嵌线，本文方法提取的镶嵌线尽量沿着图像中的道路、农田和河流等地物的边界，镶嵌线可以很好地隐藏在地物信息之中，通过纹路信息达到自然地过渡，避免穿过明显地物特征。

(a) 建筑区域原图 (b) 建筑区域分水岭分割 (c) 建筑区SLIC分割 (d) 建筑区SLIC合并

(e) 河流区域原图 (f) 河流区域分水岭分割 (g) 河流区域SLIC分割 (h) 河流区域SLIC合并图7 典型地物分割效果图Fig.7 Typical feature segmentation effect diagram

(a) 区域1(b) 区域2(c) 区域3

图8 镶嵌线提取结果对比Fig.8 Mosaic line extraction contrast

如表2所示，通过对比ENVI软件和Dijkstra模型穿越地物情况，ENVI软件穿越地物最多，提取的镶嵌线不能直接用于影像间的镶嵌； Dijkstra模型基于影像重叠区域最小差异能够一定程度降低镶嵌痕迹，但穿越明显地物数量仍然很多，在不能满足高质量镶嵌要求； 本文方法穿过障碍物数量最少，可以大幅度减少由于地物被分割而产生的错位现象，保留影像中真实地物特征，镶嵌线自动提取方法获得较好的影像拼接结果，能够满足大范围数字正射影像拼接生产的需求。

表2 各方法穿越明显地物数量Tab.2 The number of obvious objects traversed by each method (个)

3 结论

本文提出了一种基于多尺度分割的高分辨率遥感影像镶嵌线自动提取方法，使用天津地区和青岛地区3组GF-2数据进行测试，并与商业软件ENVI软件和基于重叠区域差异的Dijkstra算法对比，可以得出以下结论：

1)本文方法可以有效地处理高分辨率遥感影像镶嵌，生成的镶嵌线能够较好地避免穿过建筑，河流及农田等内部区域处理结果优于ENVI和基于重叠区域差异的Dijkstra算法。

2)对SLIC超像素区域合并能够更精准地分割地物，使用尺度集模型解决调整分割尺度参数问题。与现有镶嵌线提取技术相比，本文方法的优点在于无需辅助数据，充分利用图像地物信息，根据区域对象的差异而不是像元的差异进行镶嵌线的提取，能够避免镶嵌线穿过地物内部，提高镶嵌质量，可广泛应用于高分辨率遥感影像地图制作。