CN115880575B - 结合变化信息和建筑特征的遥感影像新增建筑提取方法 - Google Patents

Abstract

本发明公开了一种结合变化信息和建筑特征的遥感影像新增建筑提取方法，属于遥感图像处理领域。其包括：分割后时相遥感影像获得超像素图；提取线段得到线段集合；计算超像素的边缘线段垂直性；计算建筑物线段指数；提取影像基于对象的形态学特征；计算建筑物强度；计算前后时相遥感影像的纹理特征；计算前后时相遥感影像的纹理特征差和光谱特征差；求纹理和光谱特征的变化强度，获得变化强度；获取新增建筑对象斑块；采用形状指数进行后处理，获取最终的新增建筑对象。本发明充分利用光谱和纹理特征提取变化信息，并通过边缘结构特征和形态学特征提取建筑特性，为大面积遥感影像实现新增建筑识别提供了可靠的技术支撑。

Description

结合变化信息和建筑特征的遥感影像新增建筑提取方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种结合变化信息和建筑特征的遥感影像新增建筑提取方法。

背景技术

随着国民经济的发展和人口增加，我国工业化和城市化水平逐步提高，建设用地不断增加，耕地急剧减少，因此，及时掌握建设用地、耕地的更新状况，尤其是违法建筑物的情况，成为保障社会经济协调持续发展的重要任务。近年来，高分辨率遥感发展迅速，由于其可以提供地物更加详细的光谱和纹理信息，为建筑物等小尺度目标的变化检测提供了可能。传统的变化检测方法应用于高分辨率遥感影像中新增建筑物检测面临以下两点问题：一方面，基于像素级光谱信息的变化检测不能够很好的获取高分影像中地物的纹理和几何信息，提取精度较低；另一方面，高分影像中丰富的光谱和纹理信息为提取建筑物带来一定的困难，表现为其光谱与道路、裸土等地物相似性较高，结果中噪声较多。

虽然许多学者对高分辨率影像中建筑提取做了研究，但是这些研究使用的特征较为单一，无法很好体现建筑物特征及其变化信息；其次，一部分研究采用监督分类的方法，方法中面临样本难选、人工干预程度高的问题，给大面积、自动化、快速的提取变化建筑带来困难。

发明内容

本发明的目的在于提供一种结合变化信息和建筑特征的遥感影像新增建筑提取方法，其充分利用光谱和纹理特征提取变化信息，并通过边缘结构特征和形态学特征提取建筑特性，为大面积遥感影像实现新增建筑识别提供可靠的技术支撑。

本发明采用的技术方案为：

一种结合变化信息和建筑特征的遥感影像新增建筑提取方法，包括以下步骤：

步骤1，采用SEEDS算法将后时相遥感影像

分割，获得超像素图IM_seg；

步骤2，采用LSD算法对后时相遥感影像

提取线段，得到线段集合L_set；

步骤3，基于线段集合L_set计算超像素S_i的边缘线段垂直性

步骤4，利用边缘线段垂直性和相交线段数计算建筑物线段指数

步骤5，利用形态学建筑指数MBI提取影像

基于对象的形态学建筑特征

步骤6，融合形态学建筑特征

和建筑线段指数

计算影像

的建筑物强度

步骤7，分别计算前后时相遥感影像的纹理特征

和

步骤8，分别计算前后时相遥感影像的纹理特征差DT和光谱特征差DS；

步骤9，分别求纹理和光谱特征的变化强度IT和IS，通过直接求和获得变化强度CI；

步骤10，依据分割结果IM_seg求对象级的变化强度CI_seg，并与建筑物强度

加权融合得到新增建筑强度图NBI；

步骤11，计算新增建筑强度图NBI的分割阈值T，对NBI二值化后获取新增建筑对象斑块NB；

步骤12，对每一新增建筑对象斑块NB分别计算形状指数GI，过滤出GI超过阈值的斑块NB，得到最终提取的新增建筑对象NB_post。

进一步地，步骤3中，当有n条线段与超像素S_i相交时，边缘线段垂直性

的表达式为：

式中，F(Δβ_jk)表示当两条线段l_j、l_k与超像素S_i相交时的垂直性值，Δβ_jk表示线段之间的夹角，θ表示倾斜角差容许误差。

