多星形约束图割与轮廓规则化的高分遥感影像直角建筑物提取

丁亚洲，冯发杰，吏军平，胡 艳，崔卫红

1. 湖北省电力勘测设计院有限公司，湖北 武汉 430000； 2. 西安科技大学测绘学院，陕西 西安 710000； 3. 武汉大学遥感信息工程学院，湖北 武汉 430000

建筑物作为地理空间主要要素之一，在城市规划与建设、变化检测及人口密度估计等领域占据重要地位。高分遥感影像能更清楚地反映面状地物的几何、纹理等空间特征，可相对精确地描述地物，这使得利用高分遥感影像对建筑物进行精确提取成为可能[1-3]。目前，根据人工参与程度可将高分遥感影像建筑物提取方法分为自动和半自动提取。自动提取方法多基于建筑物的边缘、角点及阴影特征。文献[4]提出了一种结合边缘保护、双边滤波和直线检测提取建筑物的方法，通过检测影像上的线段，并进行连接得到建筑物轮廓。文献[5]基于建筑物的阴影特征，先检测建筑物和植被阴影，保留建筑物阴影，再根据阴影位置确定建筑物位置，最后结合超像素分割、合并及线段检测得到建筑物轮廓。文献[6]提出了一种结合图像分割和角点检测的建筑提取方法，先利用Graph Cuts方法获取建筑物图斑，再检测建筑物拐角点，连接拐角点获取建筑物。然而，这些方法仅适用于提取形状简单且轮廓清晰的建筑物。文献[7]提出了一种利用差分形态学剖面构建形态学建筑物指数自动提取高分辨率遥感影像建筑物的方法，且进一步优化了该方法[8-9]。文献[10]基于深度学习方法自动提取影像上的建筑物，但该方法在建筑物准确轮廓和复杂建筑物提取方面还有待进一步提高。

由于遥感影像包含的地物复杂多样、建筑物类型较多，以及建筑物受噪声、遮挡、阴影、低对比度的影响，自动提取建筑物的方法难以实现较高的提取精度。实际上一些学者提出的半自动提取建筑物方法取得了更好的效果。文献[11]提出了一种交互式半自动提取高分辨率遥感影像上建筑物的方法，先指定建筑物中心，再采用径向投射算法获取Snake初始化轮廓，最后得到建筑物真实轮廓。文献[12]结合Snake模型和动态规划方法，在房角处指定几个种子点，即可提取建筑物轮廓。文献[13]提出一种结合影像分割和区域合并的半自动提取建筑物的方法,先利用Mean Shift方法将包含建筑物的影像图块分割成一些小区域，再基于人工交互合并相似区域得到建筑物轮廓。文献[14]结合几何约束和影像分割实现直角平顶房屋的半自动提取。然而上述方法针对复杂建筑物，不仅交互复杂且不能保证提取精度，因此提取复杂形状和纹理的建筑物仍然是目前亟待解决的难题。随着图割优化方法的发展，许多学者将形状先验信息作为约束项融入图割优化框架，以提高分割结果的准确性[15-16]。文献[17]将椭圆形状约束加入图割模型中，文献[18]将凸形状约束加入图割模型中，文献[19]提出一种单星形先验的图割方法。文献[20]将单星形扩展为多星形，并引入Graph Cuts模型中。已有一些学者将形状约束加入图割中用于道路和建筑物的提取。文献[21]提出了一种基于先验形状约束水平集模型的建筑物提取方法，将多个先验形状竞争模型引入水平集方法中,利用先验形状来约束曲线的演化,在对图像进行分割的同时完成建筑物的检测和提取。文献[22]提出了一种结合形状先验的图割和动态外推思路自动提取高分遥感影像上道路的方法。

基于上述分析，由于高分遥感影像上绝大多数建筑物具有直角结构和星形形状特征，因此本文在人工交互获得目标前景与背景的基础上，结合多星形约束的Graph Cuts模型和轮廓规则化方法，提取高分遥感影像上的直角建筑物，以实现在建筑物上简单交互即可快速准确地提取建筑物轮廓。

1 研究方法

