CN103699900A - 卫星影像中建筑物水平矢量轮廓自动批量提取方法 - Google Patents

Abstract

本发明提供了一种卫星影像中建筑物水平矢量轮廓自动批量提取方法，首先利用K-均值算法进行图像分类，得到建筑物的主干部分，取各建筑物区域质心解决初始种子点的选取问题。生长出所有种子点所在区域后，借助图像边缘信息将建筑物区域与周围区域分开，并根据矩形度、带状指数等特征参数去除掉非建筑物区域，从而实现建筑物水平像素轮廓的自动提取。然后利用Hough变换、块处理等技术对水平像素轮廓进行直线矢量化处理，最终批量得到所有建筑物的水平矢量轮廓。本发明适用于卫星影像中俯视图为直线段结构多边形的普通建筑物的水平矢量轮廓的批量快速提取。

Description

卫星影像中建筑物水平矢量轮廓自动批量提取方法

技术领域

本发明涉及一种从单幅建筑群卫星影像中自动批量提取建筑物水平矢量轮廓的方法，特别是对俯视图为直线结构多边形的建筑物的水平矢量轮廓的自动批量提取。

背景技术

利用单幅卫星影像实现三维场景虚拟重建是一个很活跃的研究课题，它主要应用在城市建设规划、军事场景仿真、资源管理、抗震救灾模拟等方面。在真实场景的三维虚拟重建中，绝大部分为结构简单、俯视图为直线段结构多边形的普通建筑物，如何实现这类大量存在的普通建筑物的快速建模是高效重建三维建筑群的关键，而如何实现从建筑群影像中自动批量地提取建筑物的水平轮廓是高效重建三维建筑群的基础，并且决定着未来重建的虚拟场景中建筑群的分布格局与真实场景的符合程度。

当前的研究中，针对不同的应用目的，人们提出了各种各样的建筑物图像边缘提取算法，比较常见的有snake模型法、水平集曲线演化法、区域生长法等。以上方法的研究对象均为灰度图像，不适用于彩色图像，因此无法利用彩色图像中丰富的颜色信息。除此之外，snake模型法对初始位置敏感，需要依赖其他机制将初始轮廓放置在感兴趣的图像特征附近，否则轮廓提取会失败，目前大多采用手工选取的办法来设置初始边界，不仅非常繁琐，而且给轮廓线自动生成造成了困难。而较为常用的区域生长法存在着两点不足：第一，是初始种子点的选取问题。当前方法对种子点的选取大多数是手动选取，需要大量的人工干涉，耗时费力，效率很低。第二，是生长阈值的选取问题。生长阈值过大将会出现过分割，即生成的目标区域面积往往比实际面积大；而生长阈值过小又会导致分割不足，即目标区域生长不完整。因此不同建筑物需要选取不同的生长阈值，这一工作需要人工逐一根据经验完成。另外，snake模型法与区域生长法针对的是单个目标区域，也即每次只能生成一个目标区域，导致了方法的低效率。

以上方法提取的边缘轮廓，均由像素组成，这里称之为像素轮廓。本专利研究建筑物水平轮廓提取的目的，是为后续建筑群的三维虚拟重建提供建模数据，需要对提取出的建筑物像素轮廓进一步作直线矢量化处理，这里将直线矢量化处理后的轮廓称为矢量轮廓。目前的轮廓矢量化方法是将边缘像素提取出来并组成列表，然后对该列表进行直线段拟合，其中设置一个距离阈值，定义为偏离直线的最大距离，当超过该距离阈值后，直线段就会按比例一分为二。这种方法的结果受距离阈值选取的影响较大，为了得到较好的效果，影像中的每个建筑物都需要对应不同的距离阈值，以往设置距离阈值的方法往往是根据经验手工设置，操作效率低。

发明内容

为了克服现有技术的不足，本发明提供一种从单幅建筑群卫星影像中批量快速提取建筑物水平像素轮廓的自动化方法；进一步提出一种新的无需距离阈值及其他阈值选取的直线矢量化自动处理方法，对提取到的建筑物水平像素轮廓进行直线矢量化处理，得到满足后续建筑物三维虚拟重建所需的水平矢量轮廓。该方法特别适用于卫星影像中俯视图为直线段结构多边形的普通建筑物的水平矢量轮廓的批量快速提取，对航拍建筑群影像中的这类建筑群同样适用。

本发明解决其技术问题所采用的技术方案是：首先利用K-均值算法进行图像分类，得到建筑物的主干部分，取各建筑物区域质心解决初始种子点的选取问题。生长出所有种子点所在区域后，借助图像边缘信息将建筑物区域与周围区域分开，并根据矩形度、带状指数等特征参数去除掉非建筑物区域，从而实现建筑物水平像素轮廓的自动提取。然后利用Hough变换、块处理等技术对水平像素轮廓进行直线矢量化处理，最终批量得到所有建筑物的水平矢量轮廓。

本发明中，图像坐标系的原点位于图像左上角，x轴正方向竖直向下，y轴正方向水平向右。具体步骤如下：

步骤一、采用高斯平滑滤波器对建筑群卫星影像进行平滑处理，所述高斯平滑滤波器的方差σ＝0.8，窗口大小为7*7像素；然后在Hsv色彩空间内，保持每个像素的色调不变，按照步进值Δs调整饱和度，按照步进值Δv调整亮度，对卫星影像进行图像增强处理，得到源图像I_orig，其中Δs为0.06-0.08，Δv为0.05-0.07；

步骤二、轮廓自动提取，包括以下步骤：

1.不同建筑物种子点的自动生成，包括以下步骤：

1）在Lab色彩空间利用K-均值算法，在代表颜色信息的a、b色彩子空间上对源图像I_orig进行建筑物区域分割，K-均值算法执行中需要的分类类别数等于图像I_orig在a、b色彩子空间的二维直方图中的峰值个数，将分割出建筑物的图像存储为I_seg；

2）将图像I_seg转换为灰度图像，再由灰度图像转换为二值化图像I_bw，其中，建筑物区域为白色，背景为黑色，对二值化图像I_bw依次进行下列处理：

（1）采用窗口为5*5像素的中值滤波器去掉图像中孤立的像素点；

