拓扑插值与光谱特征结合的复杂城区边界提取

于莉楠，宁晓刚，王 浩，刘纪平

中国测绘科学研究院，北京 100830

城区边界提取在城市空间布局、城区界限管理等方面发挥着重要作用[1-2]。传统的手工提取城区边界存在工作量大、花费时间长等问题，因此，研究城区边界的自动化提取具有重要意义。复杂城区边界常存在建筑物大小差异较大、周边林区高低起伏、散点建筑成片出现的情况，为自动化提取带来一定困难。高空间分辨率(简称“高分”)遥感影像使地理信息要素表现为多尺度的统一[3]，使目标的局部特征得到充分展现，因此，结合高分遥感影像进行复杂城区边界提取成为一条有效途径。

城区即城市区域，是发展较快、现代化程度较高、分布较集中的居民地，主要包括城市建成区和与之紧密相连的城乡结合区。单一城区一般具有集中连片的特征，即其内部可包含建筑物、绿地、小型河流、道路等，但不被大型河流、大面积绿地等非建设用地完全隔断。本文根据城区在遥感影像上呈现的复杂程度，是否有利于城区边界提取，将城区进一步分为“理想城区”和“复杂城区”。“理想城区”即建筑物大小较均一、城区周边较平整、城区边界较清晰、有利于城区边界提取的城区(如图1(a)所示)；“复杂城区”即建筑物大小差异较大、周边林区高低起伏、散点建筑成片出现、不利于城区边界提取的城区(如图1(b)、(c)、(d)所示)。本文的复杂城区边界提取即将复杂城区作为集中连片的整体，获取其外部边缘线。城区边界可较好地体现城区的宏观形态特征。城区居民地是城区的主体部分，常被作为确定城区边界的依据。城区居民地提取方法主要分为两大类。第1类是基于分类的提取方法[4-6]，又包括基于像元的分类和面向对象的分类。文献[7]利用基于像元的分类，结合形态学处理，实现了居民地提取。文献[8]采用面向对象的分类，利用北京大兴区的高分遥感影像提取了居民地。文献[9]对比了两种分类方法的提取结果，认为面向对象的分类方法更适用于高分遥感影像。文献[10]将建筑物指数(MBI)[11]特征与支持向量机(SVM)、随机森林(RF)分类器相结合，利用面向对象的分类来提取城区。基于像元的分类方法可充分利用影像的光谱特征，但当存在大量“同谱异物、同物异谱”的情况时，效果不佳；面向对象分类可综合利用光谱、纹理、形状等特征，通过图像分割弱化“同谱异物、同物异谱”现象，但对图像分割质量要求较高，分割尺度不易确定。

第2类是基于特征的提取方法，又分为基于局部特征(如角点[12-16]、直线段[17]、边缘[18-19])和基于全局特征的(如光谱、纹理)方法。局部特征一般具有较好的局部不变性，如文献[11]结合了居民地的Harris角点；文献[17]基于直线段统计特征提取居民地；文献[18—19]充分利用了居民地的边缘密度特征；文献[20]则同时融合了直角点和直角边进行居民地提取。全局特征则是从整体上表达同类地物的共同特征。如文献[21]通过分析各地物光谱曲线，构建了居民地提取专家模型；文献[22]基于模糊准则下的灰度共生矩阵来提取居民地；文献[23]则提出一种将纹理特征与视觉注意相结合的方法。基于特征的方法在特征显著时效果较好，但当存在大量具有相同特征的非目标干扰时，或提取的特征点、线不充足的情况下，无法达到理想效果。

图1 城区复杂状况Fig.1 Complex situation of urban boundary

现有方法对于建筑物大小较均一、城区周边较平整、城区边界较清晰的“理想城区”，基本都能取得较好效果。而对于建筑物大小差异较大、周边林区高低起伏、散点建筑成片出现的“复杂城区”，会因为光谱、纹理等特征过于复杂而效果不佳。针对这一问题，本文提出一种基于拓扑插值和光谱特征的高分遥感影像复杂城区边界提取方法。一方面，分别提取大型建筑和小型建筑的特征点，并通过由点及面的拓扑插值，来避免城区边界大建筑特征点稀疏带来的城区边界缺失；另一方面，对于城区周边高低起伏的林区、散点建筑带来的冗余特征点，基于植被光谱特征将其滤除，从而抑制城区边界误提取。该方法通过融合点、线、面及光谱特征，实现了复杂城区边界提取。