进一步地，步骤4中建筑物线段指数BLI_Si采用如下公式计算：

其中，N()表示线段包含的像素个数，N(S_i∩l_m)表示线段l_m与超像素S_i相交的像素个数。

进一步地，步骤5中，基于对象的形态学建筑特征

计算方式如下：

式中，MBI(x,y)表示形态学建筑指数MBI图的(x,y)处的象元值，p和q分别表示与超像素S_i等面积的矩形的长和宽。

进一步地，步骤6中建筑物强度

的计算方法如下：

式中，

为可调节的加权因子，Nor()表示对建筑指数进行最大值最小值归一化处理。

进一步地，步骤7中，采用GLCM的方差测度作为纹理特征，其计算方法如下：

其中，d∈{1,2}，P(i,j)代表灰度级别i和j在给定方向和距离情况下的出现的可能性，μ是GLCM的均值，N代表灰度级。

进一步地，步骤8中，纹理特征差DT和光谱特征差DS的计算方法如下：

式中，S表示光谱波段，T₁和T₂分别表示遥感影像的两个时间，abs()表示求绝对值。

进一步地，步骤9中，纹理和光谱特征的变化强度IT和IS的计算方式如下：

其中，上标T表示转置运算，X_DT表示纹理特征差DT的特征矩阵，行为像素数，列为波段数，X_DS表示纹理特征差DS的特征矩阵，行为像素数，列为波段数；

变化强度CI采用如下方式计算：

CI＝Nor(IT)+Nor(IS)

其中，Nor()表示对变化特征进行最大值最小值归一化处理。

进一步地，步骤10中，依据分割结果IM_seg求对象级的变化强度CI_seg的计算方式如下：

式中，CI(x,y)表示CI图中(x,y)处的象元值，p_c和q_c分别表示与超像素S_i等面积的矩形的长和宽，

表示超像素S_i的对象级变化强度；

新增建筑强度图NBI的计算方式如下：

其中，Nor()表示对变化信息和建筑特征进行最大值最小值归一化处理。

进一步地，步骤11中，分割阈值T的计算方式如下：

T＝ms+a×std

其中，ms和std分别表示统计NBI中所有对象的均值和标准差，a为调节因子；

对NBI二值化后获取新增建筑对象斑块NB，具体方式如下：

式中，T是基于所有对象隶属度确定的阈值，NB表示新增建筑物提取结果，像元灰度值为1表示新增建筑物，像元值为0表示背景值；

步骤12中，形状指数GI的计算方法为：

其中，GI表示形状指数，是基于连通域计算的特征，P表示该连通域的周长，A表示连通域的面积；

最终提取的新增建筑对象NB_post的计算方式如下：

T_GI＝a×std(GI)+mean(GI)

其中，GI_j表示针对第j个斑块NB计算的形状指数，T_GI表示形状指数阈值，像元灰度值为1表示新增建筑物，像元值为0表示背景值；a表示调节系数，std()表示求图像灰度值的均方差，mean()表示求图像灰度值的均值。

本发明具有如下有益效果：

(1)本发明提出了基于变化强度和后时相建筑特征的新增建筑提取方法，在不需要训练样本的情况下，可以有效提取多时相遥感图像中的新增建筑。

(2)应用本发明方法可以有效提取融合纹理和光谱信息的变化信息，为其他土地类型的变化检测提供关键的技术支撑。

总之，本发明采用了全新的非监督识别方法，结合了光谱、纹理、光谱空间和边缘结构特征提取建筑特征及其变化信息，有效地实现了像素级变化信息和对象级目标信息的融合，为自动高效的新增建筑提取提供重要的技术途径。

附图说明

图1为新增建筑提取方法的原理示意图。

具体实施方式

下面结合附图及具体实施方式对本发明做进一步详细的说明。

一种结合变化信息和建筑特征的遥感影像新增建筑提取方法，包括以下步骤：

步骤1，采用SEEDS算法将后时相遥感影像

分割，获得超像素图IM_seg；

步骤2，采用LSD算法对后时相遥感影像

提取线段，得到线段集合L_set，由n条线段组成的集合可表示为：

L_set＝{l₁,l₂,l₃...l_m,m＝1,...,n}