（2）采用5*5像素的方形结构算子进行开闭运算，去掉面积小于方形结构算子覆盖范围的区域；

（3）去除面积小于s_area个像素的连通区域，s_area的取值范围为[700，2000]内的整数；

（4）去除面积大于I_area个像素的连通区域，I_area的取值范围为[8000，10000]内的整数；

3）对建筑物区域分别编号，并记录其质心作为各建筑物的初始种子点；设种子点位置坐标为(x_k,y_k)，k为建筑物编号，x_k和y_k分别为第k个建筑物种子点在图像中的行号和列号，i表示图像的行号，j表示图像的列号，f(i,j)表示(i,j)处像素的灰度值，D表示联通区域，符号

表示向下取整运算；

2.提取建筑物区域，包括以下步骤：

1）利用区域生长法，生成各候选建筑物的区域，具体步骤如下：

（1）建立一个与源图像I_orig一样大小、像素灰度全部为0的图像I_out，设生长阈值为thresh，其取值范围为8-10个像素；

（2）定位出一个初始种子点(x_k,y_k)，k为建筑物编号，设置I_out(x_k,y_k)＝1；

（3）以(x_k,y_k)为中心，考虑(x_k,y_k)的四邻域像素(x_k-1,y_k)，(x_k+1,y_k)，(x_k,y_k-1)，(x_k,y_k+1)，在Lab颜色空间内，根据欧氏距离分别计算种子点邻域内像素与种子点的颜色差 $d(x,y)=\sqrt{{(L(x,y)-L_{\mathrm{seed}})}^2+{(a(x,y)-a_{\mathrm{seed}})}^2+{(b(x,y)-b_{\mathrm{seed}})}^2},$ 其中，L(x,y)、a(x,y)、b(x,y)分别表示邻域内当前像素(x,y)的L、a、b分量值，而L_seed、a_seed、b_seed则表示当前种子点的L、a、b分量值；如果邻域内当前像素（x,y）的I_out(x,y)＝0，并且d(x,y)<thresh，则设置I_out(x,y)＝1，同时将(x,y)压入堆栈；

（4）从堆栈中弹出一个像素当作新的种子点，赋值给(x_k,y_k)，回到步骤（3）；当堆栈为空时，生长结束；

（5）回到步骤（2），定位出下一个种子点，重复步骤（3）和（4），直到所有种子点扫描完成，整个生长过程结束；

2）对源图像I_orig灰度化后用Canny算子检测出强边缘，强边缘上的像素为白色，其他像素为黑色，将此含有强边缘信息的图像存储为I_edge；

3）遍历图像I_out，当I_out中的像素处于白色的建筑物区域时，查看图像I_edge中对应位置的像素是否处于白色的边缘，若是，则将I_out中对应像素变为黑色；

4）在图像I_out中，计算矩形度

长宽比

和带状指数

其中A，l分别是各区域的面积和周长，A_e，h，w分别对应各区域最小外接矩形的面积、长度和宽度；当d>0.7且s<3且r<5时，确定该区域为建筑物，被确认为是建筑物的区域将保留下来，否则将该区域变为黑色；

步骤三、利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号；对分割出来的单个建筑物图像I_{out_k}做如下处理：

（1）在I_{out_k}中，利用Hough变换检测出建筑物水平像素轮廓中对应的最长直线段，作为建筑物的主轴p；以主轴一端为定点，将水平像素轮廓旋转至主轴与图像坐标系的一个坐标轴平行，具体做法如下：

设theta角为建筑物主轴p与图像坐标系中x轴正方向的夹角，逆时针方向为正，单位为度；当|cos(theta)|＞cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的x轴平行；当|cos(theta)|≤cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的y轴平行；

（2）将旋转后的水平像素轮廓按其外界矩形在x，y方向上等间隔分块，每块大小为5*3个像素，定义为单元块，记做unit；在旋转后主轴所在的方向上将单元块定义为3个像素，另一方向定义为5个像素；最后一行和最后一列不足一个单元块的按一个单元块计算；

（3）在图像I_{out_k}中，从上到下，从左到右逐行扫描各单元块，计算各单元块中建筑物区域面积与单元块面积的百分比

其中，A_build是单元块中建筑物区域的面积，A_unit为单元块的面积，面积单位均为像素个数；当Percentage大于等于0.45时，认为该单元块是建筑物区域，保留该单元块并全部填充为白色，否则丢弃该单元块；

（4）将得到的直线化轮廓绕主轴的定点旋转-theta角度调整回原始位置，得到单个建筑物的直线化轮廓；

（5）分别对所有的建筑物区域进行上述步骤（1）-（4）的操作，并将处理后的各建筑物直线化轮廓合并到同一图像中，并存储为矢量图，最终完成对原始影像上建筑群中所有建筑物水平矢量轮廓的批量提取。

本发明的有益效果是：解决了现有基于卫星影像的建筑群水平像素轮廓提取和像素轮廓直线矢量化过程中半自动化和手动方法的缺陷，提出了自动的智能处理技术。首先，提出了一种自动生成不同建筑物初始种子点的方法，并利用区域生长的过分割和影像边缘信息分离建筑物与周围物体，得到建筑物的水平像素轮廓，克服了现有区域生长法中，初始种子点和生长阈值需要人工逐一根据经验选取的弊端。其次，为满足后续三维虚拟重建的需要，对提取到的水平像素轮廓自动进行直线矢量化处理，利用Hough变换检测出建筑物像素轮廓中的主轴，配合块处理技术进行直线矢量化，从而得到建筑物的水平矢量轮廓，避开了传统方法中不同区域需要根据经验手工设置不同距离阈值的不足。本发明可以对卫星影像中大量存在的俯视图为直线结构多边形的建筑物的水平矢量轮廓实现批量快速提取，为后续建筑群三维重建高效地提供建模数据。对航拍建筑群影像中的这类建筑同样适用。可广泛应用于城市建设规划、军事场景仿真等诸多三维建模和虚拟仿真领域。

附图说明

图1是建筑物水平像素轮廓主轴p与x轴的夹角theta示意图，p由Hough变换检测得到。

图2是水平像素轮廓直线矢量化示意图。

