CN104217440B - 一种从遥感图像中提取建成区的方法 - Google Patents

Abstract

本发明公开了一种从遥感图像中提取建成区的方法。该方法可采用基于贝叶斯推理的视觉显著性检测方法，结合图像分割、自动阈值选取实现高分辨率遥感图像建成区的自动检测和提取，可广泛应用于城市规划、城市扩张研究、灾情评估和救灾决策等多个领域。

Description

一种从遥感图像中提取建成区的方法

技术领域

本发明涉及图像处理领域，具体地，涉及一种从遥感图像中提取建成区的方法。

背景技术

从遥感图像中提取建成区是将遥感图像应用于城市研究的一个关键环节，特别是对城市规划、城市扩张研究、城区受灾范围及损失评估等有重要意义，也是目前遥感图像处理与分析领域的研究热点之一。近年来，随着经济快速发展，我国的城市化进程也进入了全新的阶段。城市扩张的特征之一是城市面积的变化，尤其是土地利用性质的变化。合理进行城市规划和有效利用土地资源是使城市化合理有序进行以及可持续发展的前提。卫星遥感图像是目前城市研究的重要来源，如何利用遥感图像来界定城市边界、提取城市区域是利用影像进行城市研究的前提。目前，国内外很多学者进行了这方面的研究，但是城市是一个复合体，具有很多不确定性，每个城市都有其特点，针对一个城市进行的研究成果，很难适用于另一个城市。因此，建成区的准确提取仍是一个世界性技术难题，需对其做进一步研究和探索。

目前，基于遥感影像的城市建成区的提取方法主要有：目视解译法、非监督分类法、监督分类法、人工建筑指数法、基于光谱知识的城区提取模型等。目视解译法靠人工确定影像上的物体属性或特征，对影像进行识别，描述影像上的各种关系，根据影像上的类别、属性和关系进行系统的解释。这种方法由于加入了人的识别，判断精确度较高，但工作较繁琐，效率低下，同时需要解译人员具有较高的遥感图像解译识别知识与经验。非监督分类的前提是假定遥感图像上同类物体在同样条件下具有相同的光谱特征条件。非监督分类不必获取影像地物的先验知识，仅靠影像上不同类的物光谱信息(或纹理信息、几何信息)进行特征提取再通过非监督分类方法或特征分布模型进行分类，最后对已分出的各个类别的属性进行确认。与非监督分类方法相对应地，监督分类法需要大量的已知建成区数据进行训练，将训练好的模型应用于新的影像时，可以直接从影像中提取建成区。这种方法虽然较非监督分类方法更为精确，但是需要大量的人工标注。归一化建筑指数提取建成区是以Landsat TM数据为基础利用第5和第4个波段的反射差来计算建筑指数，然后再用阈值分割的方法提取建成区。多光谱或者高光谱数据波段众多，光谱精细，为遥感图像分析提供了新的思路。基于光谱知识的建成区提取模型综合利用多光谱或者高光谱数据的多个波段，分析建成区在各个波段的反射率，建立光谱反射模型，或用光谱匹配或用规则来提取建成区。

上述方法在实际应用(例如，针对特定数据类型的应用)中仍然需要人工干预和背景知识的加入。现有方法的普适性，全自动化能力仍然存在不足，需要研究新的方法。

发明内容

本发明的目的是提供一种方法，该方法能从如资源三号、高分一号、Quickbird等高分辨率遥感图像中自动提取建成区。

为了实现上述目的，本发明提供一种从遥感图像中提取建成区的方法，包括：获取图像中每个像素的显著性值；根据显著性值将图像划分为前景区域和背景区域；分别计算前景区域和背景区域的像素特征似然函数；针对图像中的每个像素，根据该像素的显著性值和像素特征似然函数得到该像素为建筑区的概率值；根据所述概率值，判断图像中的每个像素是否为建成区；以及从图像中对应于被判断为建成区的区域中提取建成区。

优选地，本发明可基于图像中的边缘像素来计算每个像素的显著性值，并对显著性值进行归一化得到先验概率。

优选地，本发明可对图像进行超像素分割，并以超像素为单位将图像划分为前景区域和背景区域。