步骤3，基于线段集合L_set计算超像素S_i的边缘线段垂直性

当与超像素S_i有n条线段相交时，边缘线段垂直性

计算方式如下：

式中，F(Δβ_jk)表示当两条线段l_j、l_k与超像素S_i相交时其垂直性值，Δβ_jk表示线段之间的夹角，θ表示倾斜角差容许误差。

步骤4，利用线段垂直性和相交线段数计算建筑物线段指数

计算方法如下：

其中，N()表示线段包含的像素个数，N(S_i∩l_m)表示l_m与超像素S_i相交的像素个数。

步骤5，利用形态学建筑指数(MBI)提取影像

基于对象的形态学特征

计算表达式如下：

式中，MBI(x,y)表示MBI图的(x,y)处的象元值，p和q分别表示与超像素S_i等面积的矩形的长和宽。

步骤6，融合形态学特征

和建筑线段指数

计算影像

的建筑物强度

其计算方法如下：

式中，

为可调节的加权因子，Nor()表示对建筑指数归一化处理。

步骤7，采用GLCM的方差测度作为纹理特征，分别计算前后时相遥感影像的纹理特征

和

其计算方法如下：

其中，d∈{1,2}，P(i,j)代表灰度级别i和j在给定方向和距离情况下的出现的可能性，μ是GLCM的均值，N代表灰度级。

步骤8，分别计算前后时相遥感影像的纹理特征差DT和光谱特征差DS，方法如下：

式中，S表示光谱波段，abs()表示求加绝对值；

步骤9，分别求纹理和光谱特征的变化强度IT和IS，其计算方式如下：

其中，*^T表示转置运算，X_DT表示纹理特征差DT的特征矩阵，行为像素数，列为波段数，X_DS表示纹理特征差DS的特征矩阵，行为像素数，列为波段数；

通过对归一化后的IT和IS直接求和，获得变化强度CI，

CI＝Nor(IT)+Nor(IS)

其中，Nor()表示对变化特征归一化处理。

步骤10，依据分割结果IM_seg求对象级的变化强度CI_seg的计算方式如下：

式中，CI(x,y)表示CI图中(x,y)处的象元值，p_c和q_c分别表示与超像素S_i等面积的矩形的长和宽。

变化强度CI_seg与建筑物强度

加权融合得到新增建筑强度NBI：

其中，Nor()表示对变化信息和建筑特征进行最大值最小值归一化处理；

步骤11，计算新增建筑强度图NBI的分割阈值T，其计算方式如下：

T＝ms+a×std

其中，ms和std分别统计NBI中所有对象的均值和标准差，a为调节因子；

然后对NBI二值化后获取新增建筑对象斑块NB，判断条件如下：

式中，T是基于所有对象隶属度确定的阈值，NB表示最终的新增建筑物提取结果，像元灰度值为1表示新增建筑物，像元值为0表示背景值。

步骤12，采用形状指数GI进行后处理，可以过滤道路等长宽比较大的斑块，形状指数GI计算方式如下：

其中，GI表示形状指数，是基于连通域计算的特征，P表示该连通域的周长，A表示连通域的面积；

最后提取的新增建筑对象NB_post，计算方式如下：

其中，GI_j表示第j个连通域计算的形状指数，T_GI表示形状指数阈值。

总之，本发明充分利用光谱和纹理特征提取变化信息，并通过边缘结构特征和形态学特征提取建筑特性，可以有效提取高分辨率遥感影像中的新增建筑物。该方法无需选择训练样本，无需人工干预识别，为便捷、快速、精确和大范围的新增建筑物提取提供了重要的技术支撑。

Claims (10)

1.一种结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，包括以下步骤：

步骤1，采用SEEDS算法将后时相遥感影像

分割，获得超像素图IM_seg；

步骤2，采用LSD算法对后时相遥感影像

提取线段，得到线段集合L_set；

步骤3，基于线段集合L_set计算超像素S_i的边缘线段垂直性

步骤4，利用边缘线段垂直性和相交线段数计算建筑物线段指数

步骤5，利用形态学建筑指数MBI提取影像

基于对象的形态学建筑特征

步骤6，融合形态学建筑特征

和建筑线段指数

计算影像

的建筑物强度