图3是本发明方法的流程图。

图中，1–建筑物水平像素轮廓，2–建筑物水平像素轮廓的主轴p（此处左上端为p的定点），3–主轴p与x轴的夹角theta，4–假想图像坐标系平移至原点与主轴定点重合，以示theta角构成，5–水平像素轮廓的外界矩形，6–单元块unit，7–5*3个像素。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明具体步骤如下：

步骤一、图像预处理

首先采用高斯平滑滤波器对建筑群卫星影像进行平滑处理，平滑过滤掉图像中存在的一些斑点噪声。这里采用方差σ＝0.8，窗口大小为7*7的高斯平滑滤波器对原影像进行滤波。然后在Hsv色彩空间内，保持每个像素的色调h值不变，按照s＝s+Δs调整饱和度，按照v＝v+Δv调整亮度，对卫星影像进行图像增强处理，使得建筑物颜色与周边景物的颜色对比度得到加强，便于区分建筑物与非建筑物区域，其中Δs的范围在0.06-0.08之间，Δv的范围在0.05-0.07之间。将预处理后用于后续处理的图像称为源图像，存储为I_orig。

步骤二、轮廓自动提取

1.不同建筑物种子点的自动生成

1）在Lab色彩空间利用K-均值算法，在代表颜色信息的a、b色彩子空间上对图像I_orig进行建筑物区域分割，K-均值算法执行中需要的分类类别数等于图像I_orig在a、b色彩子空间的二维直方图中的峰值个数。将分割出建筑物的图像存储为I_seg。之所以选择Lab色彩空间，是因为它更适合于人眼感觉，颜色更丰富。

2）通过K-均值算法完成区域的粗略分割后，卫星影像上一些与建筑物颜色相近的地面、车顶、道路等也被划分到建筑物的类别中，需要进一步处理。首先将图像I_seg由彩色图像转换为灰度图像，再由灰度图像转换为二值化图像，存储为I_bw。在二值化图像I_bw中，建筑物区域为白色，背景为黑色。对图像I_bw依次进行下列处理：

（1）中值滤波：采用窗口为5*5（单位：像素）的中值滤波器，去掉图像中孤立的像素点。

（2）开闭运算：采用5*5（单位：像素）的方形结构算子进行开闭运算，去掉面积小于方形结构算子覆盖范围的区域。

（3）去除小连通区域：去除面积小于s_area个像素的连通区域，主要目标是去除类似车顶、地面之类的连通区域。这里s_area的取值范围为[700，2000]内的整数。

（4）去除大连通区域：去除面积大于I_area个像素的连通区域；这里针对的对象主要是地面、道路等面积比较大的区域，I_area的取值范围为[8000，10000]内的整数。

通过上述一系列处理后，可以保留建筑物的主干部分，同时去除掉大部分非建筑物部分，其中残留的非建筑物部分将会在后续步骤2中的4）时去除。

3）最后对建筑物区域分别编号，并记录其质心作为各建筑物的初始种子点。设种子点位置坐标为(x_k,y_k)，k为建筑物编号，x_k和y_k分别为第k个建筑物种子点在图像中的行号和列号，这里将初始种子点的位置坐标作为(x_k,y_k)的初始值，计算方法如下：

上式中，i表示图像的行号，j表示图像的列号，f(i,j)表示(i,j)处像素的灰度值，D表示联通区域，符号

表示向下取整运算。

2.提取建筑物区域

主要包括两个部分：使用上述步骤1中得到的建筑物初始种子点发起区域生长；利用源图像I_orig的边缘信息分离建筑物区域和非建筑物区域，避开以往区域生长法需人工对不同建筑物选取不同生长阈值的不足。具体分为以下几步：

1）利用区域生长法，生成各候选建筑物的区域，通过生长阈值thresh的选择保证候选建筑物区域呈现过分割的情况，以使候选建筑物区域得以完整保留。具体步骤如下：

（1）首先建立一个与源图像I_orig一样大小、像素灰度全部为0的图像I_out，设生长阈值为thresh，其取值范围为8-10个像素。

（2）由公式（1）和公式（2）定位出一个初始种子点(x_k,y_k)，k为建筑物编号，设置I_out(x_k,y_k)＝1。

（3）以(x_k,y_k)为中心，考虑(x_k,y_k)的四邻域像素(x_k-1,y_k)，(x_k+1,y_k)，(x_k,y_k-1)，(x_k,y_k+1)，在Lab颜色空间内，根据欧氏距离分别计算种子点邻域内像素与种子点的颜色差d(x,y)，计算方法为：

$d(x,y)=\sqrt{{(L(x,y)-L_{\mathrm{seed}})}^2+{(a(x,y)-a_{\mathrm{seed}})}^2+{(b(x,y)-b_{\mathrm{seed}})}^2}---\left(3\right)$

上式中，L(x,y)、a(x,y)、b(x,y)分别表示邻域内当前像素(x,y)的L、a、b分量值，而L_seed、a_seed、b_seed则表示当前种子点的L、a、b分量值。

如果邻域内当前像素（x,y）的I_out(x,y)＝0，并且d(x,y)<thresh，则设置I_out(x,y)＝1，同时将(x,y)压入堆栈。

（4）从堆栈中弹出一个像素，把它当作新的种子点，赋值给(x_k,y_k)，回到上一步（3）。当堆栈为空时，生长结束。结果使得I_out在初始种子点附近逐渐形成了建筑物区域（白色部分）。

（5）回到步骤（2），定位出下一个种子点，重复以上步骤（3）和（4），直到所有种子点扫描完成，整个生长过程结束。期间生长阈值和种子点邻域的选取方式保持不变。

与源图像I_orig相比照，上述过程中得到的图像I_out表现为候选建筑物区域出现过分割，区域边缘的一些非建筑物区域也同时生成，这种方式保证了候选建筑物区域的完整性，可以确保后续操作对候选建筑物区域的完整提取。

2）对源图像I_orig灰度化后用Canny算子检测出强边缘，强边缘上的像素为白色，其他像素为黑色，将此含有强边缘信息的图像存储为I_edge。