本文基于分割获取建筑物图斑以及与建筑物方向一致的边缘直线，通过连接直线交点得到建筑物准确轮廓。该方法主要分为人工交互、影像预处理、建筑物图斑获取和图斑规则化4步，具体流程如图1所示。首先，通过人工交互，在原始影像的目标建筑物上画种子线确定目标建筑物的大致范围，获取包含目标建筑物的影像图块；其次，通过双边滤波对影像块进行预处理，并探测建筑物主方向；再次，用SLIC超像素分割方法对预处理后的影像块进行超像素分割，根据种子线得到目标候选前景与背景超像素，用多星形约束的图割方法获取建筑物图斑；最后，利用Harris算子检测建筑物图斑上的角点，并对角点进行分组拟合，将拟合的直线旋转至建筑物主方向，连接直线交点得到规则的直角建筑物。

图1 建筑物提取过程Fig.1 The framework of building extraction

1.1 人工交互

本文方法需人工在建筑物上指定种子线得到目标建筑物大致位置与范围，对种子线有两点要求：①种子线尽量处于建筑物的中间位置；②种子线外包矩形长度应当长于建筑物的2/3。如图2所示，绿色线为人工交互的种子线，Wseed和Hseed分别为种子线外包矩形S1的宽和高，O为S1的中心。建筑物的大致范围S2是一个以O为中心，边长为2×max(Wseed,Hseed)的正方形。根据目标建筑物的大概范围，选取以O为中心，边长为2.5×max(Wseed,Hseed)的正方形区域S3作为包含目标建筑物的影像图块。

图2 建筑物大致位置与范围Fig.2 The approximate location of the building

1.2 影像预处理

由于高分遥感影像上的大量噪声降低了影像质量，因此对影像预处理是必不可少的。双边滤波[23]是一种简单有效的非线性滤波方法，与传统的高斯滤波器相比，其考虑了图像的空域信息和灰度相似性，因此能够达到保边去噪的目的。以BF代表双边滤波的运算，其定义如式(1)所示

(1)

式中，Wp是一个标准量，如式(2)所示

(2)

式中，σd和σr分别控制空间域和亮度域特征；

Gδd是一个空间函数，用于减少远距离像素影响；Gδr是一个范围函数，用于减少灰度值不同于Ip的像素q的影响。本文利用双边滤波对影像图块进行预处理，σd和σr分别取10和30，以图3(a)为例，滤波结果如图3(b)所示。

图3 双边滤波器预处理结果Fig.3 Result of bilateral filter prepocessing

1.3 获取建筑物图斑

1.3.1 SLIC超像素分割

利用超像素取代单个像素将有效降低图像处理的复杂度。简单线性迭代聚类算法(SLIC)[24]通过构造一个由CIELAB颜色空间的[lab]值和像素坐标值[x,y]组成的五维特征参数，采用融合颜色相似度和像素空间位置邻近程度的归一化距离度量来对像素进行简单线性局部聚类，从而得到一系列具有相似特征且不破坏影像边界信息的均匀超像素。其算法简单、耗时少且需用户调节参数少。本文利用SLIC方法对预处理后的影像进行过分割得到超像素，再以人工交互所得种子线来确定目标候选前景与背景超像素。为使得过分割结果尽量准确，每个超像素大小设为10，如图4(a)为均匀超像素，图4(b)中的绿色线为种子线，将种子线所在的超像素作为前景，以1.1节中所得建筑物范围矩形S2与影像块边界S3之间的超像素为背景，如图4(c)中的红色和蓝色部分所示。

图4 目标前景与背景获取Fig.4 Obtain the foreground and background of the building object

1.3.2 星形形状约束下建筑物图斑获取

星形形状[19]即轮廓点对中心点可见，如图5(a)所示，闭合曲线y为一个单点星形形状例子。其中c是y的星形中心，p为y内的任意一点，c与p连线Γc,p上的任意点q仍在y内。即若用1和0分别标记目标前景与背景，那么根据定义可知，如果Lp=1，有Lq=1。则定义在c上的星形形状约束的能量函数表示如式(3)所示

(3)

图5 星形形状约束Fig.5 The star shape constraint

由于单星形形状约束只能处理凸形状的目标，文献[20]对该模型进行了改进，将一个星形中心改为具有多个中心的集合O，并将单条直线Γc,p扩展为O与p连线的集合ΓO,p。多点星形形状即该轮廓上的点至少被一个中心点可见，如图5(b)所示，y为一个多星形形状，其中c1、c2分别为星形形状中心，p点对于c1可见但对于c2不可见，因此该多边形形状属于以{c1,c2}为中心的星形形状。可见多星形适于形状复杂图形分割，多星形约束的能量函数如式(4)所示

(4)

图割即对像素点进行标号的过程[25]，其核心是设计能量最小函数，并基于最大流最小割优化理论，运用图映射和网络流知识，通过加权图的最小割来实现全局优化下的前景提取，能量函数为式(5)[26]