优选地，本发明可采用自动阈值分割以将图像划分为前景区域和背景区域。

优选地，可计算前景区域和背景区域中的像素特征似然函数，并且可基于贝叶斯公式并利用先验概率和像素特征似然函数计算像素为建成区的后验概率，然后可采用自动阈值分割进行建成区判断，最后从被判断为建成区的区域中提取建成区。

通过上述技术方案，可结合图像分割、自动阈值分割、特征似然概率分析技术自动从遥感图像中提取建成区。在上述过程中，不需要进行数据训练和人工干预。在优选方案中，本发明可以以超像素为单位来分割图像的前景区域和背景区域，能有效减少计算量，并且本发明可以以像素为单位提取建成区，和以更大区域为单位提取建成区的方法相比具有更高的准确性。

本发明的其它特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1示出了根据本发明的优选实施方式从图像中提取建成区的流程图；

图2示出了一幅高分辨率遥感图像；

图3示出了对应于图2的边缘图像；

图4是在图像中设置对应于像素的窗口的示意图；

图5示出了对应于图2的基于边缘密度的像素显著性图；

图6示出了对图2进行超像素分割后得到的结果示意图；

图7示出了以超像素为单位将图2所示的遥感图像划分为前景区域和背景区域后的结果示意图；

图8示出了图2所示的图像中的像素为建成区的概率示意图；

图9示出了可施加至图2的二值掩码图像；以及

图10示出了对图2所示的遥感图像进行建成区提取后得到的结果示意图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

图1示出了根据本发明的优选实施方式从遥感图像中提取建成区的流程图。下文中将以图2所示的高分辨率遥感图像作为应用对象对本发明的优选实施方式进行详细描述。

在步骤S11中，可获取图像中每个像素的显著性值。本发明可选用经过图像融合后的高分辨率(例如，每个像素点所表示的区域的长和宽在1～4米范围内)多光谱遥感图像。在遥感图像中，建成区中人工建筑物与其周围光谱差异显著，使得建成区通常具有丰富的边缘，而如耕地、林地、水体等背景区域中边缘较少，因此可基于边缘密度得到像素的显著性值。本优选实施方式中计算显著性值的具体步骤如下：

(1)计算每个像素的梯度幅值。可先将例如遥感图像的彩色图像转换为灰度图像，然后计算每个灰度像素的梯度值，进一步得到像素的梯度幅值。原图像和灰度图像中相同位置上的像素一一对应。I为用像素值表示的灰度图像，可采用下式(1)对I滤波以得到X轴方向和Y轴方向的梯度图像g_x和g_y：

g_x＝I*f₁,g_y＝I*f₂ 式(1)

其中“*”表示卷积运算，f₁和f₂为如下所示的梯度算子：

f₁＝[-1,0,1],f₂＝[-1,0,1]^T 式(2)

将图像g_x和g_y相加可得到例如遥感图像的图像中每个像素的梯度幅值。本发明可采用梯度一阶模计算梯度幅值：

g＝|g_x|+|g_y| 式(3)

根据需要，也可选用梯度二阶模等来计算梯度幅值。

(2)可采用自适应阈值选取方法计算该图像的梯度幅值阈值，以判断每个像素是否为边缘像素。本优选实施方式中，可使用大津法(Otsu法)计算梯度幅值阈值。文献OtsuN.A threshold selection method from gray-level histograms[J].Automatica,1975,11(285-296):23-27中详细介绍了大津法的具体计算步骤，本发明不再对其赘述。针对图像中的每个像素，如果该像素的梯度幅值大于计算出的梯度幅值阈值，则判断该像素为边缘像素；否则判断该像素不为边缘像素。图3示出了根据上述方法得到的对应于图2的边缘图像。图3中边缘像素用白色像素点来表示，非边缘像素用黑色像素点来表示。

除上述方法外，还可使用canny边缘像素检测算法等本领域已知技术手段来得到如图2所示的图像的边缘图像。