3）遍历图像I_out，当I_out中的像素处于白色的建筑物区域时，查看图像I_edge中对应位置的像素是否处于白色的边缘，若是，则将I_out中对应像素变为黑色，即背景的颜色。经过上述处理后，在图像I_out中候选建筑物将会与周围区域分离，在图像中表现为大片的候选建筑物区域伴随着琐碎、零散的非建筑物区域。

4）在图像I_out中，利用矩形度d、长宽比r、带状指数s等特征参数去除掉非建筑物区域和候选建筑物区域中的非建筑物。各特征参数的计算方法分别为：

$d=\frac{A}{A_e},r=\frac{h}{w},s=\frac{l}{4\sqrt{A}}$

其中A，l分别是各区域的面积和周长，A_e，h，w分别对应各区域最小外接矩形的面积、长度和宽度。

由于建筑群卫星影像中，普通建筑物的俯视图一般都是直线段结构的封闭多边形，具有规则矩形或者类似矩形的特征，因此会具有较高矩形度、较低的带状指数以及一定的长宽比，对建筑物通常满足：矩形度d大于0.7，带状指数s小于3.0，长宽比r小于5.0。根据上述三个特征参数可以去除掉非建筑物区域和候选建筑物区域中的非建筑物。判断准则为：

当d>0.7且s<3且r<5时，确定该区域为建筑物。

被确认为是建筑物的区域将保留下来，否则将该区域变为黑色，即背景的颜色。

步骤三、水平像素轮廓的直线矢量化处理

经过上面步骤后得到的建筑物区域边缘部分像素会参差不齐，原本该为直线的边缘，变成了由许多琐碎小线段构成的折线段，这样的情况不利于后续的三维建模工作的进行。在三维建模中，建筑物水平轮廓中小线段增多会导致三维建筑模型中三角形片面数的激增，直接影响虚拟仿真***运行的实时性，尤其在大场景仿真中，影响更加严重。因此本部分利用Hough变换、块处理等技术对建筑物的水平像素轮廓进行直线矢量化处理。

首先，在步骤二的基础上，利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号。对分割出来的单个建筑物图像I_{out_k}做如下处理：

（1）在I_{out_k}中，利用Hough变换检测出建筑物水平像素轮廓中对应的最长直线段，作为建筑物的主轴p。以主轴一端为定点，将水平像素轮廓旋转至主轴与图像坐标系的一个坐标轴平行，具体做法如下：

设theta角为建筑物主轴p与图像坐标系中x轴正方向的夹角，逆时针方向为正，单位为度。

当|cos(theta)|＞cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的x轴平行；

当|cos(theta)|≤cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的y轴平行。

（2）将旋转后的水平像素轮廓按其外界矩形在x，y方向上等间隔分块，每块大小为5*3个像素，定义为单元块，记做unit。由于通常主轴方向上会有阳台等突出物这一类的细节信息，因此在旋转后主轴所在的方向上将单元块定义为3个像素，以免丢失细节信息，另一方向定义为5个像素。最后一行和最后一列不足一个单元块的按一个单元块计算。

（3）在图像I_{out_k}中，从上到下，从左到右逐行扫描各单元块，计算各单元块中建筑物区域面积与单元块面积的百分比：

$\mathrm{Percentage}=\frac{A_{\mathrm{build}}}{A_{\mathrm{unit}}}$

其中，A_build是单元块中建筑物区域的面积，A_unit为单元块的面积，面积单位均为像素个数。

当Percentage大于等于0.45时，认为该单元块是建筑物区域，保留该单元块并全部填充为白色，否则丢弃该单元块。扫描结束后，最后保留的所有单元块组成了新的建筑物区域，跟踪区域的边缘就可以得到直线化的轮廓。

（4）将得到的直线化轮廓绕主轴的定点旋转-theta角度调整回原始位置，得到单个建筑物的直线化轮廓。

（5）分别对所有的建筑物区域进行上述（1）-（4）的操作，并将处理后的各建筑物直线化轮廓合并到同一图像中，并存储为矢量图，最终完成对原始影像上建筑群中所有建筑物水平矢量轮廓的批量提取。

本实例中，卫星影像建筑群的提取包括以下三个步骤：

步骤一、图像预处理

首先采用方差σ＝0.8、窗口大小为7*7的高斯平滑滤波器，平滑过滤掉建筑群卫星影像中存在的一些斑点噪声。然后在Hsv色彩空间内，对卫星影像作如下调整：每个像素色调h不变，饱和度s＝s+Δs，亮度v＝v+Δv，其中Δs的范围在0.06-0.08之间，Δv的范围在0.05-0.07之间，完成对影像的图像增强处理。将预处理后用于后续处理的图像称为源图像，存储为I_orig。

步骤二、轮廓自动提取

详细操作步骤如下：

1.不同建筑物种子点的自动生成

1）在Lab色彩空间利用K-均值算法，在代表颜色信息的a、b色彩子空间上对图像I_orig进行建筑物区域分割，将得到的图像存储为I_seg。其中K-均值算法需要的分类类别数等于图像I_orig在a、b色彩子空间的二维直方图中的峰值个数。

2）将图像I_seg由彩色图像转换为灰度图像，再由灰度图像转换为二值化图像，存储为I_bw，其中建筑物区域为白色，背景为黑色。对图像I_bw依次进行下列处理，去除掉大部分的非建筑区域：

（1）中值滤波：采用窗口为5*5（单位：像素）的中值滤波器，去掉图像中孤立的像素点。

（2）开闭运算：采用5*5（单位：像素）的方形结构算子进行开闭运算，去掉面积小于方形结构算子覆盖范围的区域。

（3）去除小连通区域：去除面积小于s_area个像素的连通区域，主要目标是去除类似车顶、地面之类的连通区域。这里s_area的取值范围为[700，2000]内的整数。

（4）去除大连通区域：去除面积大于I_area个像素的连通区域；这里针对的对象主要是地面、道路等面积比较大的区域，I_area的取值范围为[8000，10000]内的整数。

通过上述一系列处理后，少数残留的非建筑物部分将会在后续步骤2中的4）时去除。