1 本文方法

本文基于拓扑插值和光谱特征的城区边界提取方法，主要包括以下5步：①大、小建筑物特征点的分别提取；②基于拓扑插值的建筑物特征点加密；③基于光谱特征的冗余特征点滤除；④城区特征图像构建；⑤城区边界矢量生成。本文方法流程如图2所示。

图2 本文方法流程Fig.2 The work flow of the proposed method

1.1 大、小建筑物特征点的分别提取

城区在高分遥感影像上常呈现为具有直角特征的建筑顶面，本文即以直角点作为城区特征点，并根据建筑物边长阈值，将建筑物分为大建筑和小建筑。为了对大、小建筑特征点进行不同程度的加密，首先需要分别提取大、小建筑物的直角点。

文献[20]通过Harris角点结合直角边约束来获取直角点。该方法对于理想城区效果较好，但对于复杂城区，由于Harris角点检测对纹理特征非常敏感，会在复杂城区周边高低起伏的林区和散点建筑处产生大量冗余角点，且这些角点对应的直角边长度与城区小型建筑边长相近，很难基于边长约束将其滤除。这种先提取点、再检边、再筛点的流程对于复杂城区不仅耗时，且效果不佳。因此，本文通过垂直边求交直接获取直角点来简化这一流程。直线段检测的常规方法有Hough变换[24]、LSD直线段提取算法[25]，但Hough变换较难提取细小线段，LSD算法容易使线段在线与线相交处断开，不利于长线段提取，且参数较少，不易控制线段提取的具体情况。

本文采用Canny边缘检测[26]、Douglas-Peucker直线段压缩[27]、角度约束相结合的方式来检测垂直线段。其优势在于可根据城区复杂程度，通过Canny边缘检测的参数调整来控制提取线段的数量，通过Douglas-Peucker的线段简化程度来控制线段提取的长度。具体步骤如下：

(1)基于Canny边缘检测进行边缘提取。

(2)基于Douglas-Peucker算法进行直线段压缩，即通过设定不同的距离阈值，对建筑物边缘实现两种不同程度的简化(如图3所示)。当简化程度较大时，较小的边缘起伏会被忽略，从而得到较长的线段集合，记为“粗略线段集合”LA；当简化程度较小时，即使较小的边缘起伏也会将边缘分割成线段，从而得到较短的线段集合，记为“精细线段集合”LB。

图3 基于Douglas-Peucker算法的直线段提取Fig.3 Line extraction based on Douglas-Peucker algorithm

(3)对于“粗略线段集合”LA，操作如下：①设线段di是LA中的第i条线段，其端点分别为pi1、pi2(图4)，先以其中一个端点为基准，在半径为r的圆域内，搜索与该端点最邻近且夹角近似等于90°的线段dj，并设其端点分别为pj1、pj2；②当线段di的长度li、线段dj的长度lj与长度阈值lthreshold满足式(1)，且二者的夹角θij与角度阈值θthreshold满足式(2)时，近似认为这两条线段垂直且来源于大建筑边缘，求解两条线段所在直线的交点cij，并将其作为直角点；③分别将线段di、dj的端点中与cij距离较远的端点记为pi、pj，将点集Cij={cij,pi,pj}作为子点集放入点集合UA中，对线段另一个端点做相同操作；④依次遍历完LA中的所有线段，即得到大型建筑初始特征点集UA。

图4 直角点获取Fig.4 Right angle point extraction building feature points

(1)

|θij-90|<θthreshold

(2)

(3)

(4)对于“精细线段集合”LB，其操作过程中的①③④与步骤(3)中的①③④相同，但在步骤②中，li与lj需满足式(3)，最后将得到的点集合记为小型建筑初始特征点集UB。

1.2 基于拓扑插值的建筑物特征点加密

为避免大型建筑稀疏特征点带来的城区边界缺失，本文通过先拓扑、后插值来增加大型建筑特征点密度。对于小型建筑，不进行拓扑，直接简单插值。通过这种由点及面、大小建筑区别处理的方式，使特征点得到不同程度的加密，从而更好地体现城区整体特征。具体过程如下：