(3)针对图像中任意像素l，可在图像中设置包含像素l的大小为w的窗口，用该窗口的边缘像素密度和边缘像素分布来表示像素l的显著性。w为经验值，可根据需要设置。如图4所示，本实施方式中，可在设置窗口时将像素l作为中心，并可用像素l作为坐标原点建立直角坐标系。用n_e表示落在该窗口内的边缘像素数目，用n_w表示该窗口内的像素总数，则该窗口的边缘像素密度可表示为：

式(4)

用n_i,i＝1,2,3,4表示该窗口中的分别落在上述坐标系的四个象限内的边缘像素的数目，则该窗口的边缘像素分布为：

式(5)

本优选实施方式中，像素l的显著性值p(l)可被表示为：

p(l)＝density(l)×evness(l) 式(6)

(4)对图像中每个像素l的显著性值p(l)归一化，得到像素l的基于边缘密度的先验概率p(l^*)：

式(7)

其中，L表示图2所示的图像中的所有像素的集合。

图5示出了对应于图2的基于边缘密度的像素显著性图。图5中亮度越高的像素点的显著性值越大，反之亮度越低的像素点的显著性值越小。

在步骤S12中，可根据显著性值，将例如遥感图像的图像划分为前景区域和背景区域。例如，可使用自适应阈值选取方法计算图像中的像素的显著性值阈值，然后针对图像中的每个像素，如果该像素的显著性值大于计算出的像素的显著性值阈值，则判断该像素属于前景区域，否则判断该像素属于背景区域。

根据本发明的优选实施方式，还可先对遥感图像进行超像素分割，然后以超像素为单位来区分前景区域和背景区域。该优选实施方式的具体步骤如下：

(1)对图2所示的遥感图像进行超像素分割。超像素分割能将图像分割成同质且不重叠的图像块，以超像素而不是像素为处理单元进行后续的前景、背景区域划分能大大减少处理时间和计算量，加快处理速度。可采用本领域已知的任意技术手段进行超像素分割。本优选实施方式可使用SLIC(simple linear iterative clustering，简单线性迭代聚类)超像素分割方法，这是一种基于聚类算法的超像素分割，可由LAB颜色空间以及x、y像素坐标共5维组成特征空间，然后在设定聚类中心后用kmeans聚类方法生成各个超像素。文献Achanta R,Shaji A,Smith K,et al.SLIC superpixels compared to state-of-the-artsuperpixel methods[J].Pattern Analysis and Machine Intelligence,IEEETransactions on,2012,34(11):2274-2282中对SLIC超像素分割方法进行了详细介绍，本发明不再对此赘述。图6示出了对图2进行超像素分割后得到的结果示意图。

(2)计算分割得到的每个超像素的显著性值。设某一超像素s共包含N个像素，该N个像素的集合为s，则超像素s的显著性值可为其所包含的N个像素的显著性值的平均值。本优选实施方式中，将超像素s的显著性值p(s)设为其所包含的N个像素的显著性值的简单算术平均值：

式(8)

可按照式(8)计算图像中每个超像素的显著性值。

(3)判断每个超像素属于前景区域或者背景区域。本优选实施方式中，可先使用自适应阈值选取方法(例如大津法)计算图像的超像素的显著性值阈值；然后针对图像中的每一个超像素，如果该超像素的显著性值大于计算出的超像素的显著性值阈值，则判断该超像素属于前景区域，否则判断该超像素属于背景区域。图7示出了以超像素为单位将图2所示的遥感图像划分为前景区域和背景区域后的结果示意图，图中白色部分表示分割后得到的前景区域，黑色部分表示背景区域。

在步骤S13中，计算像素特征在分割得到的前景区域和背景区域中的特征似然函数。为了得到更为准确的预测结果，可针对多个统计独立的像素特征分别计算其似然函数，以提高提取的准确性。下文中将详细介绍本优选实施方式中所关注的三个特征颜色、方向、方向熵的特征似然函数。这三个特征似然函数可被视为统计独立的似然函数。