3)最后对建筑物区域分别编号，并记录其质心作为各建筑物的初始种子点。设种子点位置坐标为(x_k,y_k)，k为建筑物编号，x_k和y_k分别为第k个建筑物种子点在图像中的行号和列号，这里将初始种子点的位置坐标作为(x_k,y_k)的初始值，计算方法如下：

上式中，i表示图像的行号，j表示图像的列号，f(i,j)表示(i,j)处像素的灰度值，D表示联通区域，符号

表示向下取整运算。

2.提取建筑物区域

1）基于上述步骤1中得到的建筑物初始种子点，利用区域生长法，生成各候选建筑物的区域，通过生长阈值thresh的选择保证候选建筑物区域呈现过分割的情况，以使候选建筑物区域得以完整保留。具体步骤如下：

（1）首先建立一个与源图像I_orig一样大小、像素灰度全部为0的图像I_out，设生长阈值为thresh，其取值范围为8-10个像素。

（2）由公式（1）和公式（2）定位出一个初始种子点(x_k,y_k)，k为建筑物编号，设置I_out(x_k,y_k)＝1。

（3）以(x_k,y_k)为中心，考虑(x_k,y_k)的四邻域像素(x_k-1,y_k)，(x_k+1,y_k)，(x_k,y_k-1)，(x_k,y_k+1)，在Lab颜色空间内，根据欧氏距离分别计算种子点邻域内像素与种子点的颜色差d(x,y)，计算方法为：

$d(x,y)=\sqrt{{(L(x,y)-L_{\mathrm{seed}})}^2+{(a(x,y)-a_{\mathrm{seed}})}^2+{(b(x,y)-b_{\mathrm{seed}})}^2}---\left(3\right)$

上式中，L(x,y)、a(x,y)、b(x,y)分别表示邻域内当前像素(x,y)的L、a、b分量值，而L_seed、a_seed、b_seed则表示当前种子点的L、a、b分量值。

如果邻域内当前像素（x,y）的I_out(x,y)＝0，并且d(x,y)<thresh，则设置I_out(x,y)＝1，同时将(x,y)压入堆栈。

（4）从堆栈中弹出一个像素，把它当作新的种子点，赋值给(x_k,y_k)，回到上一步（3）。当堆栈为空时，生长结束。结果使得I_out在初始种子点附近逐渐形成了建筑物区域（白色部分）。

（5）回到步骤（2），定位出下一个种子点，重复以上步骤（3）和（4），直到所有种子点扫描完成，整个生长过程结束。期间生长阈值和种子点邻域的选取方式保持不变。

与源图像I_orig相比照，上述过程中得到的图像I_out表现为候选建筑物区域出现过分割，以确保后续操作对候选建筑物区域的完整提取。

2）对源图像I_orig灰度化后用Canny边缘检测出强边缘，强边缘上的像素为白色，其他像素为黑色。将此含有强边缘信息的图像存储为I_edge。

3）遍历图像I_out，当I_out中是白色的建筑物区域时，查看图像I_edge中对应位置是否为白色的边缘，当为白色边缘时，将I_out中对应像素变为黑色，即背景的颜色。这样，在图像I_out中候选建筑物将会与周围区域分离，在图像中表现为大片的候选建筑物区域伴随着琐碎、零散的非建筑物区域。

4）在图像I_out中，利用矩形度d、长宽比r、带状指数s等特征参数去除掉非建筑物区域和候选建筑物区域中的非建筑物。各特征参数的计算方法分别为：

$d=\frac{A}{A_e},r=\frac{h}{w},s=\frac{l}{4\sqrt{A}}$

其中A，l分别是各区域的面积和周长，A_e，h，w分别对应各区域最小外接矩形的面积、长度和宽度。

根据准则“当连通域的d>0.7且s<3且r<5时，确定该区域为建筑物。”判断I_out中的各区域是否为建筑物，如果是建筑物，则保留该区域；否则将该区域变为黑色，即背景的颜色。

步骤三、水平像素轮廓的直线矢量化处理

本部分利用Hough变换、块处理等技术对建筑物的水平像素轮廓进行直线矢量化处理。首先，在图像I_out中，利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号。对图像I_{out_k}分别做如下处理：

（1）在I_{out_k}中，利用Hough变换检测出建筑物水平像素轮廓中对应的最长直线段，作为建筑物的主轴p。以主轴一端为定点，将水平像素轮廓旋转至主轴与图像坐标系的一个坐标轴平行。设theta角为建筑物主轴p与图像坐标系中x轴正方向的夹角，逆时针方向为正，单位为度。theta等于90°表示该建筑物在图像中水平放置，theta等于0°表示垂直放置。当theta既不是0°，也不是90°时，建筑物在图像中表现为倾斜放置，因此需要对其进行旋转操作，使其水平或垂直放置：

当|cos(theta)|≥cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的x轴平行；

当|cos(theta)|≤cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的y轴平行。

（2）将旋转后的水平像素轮廓按其外界矩形在x，y方向上等间隔分块，每块大小为5*3个像素，定义为单元块，记做unit。由于通常主轴方向上会有阳台等突出物这一类的细节信息，因此在旋转后主轴所在的方向上将单元块定义为3个像素，以免丢失细节信息，另一方向定义为5个像素。最后一行和最后一列不足一个单元块的按一个单元块计算。

（3）在图像I_{out_k}中，从上到下，从左到右逐行扫描各单元块，计算各单元块中建筑物区域面积与单元块面积的百分比：

$P_{\mathrm{ercentage}}=\frac{A_{\mathrm{build}}}{A_{\mathrm{unit}}}$

其中，A_build是单元块中建筑物区域的面积，A_unit为单元块的面积，面积单位均为像素个数。

当Percentage大于等于0.45时，认为该单元块是建筑物区域，保留该单元块并全部填充为白色，否则丢弃该单元块。扫描结束后保留的所有单元块组成了新的建筑物区域，跟踪区域的边缘就可以得到直线化的轮廓。

（4）将得到的直线化轮廓绕主轴的定点旋转-theta角度调整回原始位置，得到单个建筑物的直线化轮廓。

（5）分别对所有的建筑物区域进行上述（1）-（4）的操作，并将处理后的各建筑物直线化轮廓合并到同一图像中，并存储为矢量图，最终完成对原始影像上建筑群中所有建筑物水平矢量轮廓的批量提取。