(5)

式中，L={Lp|p∈P}表示图像P的标记；Lp∈{0，1}；N为图像中所有邻接像素对组成的边集合；Dp(Lp)为数据项，表示像素p取标记Lp的代价；Vp,q(Lp,Lq)为边界项，表示相邻像素p、q分别取Lp、Lq的代价。

对于绝大多数建筑物而言，通过简单交互，很容易保证建筑物是关于种子线的星形形状，即使是对于特别复杂的建筑物也可通过添加多个中心保证建筑物具有多星形形状。因此本文采用基于多星形约束的图割模型，能量函数表达为式(6)

λE*(y|O)

(6)

式(6)是在式(5)的基础上增加多星形约束项E*(y|O)得到的，其中，γ、λ为平衡系数。

1.4 建筑物主方向探测

EDLines[27]是一种自动、快速且无需任何参数调整，能够准确提取影像上直线的直线检测算法。本文利用EDLines检测到的直线线段方向直方图确定建筑物的主方向。由于建筑物有复杂背景，导致很多无关直线对建筑物主方向的确定造成干扰，因此需要进一步确定建筑物的准确位置。本文在获取建筑物图斑后，利用图斑的外接矩形获取建筑物更准确的位置，如图6(a)上的蓝色矩形，(b)上蓝色矩形以内的图块作为建筑物的准确位置。然后利用EDLines提取(b)中矩形内部影像上的线段L=(l1,l2,…,ln)，如图6(c)上的红色线段；再将种子线中点作为建筑物的中心，计算影像上所有线段中点到建筑物中心的距离，根据提取的线段距离建筑物中心越近则位于建筑物上的概率越大的原则，利用距离值对线段赋权值，如图6(d)为影像上所有像素点到建筑物中心的归一化距离，称其为种子线距离图。线段方向直方图取值方式和建筑物主方向计算方式分别如式(7)和式(8)所示

(7)

θM=arg max(nθ+nθ+90)θM∈[0,90]

(8)

式中，θ∈[0,180)表示线段方向，同时也代表直方图横轴；nθ是直方图方向值为θ处的强度；li表示一条线段；Len(li)表示线段li的长度；D(xli,yli)表示线段li的中点到建筑物中心的归一化距离值；θM为建筑物主方向。通过上述计算，分别得到图6(e)所示的线段方向直方图与图6(f)所示的建筑物主方向。

图6 建筑物主方向探测Fig.6 Detection of main direction of building

1.5 图斑规则化

由于仅通过图割得到的建筑物并不是规则的建筑物，其边缘线条不够平滑，不满足建筑物为规则形状的特点，因此有必要对获取的建筑物图斑进行规则化处理，以得到准确的建筑物轮廓。本文利用Harris角点检测算法[28]提取建筑物图斑上的角点，然后依次经过角点分组、直线拟合与直线旋转相交3个步骤得到建筑物轮廓。首先进行角点检测(图7(a)中的绿色点)，采用每个角点与其距离最近的初始轮廓点的次序，对角点进行排序；然后利用角度阈值对检测到的角点分组，再对每组角点进行直线拟合(图7(b)中的绿色线)，最后将拟合直线旋转至与建筑物主方向一致或垂直方向(图7(c))，连接相邻直线交点(图7(c)中的红色点)得到建筑物轮廓(图7(d))。具体步骤如下：

(1) 对建筑物图斑轮廓上的点Ni按照顺时针进行编号，每个点的编号为i(i=0,1,2,…,N-1)；

(2) 计算每个角点到建筑物图斑轮廓上点Ni的欧氏距离di，若dj最小,则角点标号为j；

(3) 按照角点标号顺序对所有的角点进行排列得到有序的角点；

(4) 利用余弦定理计算相邻3个角点连线的夹角，根据角度阈值得到建筑物的近似拐角点，将相邻拐角点之间的点分为一组，完成对角点分组；

图7 图斑规则化过程Fig.7 The procedure of regulating building

(5) 采用最小二乘方法对每组角点拟合得到直线，将直线旋转至建筑物主方向或垂直于建筑物主方向；

(6) 依次连接相邻直线的交点，获得建筑物外轮廓。

2 试验结果与分析