(1)设Cij={cij,pi,pj}为大型建筑初始特征点集UA中的一个子点集，即由相互垂直的线段di、dj得到的特征点集，cij、pi、pj的图像坐标分别为(xc,yc)、(xa,ya)、(xb,yb)。由子点集获取过程可知，以cij、pi、pj为顶点可构建直角三角形，以直角边为邻边，以斜边为对角线，将其拓扑为矩形(图5)。只要插入第4点t1即可构建一个矩形面，t1(x1,y1)的坐标，计算如式(4)

(4)

(2)在以Cij={cij,pi,pj,t1}中的4个点为顶点构建的矩形上，并在每条边的四等分点处插值，如图5所示，得到t2、t3、t4、t5、t6、t7、t8、t9、t10、t11、t12、t13共12个边缘插值点。

图5 大型建筑特征点拓扑插值Fig.5 Topological interpolation of big feature points

(3)在以Cij={cij,pi,pj,t1}中的4个点为顶点的矩形内部进行插值，在对角线pipj的中点位置插值得到矩形的中心点q1，在q1与cij、pi、t1、pj连线的中点上分别插值得到q2、q3、q4、q5，即得到5个内部插值点。

(5)

图6 小型建筑特征点插值Fig.6 Interpolation of small building

1.3 基于光谱特征的冗余特征点滤除

由于城区周边高低起伏的林区和散点建筑附近常存在绿色植被、林地或田地，使这些非城区部分产生的冗余特征点具有了植被光谱特征。因此，本文采用归一化植被指数(NDVI)，并定义绿色程度指数(GI)来滤除这些冗余特征点，其计算方法分别如式(6)、式(7)所示。

(6)

(7)

NDVI适用于含有红(R)、绿(G)、蓝(B)、近红外(NIR)波段的高分遥感影像，而对于不含NIR的RGB遥感影像，可结合GI进行滤除判定。结合图7，冗余特征点滤除步骤说明如下：

(1)设p(x0,y0)为高分遥感影像上的特征点，若该影像含有R、G、B、NIR 4个波段，设植被指数阈值为Vindex_threshold。当NDVI(x0,y0)>Vindex_threshold时，直接将其从特征点集中去除；当NDVI(x0,y0)≤Vindex_threshold时,在以(x0,y0)为中心，大小为N×N的邻域边界上(如图7中的斜线区)，累计满足NDVI(x,y)>Vindex_threshold的像素个数n,当满足式(8)时(Rration为比例阈值)，将其从特征点集中去除

(8)

图7 冗余特征点滤除窗口Fig.7 Process window for removing redundant feature points

(2)若p(x0,y0)所属的高分遥感影像只含有R、G、B3个波段，采用绿色程度指数GI来滤除冗余特征点，方法与基于NDVI的特征点滤除方法相同。

1.4 城区特征图像构建

在高分遥感影像上，像素点距离城区特征点越近，其属于城区的概率越大，而二维高斯函数可较好地描述这一概率分布特征。本文结合二维高斯函数，采用不同标准偏差，通过局部空间投票来生成与原遥感影具有相同大小和分辨率的城区特征图像。生成过程如下：

(1)设p(x,y)为遥感影像上某像素点，在以p(x,y)为中心，半径为R的矩形局部窗口内，搜索城区特征点，并将搜索到的边缘特征点(包括直角点、边缘插值点、边缘线段端点)个数设为m，将搜索到的内部插值特征点(即不位于边缘线段的插值特征点)个数设为n。

(2)通过局部空间投票，计算p(x,y)属于居民区的概率PR(x,y)，并将其作为特征图像在坐标(x,y)位置的特征值，PR(x,y)的计算方法如式(9)所示。其中xi、yi(i=1,2,3，…，m)为边缘特征点坐标，xj、yj(j=1,2,3，…，n)为内部插值特征点的坐标，σ1、σ2分别为边缘特征点和内部插值特征点对应的高斯函数标准偏差，lthreshold为用于划分大小建筑物的长度阈值。

图8为特征点加密过程及效果，图9为冗余特征点滤除过程及效果。

图8 基于拓扑插值的建筑物特征点加密Fig.8 Topological interpolation for getting density building feature points

图9 基于光谱特征的冗余特征点滤除Fig.9 Removing redundant feature points according to spectral feature

(9)

1.5 城区边界矢量生成

1.5.1 特征图像二值化