实施例1：对原始图像为一张702*902像素的北京某地区卫星影像进行建筑物轮廓的提取。

采用方差σ＝0.8，窗口大小为7*7的高斯平滑滤波器对原影像进行平滑滤波；在Hsv空间内，保持每个像素的h值不变，v=v+0.08，s=s+0.07进行图像增强处理，使得建筑物与背景的颜色对比度加强。将处理后的图像存储为I_orig，作为后续处理的源图像。

图像I_orig在Lab色彩空间的a、b色彩子空间中峰值个数为3，作为用K-均值算法对其进行建筑物区域分割时的类别数，分割后的图像保存为I_seg，此实施例中，建筑物区域存在于第3类。然后对图像I_seg进行灰度化和二值化处理，得到其二值化图像I_bw，此时I_bw中的连通区域个数为100，其中包括了建筑物与非建筑物区域。接着遍历每个连通区域，去除掉面积小于s_area=2000和面积大于I_area=10000的区域，利用大小为5*5（单位：像素）的方形结构算子对图像进行开运算，去掉面积小于方形结构算子覆盖范围的区域，此时二值图像I_bw中连通区域的个数变成了9，即建筑物的主干部分有9个。分别计算每个连通区域的质心坐标作为各建筑物的初始种子点坐标，如表1所示：

表1各建筑物的初始种子点坐标

区域编号	种子点x坐标	种子点y坐标
			1	100	352
2	245	479
			3	319	107
4	345	700
			5	352	217
6	463	230
			7	540	739
8	563	314
			9	655	366

利用上面求得的9个初始种子点，选取生长阈值thresh=10，利用区域生长法得到9个候选建筑物区域，图像存为I_out。结合Canny算子从源图像I_orig中检测出的强边缘，将图像I_out中的候选建筑物区域与非建筑区域分离。在图像I_out中，分别计算各区域的矩形度d、带状指数s、长宽比r三个特征参数，其中候选建筑物区域的特征参数如表2所示（由于非建筑物区域琐碎、零散、无规律且易于区分，这里不再列举。）：

表2九个候选建筑物区域的特征参数

区域编号	矩形度d	带状指数s	长宽比r
				1	0.75	0.0029	1.9
2	0.77	0.0027	2.26
				3	0.83	0.0035	1.98
4	0.71	0.003	2.83
				5	0.76	0.004	2.49
6	0.82	0.0031	1.91
				7	0.76	0.0028	2.53
8	0.85	0.003	2.06
				9	0.55	0.0093	1.12

根据准则“当d>0.7且s<3且r<5时，确定该区域为建筑物。”判断各区域是否为建筑物。其中区域9不满足指标，因此舍去这个区域，至此初步得到了图像中所有8个建筑物的水平像素轮廓。

利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号。利用Hough变换检测出每个I_{out_k}中建筑物水平像素轮廓的主轴p，分别计算水平像素轮廓主轴与x轴的夹角theta，根据条件|cos(theta)|≥cos(45°)和|cos(theta)|≤cos(45°)将水平像素轮廓旋转至水平或者垂直位置，方便下一步操作。

表3各建筑物水平像素轮廓主轴与x轴的夹角theta

区域编号	主轴与x轴的夹角theta
		1	-42.02°
2	-41.28°

3	-41.11°
		4	-30.37°
5	-22.38°
		6	-40.99°
7	-30.04°
		8	-41.66°

此实例中，表3中建筑物都处于倾斜位置，因此均需要对其进行旋转。利用大小为5*3个像素的单元块unit分别对单个旋转后的建筑物区域进行划分，判断每个单元块中建筑物区域（白色部分）占整个单元块区域面积的百分比Percentage，当Percentage>0.45时，即认为该单元块是建筑物区域并全部填充为白色，处理完全部单元块后，将建筑物水平像素轮廓旋转回初始位置。依次对所有建筑物区域做如上操作，跟踪区域的边缘就可以得到直线化后的建筑物轮廓，最后将各建筑物直线化轮廓合并到同一图像中，并存储为矢量图。

实施例2：对原始图像为一张775*401像素的普渡大学宿舍区的卫星影像进行建筑物轮廓的提取。

采用方差σ＝0.8，窗口大小为7*7的高斯平滑滤波器对原影像进行平滑滤波；在Hsv空间内，保持每个像素的h不变，v=v+0.06，s=s+0.05进行图像增强处理，使得建筑物与背景的颜色对比度加强。将处理后的图像存储为I_orig，作为后续处理的源图像。

图像I_orig在Lab色彩空间的a、b色彩子空间中峰值个数为3，作为用K-均值算法对其进行建筑物区域分割时的类别数，分割后的图像保存为I_seg，此实施例中，建筑物区域存在于第2类。然后对图像I_seg进行灰度化和二值化处理，得到其二值化图像I_bw，此时I_bw中的连通区域个数为16，其中包括了建筑物与非建筑物区域。接着遍历每个连通区域，去除掉面积小于s_area=700和面积大于I_area=8000的区域，利用大小为5*5（单位：像素）的方形结构算子对图像进行开运算，去掉面积小于方形结构算子覆盖范围的区域，此二值图像I_bw中连通区域的个数变成了9，即建筑物的主干部分有9个。此时分别计算每个连通区域的质心坐标作为各建筑物的初始种子点坐标，如表4所示：

表4各建筑物的初始种子点坐标

区域编号	种子点x坐标	种子点y坐标
			1	108	45
2	96	194
			3	182	366
4	207	209
			5	359	179
6	437	58
			7	542	333
8	513	179
			9	665	118

利用上面求得的9个初始种子点，选取生长阈值thresh=8，利用区域生长法得到9个候选建筑物区域，图像存为I_out。结合Canny算子从源图像I_orig中检测出的强边缘，将图像I_out中的候选建筑物区域与非建筑区域分离。在图像I_out中，分别计算各区域的矩形度d、带状指数s、长宽比r三个特征参数，其中候选建筑物区域的特征参数如表5所示（由于非建筑物区域琐碎、零散、无规律且易于区分，这里不再列举。）：