步骤7，分别计算前后时相遥感影像的纹理特征

和

步骤8，分别计算前后时相遥感影像的纹理特征差DT和光谱特征差DS；

步骤9，分别求纹理和光谱特征的变化强度IT和IS，通过直接求和获得变化强度CI；

步骤10，依据分割结果IM_seg求对象级的变化强度CI_seg，并与建筑物强度

加权融合得到新增建筑强度图NBI；

步骤11，计算新增建筑强度图NBI的分割阈值T，对NBI二值化后获取新增建筑对象斑块NB；

步骤12，对每一新增建筑对象斑块NB分别计算形状指数GI，过滤出GI超过阈值的斑块NB，得到最终提取的新增建筑对象NB_post。

2.根据权利要求1所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤3中，当有n条线段与超像素S_i相交时，边缘线段垂直性

的表达式为：

式中，F(Δβ_jk)表示当两条线段l_j、l_k与超像素S_i相交时的垂直性值，Δβ_jk表示线段之间的夹角，θ表示倾斜角差容许误差。

3.根据权利要求2所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤4中建筑物线段指数

采用如下公式计算：

其中，N()表示线段包含的像素个数，N(S_i∩l_m)表示线段l_m与超像素S_i相交的像素个数。

4.根据权利要求3所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤5中，基于对象的形态学建筑特征

计算方式如下：

式中，MBI(x,y)表示形态学建筑指数MBI图的(x,y)处的象元值，p和q分别表示与超像素S_i等面积的矩形的长和宽。

5.根据权利要求4所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤6中建筑物强度

的计算方法如下：

式中，

为可调节的加权因子，Nor()表示对建筑指数进行最大值最小值归一化处理。

6.根据权利要求5所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤7中，采用GLCM的方差测度作为纹理特征，其计算方法如下：

其中，d∈{1,2}，P(i,j)代表灰度级别i和j在给定方向和距离情况下的出现的可能性，μ是GLCM的均值，N代表灰度级。

7.根据权利要求6所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤8中，纹理特征差DT和光谱特征差DS的计算方法如下：

式中，S表示光谱波段，T₁和T₂分别表示遥感影像的两个时间，abs()表示求绝对值。

8.根据权利要求7所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤9中，纹理和光谱特征的变化强度IT和IS的计算方式如下：

其中，上标T表示转置运算，X_DT表示纹理特征差DT的特征矩阵，行为像素数，列为波段数，X_DS表示纹理特征差DS的特征矩阵，行为像素数，列为波段数；

变化强度CI采用如下方式计算：

CI＝Nor(IT)+Nor(IS)

其中，Nor()表示对变化特征进行最大值最小值归一化处理。

9.根据权利要求8所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤10中，依据分割结果IM_seg求对象级的变化强度CI_seg的计算方式如下：

式中，CI(x,y)表示CI图中(x,y)处的象元值，p_c和q_c分别表示与超像素S_i等面积的矩形的长和宽，

表示超像素S_i的对象级变化强度；

新增建筑强度图NBI的计算方式如下：

其中，Nor()表示对变化信息和建筑特征进行最大值最小值归一化处理。

10.根据权利要求9所述的结合变化信息和建筑特征的遥感影像新增建筑提取方法，其特征在于，步骤11中，分割阈值T的计算方式如下：

T＝ms+a×std

其中，ms和std分别表示统计NBI中所有对象的均值和标准差，a为调节因子；

对NBI二值化后获取新增建筑对象斑块NB，具体方式如下：

式中，T是基于所有对象隶属度确定的阈值，NB表示新增建筑物提取结果，像元灰度值为1表示新增建筑物，像元值为0表示背景值；

步骤12中，形状指数GI的计算方法为：

其中，GI表示形状指数，是基于连通域计算的特征，P表示该连通域的周长，A表示连通域的面积；

最终提取的新增建筑对象NB_post的计算方式如下：

T_GI＝a×std(GI)+mean(GI)

其中，GI_j表示针对第j个斑块NB计算的形状指数，T_GI表示形状指数阈值，像元灰度值为1表示新增建筑物，像元值为0表示背景值；a表示调节系数，std()表示求图像灰度值的均方差，mean()表示求图像灰度值的均值。

Images