复杂城区特征图的特征强度变化较多，常规自动阈值确定方法较难达到理想二值化效果。利用笔者编写的交互界面(图10(a))，基于特征图的像素值，采用半自动方法确定二值化最佳阈值，即通过调整交互界面上的阈值滑块，将该阈值下的二值化结果对应的白色区域(即城区图斑)与遥感影像进行实时叠加显示，并参考遥感影像调整阈值。当白色区域刚好覆盖遥感影像上的城区部分时，即找到最佳二值化阈值，进而实现特征图像二值化。

1.5.2 城区图斑修整

城区内非建筑区(如水体、道路等)因缺少建筑特征点，易使特征图二值化结果出现空洞，而城区周边未被完全滤除的冗余特征点则易产生冗余小图斑。因此，本文基于连通域分析法，对图斑进行修整。首先，将被白色区域包围的黑色空洞填充为白色，即在交互界面上用鼠标框选欲填充的最大黑色空洞(图10(b))，所有面积小于该选框面积的黑色空洞均会被填充，从而将城区图斑补充完整；然后，框选欲去除的最大白色小图斑(图10(c))，所有面积小于该选框的白色小图斑将被去除，从而保留城区主体(如图10(d))，完成城区图斑修整。

图10 城区图斑生成Fig.10 Urban map spot generation

1.5.3 城区图斑矢量化

基于栅格转矢量功能实现二值图像矢量化，矢量结果中属性值为“0”的矢量图斑代表非城区，将其去除后即可得到城区对应的矢量边界。

2 试验结果与对比分析

为验证本文算法的有效性，基于Visual Studio 2010 C++开发环境实现了该算法，然后结合复杂城区高分遥感影像进行了试验，并将本文方法结果与文献[20]、文献[22]方法结果进行了对比分析，最后基于手工勾绘结果进行了定量精度评价。

2.1 数据来源及结果

本文选取了两幅城区边界较复杂的高分遥感影像作为试验数据，具体情况如表1所示。数据1(图11(a))的城区边界分布着大小差异较大的建筑，用该数据验证基于拓扑插值的特征点加密算法的有效性(图12)。数据2(图13(a))的城区周边分布着大片高低起伏的林区和散点建筑，用该数据验证基于光谱特征的冗余特征点滤除效(图14)。

图11给出了数据1试验过程中生成的特征图，图中颜色越偏于红色，表明特征强度越大，其为城区的概率越高。由本文方法特征图(图11(b))可看出，经拓扑插值后，大型建筑区与小型建筑区均体现较大的特征强度。文献[20]特征图(图11(c))中，大部分建筑区特征强度明显，但两个大建筑和部分小建筑密集区特征强度较弱。文献[22]特征图(图11(d))中，由于大型建筑区纹理特征不明显，其特征强度明显小于小型建筑区。

图11 数据1及生成的特征图Fig.11 Test data 1 and feature image generated

图12 数据1试验结果Fig.12 Test results of the test data 1

图13 数据2及生成的特征图Fig.13 Test data 2 and feature image generated

图12(b)、(c)、(d)分别为本文算法和文献[20，22]方法基于数据1提取的城区边界矢量。通过对比看出，本文算法可使大型建筑区对应的城区边界得到较好的保留，文献[20]方法在个别大建筑和部分小建筑处出现了部分边界缺失，而文献[22]方法则在大型建筑区对应的城区边界处出现了明显缺失。参照图12(a)所示的手工勾绘结果来看，本文算法可有效避免建筑物大小差异较大带来的城区边界缺失。

表1 试验数据Tab.1 Test data

图13为几种方法针对数据2生成的特征图，图14为数据2的试验结果。图14(b)为本文算法针对数据2的试验结果，参照特征图(图13(b))可看出，城区周边高低起伏的林区和散点建筑带来的冗余特征点被较好地滤除了，使非城区部分的特征强度明显低于城区，从而呈现出较明显的城区边界；图14(c)为文献[20]的试验结果，其中大部分建筑区提取较准确，但为了基于长度阈值滤除高低起伏的林区对应的Harris角点，该长度阈值与城区内部小建筑边长相近，因而损失了部分城区内特征点，且道路特征点稀疏，导致城区提取不够完整。结合图13(d)来看数据2对应的文献[22]处理结果(图14(d))，由于城区部分和周边林区均具有较强的纹理特征，从而呈现出相近的特征强度，因此无法确定边界。结合手工勾绘的结果(图14(a))可看出，本文算法有效抑制了冗余特征点造成的城区边界误提取。

2.2 精度评价