表5九个候选建筑物区域的特征参数

区域编号	矩形度d	带状指数s	长宽比r
				1	0.85	0.003	4.75
2	0.83	0.0036	4.35
				3	0.92	0.0035	5.32

4	0.87	0.0032	4.46
				5	0.82	0.0032	4.11
6	0.97	0.006	1.93
				7	0.84	0.0033	4.7
8	0.86	0.0032	4.37
				9	0.94	0.0032	4.83

根据准则“当d>0.7且s<3且r<5时，确定该区域为建筑物。”判断各区域是否为建筑物。上述9个连通区域均符合建筑区域的参数指标，至此初步得到了图像中所有9个建筑物的水平像素轮廓。

利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号。利用Hough变换检测出每个I_{out_k}中建筑物水平像素轮廓的主轴p，分别计算水平像素轮廓主轴与x轴的夹角theta，根据条件|cos(theta)|≥cos(45°)和|cos(theta)|≤cos(45°)将水平像素轮廓旋转至水平或者垂直位置，方便下一步操作。

表6各建筑物水平像素轮廓主轴与x轴的夹角theta

区域编号	主轴与x轴的夹角theta
		1	90°
2	0°
		3	90°
4	0°
		5	90°
6	90°
		7	89.5°

8	0°
		9	90°

此实例中，表6中建筑物都处于水平或者垂直位置，因此不需要对其进行旋转。利用大小为5*3个像素的单元块unit分别对单个旋转后的建筑物区域进行划分，判断每个单元块中建筑物区域（白色部分）占整个单元块区域面积的百分比Percentage，当Percentage>0.45时，即认为该单元块是建筑物区域并全部填充为白色，处理完全部单元块后，将建筑物水平像素轮廓旋转回初始位置。依次对所有建筑物区域做如上操作，跟踪区域的边缘就可以得到直线化后的建筑物轮廓，最后将各建筑物直线化轮廓合并到同一图像中，并存储为矢量图。

实施例3：对原始图像为一张874*383像素的西安某地区的卫星影像进行建筑物轮廓的提取。

采用方差σ＝0.8，窗口大小为7*7的高斯平滑滤波器对原影像进行平滑滤波；在Hsv空间内，保持每个像素的h不变，v=v+0.07，s=s+0.06进行图像增强处理，使得建筑物与背景的颜色对比度加强。将处理后的图像存储为I_orig，作为后续处理的源图像。

图像I_orig在Lab色彩空间的a、b色彩子空间中峰值个数为5，作为用K-均值算法对其进行建筑物区域分割时的类别数，分割后的图像保存为I_seg，此实施例中，建筑物区域存在于第1类。然后对图像I_seg进行灰度化和二值化处理，得到其二值化图像I_bw，此时I_bw中的连通区域个数为84，其中包括了建筑物与非建筑物区域。接着遍历每个连通区域，去除掉面积小于s_area=1000和面积大于I_area=9000的区域，利用大小为5*5（单位：像素）的方形结构算子对图像进行开运算，去掉面积小于方形结构算子覆盖范围的区域，此二值图像I_bw中连通区域的个数变成了12，即建筑物的主干部分有12个。分别计算每个连通区域的质心坐标作为各建筑物的初始种子点坐标，如表7所示：

表7各建筑物的初始种子点坐标

区域编号

种子点x坐标

种子点y坐标

1	15	39
			2	37	146
3	93	72
			4	152	168
5	167	346
			6	239	195
7	404	192
			8	406	310
9	356	10
			10	566	253
11	778	308
			12	764	191

利用上面求得的12个初始种子点，选取生长阈值thresh=9，利用区域生长法得到12个候选建筑物区域，图像存为I_out。结合Canny算子从源图像I_orig中检测出的强边缘，将图像I_out中的候选建筑物区域与非建筑区域分离。在图像I_out中，分别计算各区域的矩形度d、带状指数s、长宽比r三个特征参数，其中候选建筑物区域的特征参数如表8所示（由于非建筑物区域琐碎、零散、无规律且易于区分，这里不再列举。）：

表8十二个候选建筑物区域的特征参数

区域编号	矩形度d	带状指数s	长宽比r
				1	0.71	0.0083	9.81
2	0.59	0.011	5.99
				3	0.65	0.005	2.35

4	0.85	0.0036	1.93
				5	0.7	0.0093	1.18
6	0.45	0.0046	2.5
				7	0.81	0.004	4.65
8	0.83	0.0039	4.62
				9	0.55	0.0093	1.12
10	0.89	0.0037	4.74
				11	0.89	0.0036	4.5
12	0.79	0.0046	3.36

根据准则“当d>0.7且s<3且r<5时，确定该区域为建筑物。”判断各区域是否为建筑物。其中区域1、2、3、5、6、9不满足指标，因此舍去这6个区域，剩余的6个连通区域均符合建筑区域的参数指标，至此初步得到了图像中所有6个建筑物的水平像素轮廓。

利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号。利用Hough变换检测出每个I_{out_k}中建筑物水平像素轮廓的主轴p，分别计算水平像素轮廓主轴与x轴的夹角theta，根据条件|cos(theta)|≥cos(45°)和|cos(theta)|≤cos(45°)将水平像素轮廓旋转至水平或者垂直位置，方便下一步操作。

表9各建筑物水平像素轮廓主轴与x轴的夹角theta

区域编号	主轴与x轴的夹角theta
		1	45°
2	-80.5°
		3	90°

4	0°
		5	-48.8°
6	90°

此实例中，表9中建筑物1、2、5处于倾斜位置，需要对其进行旋转，其余建筑物则不需要进行旋转处理。利用大小为5*3个像素的单元块unit分别对单个旋转后（如果需要）的建筑物区域进行划分，判断每个单元块中建筑物区域（白色部分）占整个单元块区域面积的百分比Percentage，当Percentage>0.45时，即认为该单元块是建筑物区域并全部填充为白色，处理完全部单元块后，将旋转过的建筑物水平像素轮廓再旋转回初始位置。依次对所有建筑物区域做如上操作，跟踪区域的边缘就可以得到直线化后的建筑物轮廓，最后将各建筑物直线化轮廓合并到同一图像中，并存储为矢量图。