为比较本文方法的有效稳健性，分别选取两个不同地区不同分辨率的高分辨率航空影像来进行试验分析，原始影像如图8所示，其中，数据1的空间分辨率为0.2 m，该地区位于广东省的惠州市，区域内房屋相对独立，且大部分为直角建筑物(图8(a))。数据2来源于国际摄影测量与遥感协会(ISPRS)，其空间分辨率为9 cm(图8(b))。两幅影像分别选取了82和57栋直角建筑物来定量评价本文方法，通过人工手动勾绘获取建筑物轮廓真值。

图8 试验原始影像Fig.8 Original imageries

选取常用的两种交互式算法GrabCut[29]及OneCut[30]与本文方法作定性和定量评价对比，为获取更好的结果，GrabCut算法迭代次数设为5次，OneCut算法中的直方图设置为64级。本文采用常用的评价指标，准确率(Precision)、召回率(Recall)和F1值(F1-measure)来定量评价3种方法的提取结果，3个指标的定义为

(9)

(10)

(11)

式中，Sauto为算法提取结果的总像素；Smanual为人工勾选的建筑物真值像素总数；Sauto∩manual为Sauto与Smanual匹配的总像素。

对两种数据上的建筑物提取结果分别如图9、图10所示。其中图9(a)、图10(a)中绿色标记为建筑物真值，(b)、(c)和(d)分别为GrabCut、OneCut及本文方法提取结果，两幅影像上提取差异较为明显的建筑物分别用蓝色和黄色圈标记。由图9的试验结果可知对于矩形形状、与背景差异明显的建筑物，3种方法提取结果均较好，但对于复杂形状，与周围背景差异小的建筑物，本文方法明显优于其他两种方法。由于数据2的分辨率较高，绝大多数建筑物为矩形形状且其轮廓较为清晰，因此图10中GrabCut与OneCut方法提取结果也较好。

为使得本文方法与其他方法具有更加明显的对比性，图11列出了3种方法对数据1中部分建筑物提取结果的对比示意图。由于OneCut分割模型应用直方图取代GrabCut的高斯混合模型(GMM)计算前景背景，同时在图割模型中引入了一些辅助节点使得分割结果更加准确，OneCut提取结果优于GrabCut。从试验结果上可以看出，多星形形状约束分割与本文规则化方法在建筑物提取上的有效性。

图9 数据1的3种方法提取结果对比Fig.9 The extraction results comparison of the three methods on the first aerial imagery

表1为两幅影像上所有建筑物提取结果的平均F1值和平均处理时间，其中本文方法对两幅影像上建筑物提取精度分别为92.4%和92.7%，相对于GrabCut和OneCut方法提取精度更高，这表明将形状约束集成在图割中能够提高分割精度。在多星形形状约束分割前对影像进行过分割，大大减少了计算量。本文方法提取每个建筑物的平均时间在0.6s左右，对于交互式分割而言这样的耗时是可以接受的。从表1可以看出，GrabCut提取建筑物需要的时间最长，这是因为GrabCut需要迭代更新。而本文算法由于增加了形状约束，使得速度比OneCut慢。图12和图13分别显示了3种方法提取每个建筑物结果的F1值，显然，本文方法提取结果的F1值波动较平缓，而GrabCut和OneCut方法波动幅度较大，这表明本文方法更稳定。

表1 不同方法提取两种数据建筑物的精度对比

图10 数据2的3种方法提取结果对比Fig.10 The extraction results comparison of the three methods on the second aerial imagery

图11 3种方法提取结果对比(数据1)Fig.11 The comparison of the three methods for some buildings on the first aerial imagery

图12 数据1上建筑物提取结果的F1值Fig.12 The F1-measure of building extraction result of the first aerial imagery

图13 数据2上建筑物提取结果的F1值Fig.13 The F1-measure of building extraction result of the second aerial imagery

3 结 论

本文在人工交互画种子线的基础上，结合多星形约束的图割和轮廓规则化方法，实现了对高分辨率遥感影像上直角建筑物的交互式半自动提取。分别对两幅不同地区和空间分辨率的遥感影像上的建筑物进行了提取试验。结果表明，多星形形状约束的图割适用于建筑物提取，且能够有效地提高分割精度，从而提高规则化精度。但本文方法基于图割，且受建筑物背景干扰影响较大，若建筑物与周围差异较小，轮廓不清晰，则导致建筑物图斑提取不准确。为了获取精确的建筑物轮廓，建议采用0.5 m以上分辨率的影像。同时，提出的规则化方法仅针对各种形状的直角建筑物，故对于非直角建筑物，本文方法无法对其规则化。如何实现各种非直角建筑物的精确提取及规则化将是下一步开展研究的内容。