(1)颜色特征似然函数。图2所示的遥感图像包括红、绿、蓝和近红外四个波段，本优选实施方式中选择红(R)、绿(G)、蓝(B)三个通道提取颜色特征。可分别计算上述三个分量在前景区域和背景区域中的直方图。假设每个颜色分量用0～255，则可建立6个长度为50、组距为5.1的直方图。将直方图归一化，即可得到对应的似然函数。用p(R|f)、p(G|f)、p(B|f)分别表示R、G、B三个分量在前景区域中的似然函数，用p(R|b)、p(G|b)、p(B|b)分别表示R、G、B三个分量在背景区域中的似然函数，其中f表示前景(foreground)，b表示背景(background)，则颜色特征在前景区域和背景区域的似然函数p(c|f)和p(c|b)分别为：

p(c|f)＝p(R|f)p(G|f)p(B|f) 式(9)

p(c|b)＝p(R|b)p(G|b)p(B|b)

(2)方向特征似然函数。在图像中，建筑区域方向主要是建筑物轮廓线段所在的方向，而非建筑区域中方向杂乱无规则，因此方向是一种比较强的区分特征。可按照下式计算任意像素l的方向值θ(l)：

式(10)

其中，g_y(l)表示像素l的Y轴方向的梯度值，g_x(l)表示像素l的X轴方向的梯度值，g_x和g_y可参见式(1)。计算出所有像素的方向值后，可分别统计该方向值在前景区域和背景区域的似然函数p(θ|f)和p(θ|b)。例如，可将式(10)计算出的方向值映射到[0°,180°]，然后建立该方向值在前景区域和背景区域中的直方图(例如组距为18°、长度为10的直方图)，然后通过归一化得到p(θ|f)和p(θ|b)。

(3)方向熵特征似然函数。可以采用本领域已知的任意方法计算图像中的像素的方向熵特征似然函数。例如，优选地，如果图像被分割为多个超像素，可以根据像素l所在的超像素s内的各个像素的方向值，计算像素l的方向熵。设像素l位于超像素s中，可根据式(10)计算超像素s中每个像素的方向值，并建立超像素s中的像素的方向值的直方图，进一步得到超像素s内的方向值分布函数p_s(θ)。设超像素s内的所有像素的方向值的集合为D，则位于超像素s内的任意像素l的方向熵H(l)为：

式(11)

可分别计算前景区域和背景区域中的方向熵似然函数p(H|f)和p(H|b)。例如，可分别建立前景区域和背景区域中的方向熵直方图，然后通过归一化方向熵直方图得到对应的似然函数。

步骤S14，可计算例如遥感图像的图像中每个像素为建成区的概率值。对于图像中的任意像素l，根据贝叶斯公式可得像素l为建成区的概率值p(buildings|l)：

式(12)

将式(12)展开得到：

式(13)

其中，p(f)表示像素l属于前景区域的先验概率，此处为像素l的显著性值的归一化值，即p(f)＝p(l^*)，可参见式(1)；p(b)表示像素l属于背景区域的概率，p(b)＝1-p(f)；p(ch|f)为前景区域的像素特征似然函数，p(ch|b)为背景区域的像素特征似然函数，p(c|f)、p(θ|f)、p(H|f)、p(c|b)、p(θ|b)、p(H|b)为上述步骤S15中得到的像素特征的似然函数，因为颜色、方向和方向熵可被视为统计独立的像素特征，所以有p(ch|f)＝p(c|f)p(θ|f)p(H|f)和p(ch|b)＝p(c|b)p(θ|b)p(H|b)。可应用式(13)计算图像中每个像素为建成区的概率值。图8示出了图2所示的图像中的像素为建成区的概率示意图，图8中亮度越高的像素点为建成区的概率越高，反之，亮度越低的像素点为建成区的概率越低。

在步骤S15中，判断图像中的每个像素是否为建成区。可使用例如大津法的自适应阈值选取方法计算在图像中的像素为建成区的概率阈值，然后针对每个像素，如果该像素为建成区的概率值大于概率阈值，则判断该像素为建成区，否则判断该像素不为建成区，得到的二值掩码图像如图9所示。图9中，白色的像素点表示该点被判断为建成区，黑色的像素点表示该点被判断不为建成区。