Claims (2)

1.一种卫星影像中建筑物水平矢量轮廓自动批量提取方法，其特征在于包括下述步 

骤： 

步骤一、采用高斯平滑滤波器对建筑群卫星影像进行平滑处理，所述高斯平滑滤波器的方差σ＝0.8，窗口大小为7*7像素；然后在Hsv色彩空间内，保持每个像素的色调不变，按照步进值Δs调整饱和度，按照步进值Δv调整亮度，对卫星影像进行图像增强处理，得到源图像I_orig，其中Δs为0.06-0.08，Δv为0.05-0.07； 

步骤二、轮廓自动提取，包括以下步骤： 

1.不同建筑物种子点的自动生成，包括以下步骤： 

1）在Lab色彩空间利用K-均值算法，在代表颜色信息的a、b色彩子空间上对源图像I_orig进行建筑物区域分割，K-均值算法执行中需要的分类类别数等于图像I_orig在a、b色彩子空间的二维直方图中的峰值个数，将分割出建筑物的图像存储为I_seg； 

2）将图像I_seg转换为灰度图像，再由灰度图像转换为二值化图像I_bw，其中，建筑物区域为白色，背景为黑色，对二值化图像I_bw依次进行下列处理： 

（1）采用窗口为5*5像素的中值滤波器去掉图像中孤立的像素点； 

（2）采用5*5像素的方形结构算子进行开闭运算，去掉面积小于方形结构算子覆盖范围的区域； 

（3）去除面积小于s_area个像素的连通区域，s_area的取值范围为[700，2000]内的整数； 

（4）去除面积大于I_area个像素的连通区域，I_area的取值范围为[8000，10000]内的整数； 

3）对建筑物区域分别编号，并记录其质心作为各建筑物的初始种子点；设种子点位置坐标为(x_k,y_k)，k为建筑物编号，x_k和y_k分别为第k个建筑物种子点在图像中的行号和列号，

i表示图像的行号，j表示图像的列号，f(i,j)表示(i,j)处像素的灰度值，D表示联通区域，符号

表示向下取整运算。 

2.提取建筑物区域，包括以下步骤： 

1）利用区域生长法，生成各候选建筑物的区域，具体步骤如下： 

（1）建立一个与源图像I_orig一样大小、像素灰度全部为0的图像I_out，设生长阈值为thresh，其取值范围为8-10个像素； 

（2）定位出一个初始种子点(x_k,y_k)，k为建筑物编号，设置I_out(x_k,y_k)＝1； 

（3）以(x_k,y_k)为中心，考虑(x_k,y_k)的四邻域像素(x_k-1,y_k)，(x_k+1,y_k)，(x_k,y_k-1)，(x_k,y_k+1)，在Lab颜色空间内，根据欧氏距离分别计算种子点邻域内像素与种子点的颜色差其中，L(x,y)、a(x,y)、b(x,y)分别表示邻域内当前像素(x,y)的L、a、b分量值，而L_seed、a_seed、b_seed则表示当前种子点的L、a、b分量值；如果邻域内当前像素（x,y）的I_out(x,y)＝0，并且d(x,y)<thresh，则设置I_out(x,y)＝1，同时将(x,y)压入堆栈； 

（4）从堆栈中弹出一个像素当作新的种子点，赋值给(x_k,y_k)，回到步骤（3）；当堆栈为空时，生长结束； 

（5）回到步骤（2），定位出下一个种子点，重复步骤（3）和（4），直到所有种子点扫描完成，整个生长过程结束； 

2）对源图像I_orig灰度化后用Canny算子检测出强边缘，强边缘上的像素为白色，其他像素为黑色，将此含有强边缘信息的图像存储为I_edge； 

3）遍历图像I_out，当I_out中的像素处于白色的建筑物区域时，查看图像I_edge中对应位置的像素是否处于白色的边缘，若是，则将I_out中对应像素变为黑色； 

4）在图像I_out中，计算矩形度

长宽比

和带状指数

其中A，l分别是各区域的面积和周长，A_e，h，w分别对应各区域最小外接矩形的面积、长度和宽度；当d>0.7且s<3且r<5时，确定该区域为建筑物，被确认为是建筑物的区域将保留下来，否则将该区域变为黑色； 

步骤三、利用标记连通区域的方法将图像I_out中的单个建筑物分割出来，分别存储为I_{out_k}，k为建筑物编号；对分割出来的单个建筑物图像I_{out_k}做如下处理： 

（1）在I_{out_k}中，利用Hough变换检测出建筑物水平像素轮廓中对应的最长直线段，作为建筑物的主轴p；以主轴一端为定点，将水平像素轮廓旋转至主轴与图像坐 标系的一个坐标轴平行，具体做法如下： 

设theta角为建筑物主轴p与图像坐标系中x轴正方向的夹角，逆时针方向为正，单位为度；当|cos(theta)|＞cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的x轴平行；当|cos(theta)|≤cos(45°)时，将建筑物水平像素轮廓旋转至主轴与图像坐标系的y轴平行； 

（2）将旋转后的水平像素轮廓按其外界矩形在x，y方向上等间隔分块，每块大小为5*3个像素，定义为单元块，记做unit；在旋转后主轴所在的方向上将单元块定义为3个像素，另一方向定义为5个像素；最后一行和最后一列不足一个单元块的按一个单元块计算； 

（3）在图像I_{out_k}中，从上到下，从左到右逐行扫描各单元块，计算各单元块中建筑物区域面积与单元块面积的百分比

其中，A_build是单元块中建筑物区域的面积，A_unit为单元块的面积，面积单位均为像素个数；当Percentage大于等于0.45时，认为该单元块是建筑物区域，保留该单元块并全部填充为白色，否则丢弃该单元块； 

（4）将得到的直线化轮廓绕主轴的定点旋转-theta角度调整回原始位置，得到单个建筑物的直线化轮廓； 

（5）分别对所有的建筑物区域进行上述步骤（1）-（4）的操作，并将处理后的各建筑物直线化轮廓合并到同一图像中，并存储为矢量图，最终完成对原始影像上建筑群中所有建筑物水平矢量轮廓的批量提取。 

Images