为定量评估本文算法所提边界的准确性，以手工勾绘的城区边界为基准，采用正确率Pc、完整率Pe、质量Pq来对结果进行精度评价。这3个评价指标的计算公式如式(10)所示。其中Sauto为算法提取城区矢量面积，Smanual为手工勾绘城区矢量面积，Sauto&manual为算法提取结果与手工提取结果交集的面积，Sauto‖manual为算法提取结果与手工勾绘结果并集的面积。表2给出了本文方法和文献[20，22]方法的精度评价结果。通过对比可以看出，本文方法相对于文献[20，22]方法具有更高的正确率、完整率和质量，更接近于手工勾绘结果。

(10)

表2 城区边界提取精度评价Tab.2 Accuracy statics of urban boundary extraction (%)

2.3 分析与讨论

通过目视对比和定量精度评价可看出，本文方法总体上优于文献[20，22]方法。文献[20]方法通过融合直角点、直角边特征，可以较好地提取大部分城区，但仍然无法避免复杂城区中大型建筑造成的边界缺失,在滤除高低起伏的林区对应的Harris角点时，损失的部分城区特征点也造成了城区的不完整。文献[22]方法借助单一纹理特征无法应对纹理特征缺失或纹理特征过于繁杂的状况。本文算法通过拓扑插值有效避免了大型建筑稀疏特征点造成的不良影响，结合光谱特征有效抑制了林区冗余特征点造成的干扰，在一定程度上提高了城区边界的准确性。

本文方法在参数设置方面需结合遥感影像实际情况和经验知识。Canny边缘检测阈值应使城区内大部分建筑边缘被准确提取，且使城区建筑边缘数量与所提取总边缘数的比值尽可能大。精细线段集的Douglas距离阈值尽可能保证小建筑直角不被简化掉，粗略线段集Douglas距离阈值应使大建筑边缘上的锯齿起伏尽可能少。边长阈值lthreshold建议取大建筑平均边长。目前本文方法主要适用于分辨率不低于2 m的高分遥感影像。

本文方法虽然主要根据城区内建筑物特征来确定城区范围，但其内部水体、道路、绿地也并未被完全遗漏，因为这些缺少建筑特征的区域产生的空洞在进行“图斑修整”时会被填充，从而保证了整个城区的完整性，但当这些区域与城区边界相接时，会造成一定的边界误差。在进行拓扑插值时，插值点有可能被插值到城区边界的外围(如数据1中城区周边纵横交错的耕地)。当这些“误插值”点落于植被区，会基于植被光谱特征被滤除；当落于非植被区，由于这些“误插值”特征点产生的特征强度相对于主体城区来说一般较弱，且相对零散，在进行特征图像二值化时会因特征强度达不到分割阈值而被滤除；部分“误插值”特征点产生的零散图斑会在“图斑修整”时，基于面积阈值被滤除，从而在一定程度上削弱了“误插值”带来的不良影响。在滤除城区周边高低起伏的林区造成的冗余特征点时，城区内的凹凸植被产生的特征点也会被滤除，这种“误滤除”产生的部分图斑空洞也会在“图斑修整”时被填充，使城区内的植被尽可能不被遗漏。但当这些城区内植被与城区边界相接时，也会在一定程度上影响城区边界的准确性。本文试验结果中的误提取主要存在于城区周边的凹凸裸地，针对该种复杂状况的城区提取后续可进行进一步研究。

3 结 论

本文提出了一种基于拓扑插值和光谱特征的高分辨率遥感影像复杂城区边界提取方法。该方法充分融合了点、线、面及光谱特征，不仅避免了大型建筑特征点稀疏造成的城区边界缺失，而且有效抑制了高低起伏的林区和散点建筑带来的城区边界误提取。试验结果表明，该方法整体上优于文献[20，22]方法，可较好地适用于复杂城区边界提取，且具有较高的精度。本文创新点主要体现在三方面：①针对复杂情况下的城区边界提取提出了解决方案；②将拓扑插值理论应用于城区边界提取中；③同时融合了点、线、面及光谱特征。该方法的不足之处在于：①参数及阈值的选择需要一定经验，无法达到全自动化；②计算效率还不够高。今后将着力解决这些不足之处，通过算法优化，简化参数设置，降低算法复杂程度，提高算法自动化程度，并通过GPU并行加速等方式进一步提升效率，从而达到更好的效果。