在步骤S16中，从被判断为建成区的区域中提取建成区。将上述得到的二值掩码图像应用于图2所示的遥感图像，从对应于被判断为建成区的区域中提取建成区，提取结果如图10所示。可使用本领域已知的任意技术手段，如支持向量机等，从被判断为建成区的区域中提取建成区。从而，可快速、准确、全自动地从例如遥感图像的图像中提取建成区。

以上结合附图详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (9)

1.一种从遥感图像中提取建成区的方法，包括：

获取图像中每个像素的显著性值；

根据所述显著性值将所述图像划分为前景区域和背景区域；

分别计算所述前景区域和所述背景区域的像素特征似然函数；

针对所述图像中的每个像素，根据该像素的所述显著性值和所述像素特征似然函数得到该像素为建筑区的概率值；

根据所述概率值，判断所述图像中的每个像素是否为建成区；以及

从所述图像中对应于被判断为建成区的区域中提取建成区，

其中，获取图像中每个元素的显著性值包括：

判断所述图像中的每个像素是否为边缘像素；

针对所述图像中的每个像素，在所述图像中设置对应于该像素的窗口，该像素位于该窗口中，并计算该窗口的边缘像素密度和边缘像素分布；以及

所述图像中每个像素的所述显著性值等于对应于该像素的所述窗口的所述边缘像素密度和所述边缘像素分布的乘积。

2.根据权利要求1所述的从遥感图像中提取建成区的方法，其中，所述窗口的所述边缘像素密度为其中，n_e表示该窗口内的所述边缘像素的数目，n_w表示该窗口内的像素总数。

3.根据权利要求1所述的从遥感图像中提取建成区的方法，其中，每个像素位于对应于该像素的所述窗口的中心，以该元素为坐标原点建立直角坐标系，用n_i,i＝1,2,3,4表示对应于该像素的所述窗口中分别落在该坐标系的四个象限内的所述边缘像素的数目，则该窗口的所述边缘像素分布为

4.根据权利要求1所述的从遥感图像中提取建成区的方法，其中，根据所述显著性值将所述图像划分为前景区域和背景区域包括：

对所述图像进行超像素分割，并且计算每个超像素的显著性值；以及

根据所述超像素的所述显著性值将所述图像以超像素为单位划分为所述前景区域和所述背景区域。

5.根据权利要求4所述的从遥感图像中提取建成区的方法，其中，每个所述超像素的显著性值为包含在该超像素中的所有像素的所述显著性值的平均值。

6.根据权利要求4所述的从遥感图像中提取建成区的方法，其中，将所述图像以超像素为单位划分为所述前景区域和所述背景区域包括：

采用自适应阈值选取方法得到所述图像中的超像素的显著性值阈值，针对所述图像中的每个超像素，如果该超像素的所述显著性值大于超像素的所述显著性值阈值，则判断该超像素属于所述前景区域，否则判断该超像素属于所述背景区域。

7.根据权利要求1所述的从遥感图像中提取建成区的方法，其中，所述像素特征似然函数包括颜色特征似然函数、方向特征似然函数和方向熵特征似然函数中的至少一者。

8.根据权利要求1或7所述的从遥感图像中提取建成区的方法，其中，基于贝叶斯公式根据下式计算所述图像中的每个像素为建成区的概率值p(bulidings|l)：

$p\left(bulidings\right|l)=\frac{p\left(f\right)\×p\left(ch\right|f)}{p\left(f\right)\×p\left(ch\right|f)+p\left(b\right)*p\left(ch\right|b)},$

其中，用l表示所述图像中的任意一个像素，p(f)表示该像素l属于所述前景区域的概率，此处p(f)为该像素l的所述显著性值归一化后的结果；p(b)表示该像素l属于所述背景区域的概率，p(b)＝1-p(f)；p(ch|f)为所述前景区域的所述像素特征似然函数，p(ch|b)为所述背景区域的所述像素特征似然函数。

9.根据权利要求1所述的从遥感图像中提取建成区的方法，其中，判断所述图像中的每个像素是否为建成区包括：

采用自适应阈值选取方法得到所述图像中的像素为建成区的概率阈值，针对所述图像中的每个像素，如果该像素的所述概率值大于所述概率阈值，则判断该像素为建成区，否则判断该像素不为建成区。