CN111339948A - 一种高分辨率遥感影像新增建筑自动识别方法 - Google Patents

Abstract

本发明涉及一种高分辨率遥感影像新增建筑自动识别方法，主要步骤包括：1)对不同时相的高分辨率影像进行几何精纠正和相对辐射校正；2)计算获取前后两个时相影像的MBI(形态学建筑指数)、Harris、NDVI(归一化植被指数)特征；3)建筑特征中去除植被区域；4)计算特征热点图；5)特征热点图取阈值；6)对二值化后的特征图取连通域得到对象，采用对象投票处理得到新增建筑斑块；7)新增建筑后处理，包括：前期建设用地现状图和最小制图单元滤除。本发明能基于遥感图像自动提取新增建筑，为土地监察和管理提供有效的技术手段。

Description

一种高分辨率遥感影像新增建筑自动识别方法

技术领域

本发明属于遥感影像处理技术领域，主要涉及土地利用变化检测中的新增建筑变化检测方法，具体涉及一种基于MBI(形态学建筑指数)和Harris的高分辨率遥感影像新增建筑自动识别方法。

背景技术

由于建设用地扩展需占用其它非建筑区域(如：耕地、林地、裸地等)，土地利用变化频繁，我国原有的依靠大量人力目视解译影像为主的土地调查方法已不能满足新形势下对土地管理的要求。迅速发展的遥感技术和计算机技术已成为进行土地利用监测，发现违章建筑的重要手段。随着国土资源管理工作的深入，要求及时发现建筑、重点关注专题用地的变化，因此对高效率、短时间、大范围的自动获取土地利用动态变化提出了迫切要求，依靠大量人力结合遥感影像技术的土地检测方案已不能满足需求。所以自动发现新增建筑，结合少量的人工干预，实现高效精确的土地动态监测，是提高作业效率和实现大规模监测的重要技术途径，为切实保护耕地、加强国土资源管理提供了技术保障。

近年来，由于高分辨率遥感影像可以提供更加丰富的地物细节信息，许多学者使用高分辨率遥感影像来提取建筑物等人工地物的变化。相比于中低分辨率遥感影像，高分辨率遥感影可以提供更加丰富的信息，同时也给传统的遥感分析和目标变化检测带来了挑战。现有的建筑物变化检测思路大致可以分为两类：1)根据建筑物的特点，使用几何、纹理和阴影等特征直接获取建筑物的变化信息。有的学者综合利用人工建筑的空间特征和纹理特征信息，提出基于特征级与像素级的变化检测方法，用以提取城区建筑物变化；有的学者利用背景模型的目标检测方法对阴影进行检测，然后通过阴影补偿的方法来检测建筑物的变化；有的学者采用特征点的相关运算，对人工目标的变化进行检测；还有一些学者根据特定的建筑指数，因其对高分辨率遥感影像建筑有较好的响应，利用特征层和决策层进行融合提取建筑物变化信息。2)首先检测出变化区域，然后再根据监督分类方式识别出变化建筑。有的学者为了调查了两个时相城市建筑的更新情况，前时相参考数据为辅助地图矢量数据，后时相是采用面向对象分类方法得到建筑，两个时相结果对比得到新增建筑；还有的学者采用“from-to”的监督分类方式，标记由其余类别变为建筑物的变化类别标签，然后输入分类器再由监督分类得到变化建筑物。但是，现有方法使用的特征较为单一，不能完全的表达建筑物特点，识别建筑物的完整性和精度较低，噪声也较多；其次，现有的监督分类方法识别建筑变化仍需要选择训练样本，人工干预较多，自动化程度低。

依据国情监测和土地管理的需要，利用遥感数据快速、准确、高效的识别新增建筑，结合少量的人工干预，已成为研究的热点。本技术方案主要依据高分辨率遥感影像，采用多特征组合的方式，不需要人工选择样本的情况下实现新增建筑自动提取。

发明内容

本发明考虑高分辨率遥感影像上房屋具有高亮的表面和形状较为规则的特点，提供了一种基于MBI和Harris的高分辨率遥感影像新增建筑自动提取方法。

本发明采用的技术方案如下：

步骤1、对不同时相的高分辨率影像进行几何精纠正和相对辐射校正；

步骤2、获取前后两个时相影像的特征：形态学建筑指数MBI、角点特征Harris、归一化差分植被指数NDVI；

步骤3、建筑特征中去除植被区域：

分别对两个时相的NDVI图进行二值化处理得到前后时相植被覆盖区域，然后分别使用得到的植被区域过滤两个时相相应的建筑特征MBI和Harris，将植被区域赋值为背景值；

步骤4、计算特征热点图：

两个时相对应建筑特征值分别求差，然后对求差后的特征变化图进行半重叠开窗计算求窗口内平均亮度值，最后得到两个特征的热点图：IMBI和IHarris；

步骤5、特征热点图取阈值：

对两个特征热点取阈值得到其对应二值图，即每个特征确定的新增建筑；

步骤6、对二值化后的特征图取连通域得到对象，采用对象投票处理得到新增建筑斑块；

步骤7、新增建筑后处理，包括基于前期建设用地现状图和最小制图单元滤除。

进一步的，步骤2中形态学建筑指数MBI是通过建立建筑的光谱和空间特征与形态学算子之间的关系，实现建筑物的提取，具体的计算步骤如下：

2.1.1，获取遥感影像上各像元在可见光波段中光谱值的最大值，作为该像元的亮度值，从而获得亮度影像：

式(1)中，M_k(x)为像元x在波段k中的光谱值，b(x)则代表像元x的亮度值，K为总波段数；

2.1.2，利用差分形态学属性的白顶帽值表示建筑的光谱和结构特征，基于多尺度和多方向线性结构元素的差分形态学属性白顶帽值定义为：

TH-DMP(d，s)＝|TH_b(d，s)-TH_b(d，s-Δs)| (2)

其中

代表基于重建的开运算后的亮度影像b，s和d分别表示线性结构元素SE的尺度和方向，Δs是尺度的间隔；

2.1.3，计算形态学房屋指数MBI，即不同尺度和方向的基于差分形态学属性白顶帽值的均值：

其中D和S分别是方向和尺度的数目。

进一步的，步骤2中角点特征Harris是一种基于信号的角点检测算法，其主要思想是使用一个检测窗口，当该窗口沿各个方向作微小移动时，检测窗口内的平均能量是否发生变化，当该能量变化值超过设定的阈值时，就将窗口的中心像素点提取为角点，具体计算过程如下，

2.2.1，计算亮度影像：

M_k(x)为像元x在波段k中的光谱值，b(x)则代表像元x的亮度值，K为总波段数；

2.2.2，对亮度影像b分别求水平和垂直方向上的梯度：

式中，X和Y分别表示对亮度影像b求水平和垂直方向上的梯度，梯度滤波器f＝(-2，-1，0，1，2)；

2.2.3，对梯度积高斯滤波，

式中，w表示高斯滤波器，滤波窗口内的元素权重不同，中心位置高，边缘位置低；X²、Y²和XY梯度运算均是基于对应的梯度像素进行的；

2.2.4，建立度量函数，

R＝AB-C²-k′(A+B)² (10)

式中的k′是一个系数，AB和C²表示求积，指整个影像对应位置像素求积，与梯度积类似，而A+B则是表示影像对应位置上的像素求和；R是度量角点特征的响应指标，其值大于某个阈值则认为是角点。

进一步的，步骤2中归一化植被指数NDVI的计算方式如下，

式中的ρ_nir和ρ_red分别表示遥感影像中的近红外波段和红色波段的像元值。

进一步的，步骤3的具体实现方式如下，

首先采用一种自适应的阈值进行二值化处理，两个时相NDVI图均为单波段影像，对像元值进行统计，得到所有有效像元的均值μ和标准差σ，则自适应阈值T＝μ+σ，经过以下处理则得到二值图：

g₁′(x，y)和g₁(x，y)分别代表了二值化后的像元值和原始NDVI值，即NDVI影像中大于阈值的像元为1，表示植被区域，其余的值为0，代表背景值。

NDVI影像取阈值后，得到只含有非植被和植被的二值图，将二值图中像元值为1的区域提取出来，在建筑指数(MBI和Harris)中找到相应的位置，并赋为背景值，即完成了建筑指数中植被剔除：

式中g₁′(x，y)表示NDVI二值图的像元值，h₁(x，y)和h₁(x，y)分别表示建筑指数(MBI和Harris)在处理前后的像元值。

进一步的，步骤4的具体实现方式如下，

4.1，两个时相建筑指数加权求差，

其中，j∈{MBI，Harris}，w^j表示由两个时相的特征计算得到的加权值，

和

分别表示两个时相的特征，D^j表示两个特征加权求差后的影像值；

4.2，半重叠开窗计算得到热点图，

首先，确定窗口大小，设最小检测单元的面积为Area_min，窗口大小m：

式中，res为影像的分辨率，int是对结果求整；

半重叠滤波：

式中，j∈{MBI，Harris}，I^j为半重叠后的加权差值特征图，记为IMBI和IHarris，

表示使用滤波器对全图进行半重叠滤波运算，f_mean为平均滤波算子，由开窗的大小决定，大小为m×m，每个元素值为1/m；

经过以上的处理后，突出了两个时相建筑信息的差异。

进一步的，步骤5的具体实现方式如下，

5.1，IMBI二值化，

先将IMBI影像中最小值和最大值线性拉伸到0-1：

式中，I_max和I_min表示IMBI影像中像素的最大值和最小值，g₂(x，y)和g₂′(x，y)分别表示线性拉伸前后的像素值，对拉伸后的影像取固定的阈值(推荐阈值为0.59)：

NL1(x，y)表示IMBI取阈值后的二值图，该结果就是基于MBI得到的新增建筑。

5.2，IHarris二值化，

IHarris阈值按如下的方法确定：

T_top＝sort_0.98(IH) (20)

T_low＝sort_0.94(IH) (21)

式中，IH为IHarris强度图像中的所有有效的像元值，即非零值，sort_0.98(IH)表示将所有像元按照升序的方式排列，取第98％位置处的值，sort_0.94(IH)则表示取第94％位置处的值，T_h则是最后的Harris阈值，IHarris的二值化方法为：

NL2(x，y)为IHarris取阈值后的二值图，该结果就是基于Harris特征得到的新增建筑斑块，h₂(x，y)是指IHarris特征。

进一步的，步骤6的具体实现方式如下，

6.1，提取连通域

取NL1和NL2的4方向连通域，其中NL1是由MBI确定的新增建筑，NL2是由Harris确定的新增建筑，则有：

其中，i∈[1，N_M]和j∈[1，N_H]，N_M和N_H分别表示NL1和NL2的连通域数目，

和

分别表示利用MBI和Harris特征识别的新增建筑对象斑块，它们是由相邻关系的点构成的对象，Δc和Δr由连通域的大小所决定，c0和r0表示连通域的起点位置；

6.2，对象投票

采用对象投票方式来完成

对

中伪新增建筑的剔除：

式中，Num(*)表示对象的像素数目，

表示两个按照位置求交集，Tp为对象投票比例参数；经过对象投票剔除后，得到新增建筑

i∈[1，N1_M]，N1_M是经过对象投票后的对象数。

进一步的，步骤7的具体实现方式如下，

1)使用前期建设用地现状图过滤结果：

使用前期建设用地现状图过滤掉新增建筑与建设用地重叠的部分，得到更新后的新增建筑斑块，现状图是只含有0和1的二值图，其中，0表示非建设用地，1表示建设用地，过滤方法如下：

式中f₁(x，y)表示前期现状图中的像元值，h₃(x，y)和h₃(x，y)分别表示新增建筑斑块

处理前后的像元值；

经过以上步骤的处理，得到的新增建筑结果记为

i∈[1，N2_M]，N2_M是经过前期建设用地过滤后的对象数；

2)最小制图单元过滤：

此处的过滤采用面向对象的处理方式：

式中，res为影像的分辨率，Area_min为最小检测单元；

经过上述步骤处理后，得到最终的新增建筑

i∈[1，N3_M]，N3_M是经过剔除对象后的对象数。

本发明的主要思路是先提取较为可靠的建筑变化信息，然后根据一系列约束条件筛除伪新增建筑，得到最终的新增建筑。针对高分辨率遥感影像中房屋亮度大、各向同性、类矩度的特点，采用形态学房屋指数自动提取遥感影像房屋信息，可以取得较好的效果；同时，为了避免出现像素级检测的噪声现象和面向对象分割的参数设置的复杂性，采用特征二值化取连通域得到对象，然后采用对象投票的方式确定新增建筑。基于MBI和步骤4中半重叠计算可以较完整的识别新增建筑轮廓，但会有伪检测结果，而角点特征Harris在建筑物上可以产生较强的响应，但是轮廓信息不完整，故本发明采用了对象投票的方式融合Harris特征和MBI的信息，剔除伪新增建筑，得到较为理想的实验结果。

和现有技术相比，本发明具有如下优势：(1)无需设置对象分割参数，无需人工选择训练样本，可实现高分辨率遥感影像新增建筑全自动化提取；(2)整个技术流程直接计算了建筑的变化信息，不涉及其它纹理指数的计算；(3)本发明需要调节的参数较少，计算时间较其他方法较短，自动化程度高，可以广泛的使用于大区域新增建筑提取领域。

附图说明

图1新增建筑提取流程图。

具体实施方式

步骤1、对不同时相的高分辨率影像进行几何精纠正和相对辐射校正；

步骤2、获取前后两个时相影像的特征：MBI、Harris、NDVI；

步骤2中直接通过两个建筑指数提取多时相建筑信息：形态学建筑指数(MBI)、Harris角点特征，Harris是依据亮度图像(可见光波段最大值)计算得到，是一种计算角点特征的指数，使用一个检测窗口，沿各个方向作微小移动时检测窗口内平均能量是否发生变化，是一种比较稳定的特征提取算法；MBI也是由亮度图像计算得到的，是一种自动提取高分辨率影像上建筑的形态学指数，根据影像上建筑的光谱和纹理信息进行形态学运算突出建筑的信息。

MBI(Morphological Building Index)：一种自动提取高分辨率影像上建筑的形态学指数，根据影像上建筑的光谱结构信息进行形态学运算突出建筑的信息，是一种无需选择样本的非监督建筑提取指数，具有较广泛的适用性。

Harris特征：一种计算角点特征的指数，使用一个检测窗口，当该窗口沿各个方向作微小移动时，检测窗口内的平均能量是否发生变化，当该能量变化值超过设定的阈值时，就将窗口的中心像素点提取为角点。

NDVI(Normalized Difference Vegetation Index)：归一化差分植被指数，等于近红外波段的反射值与红光波段的反射值之差比上两者之和，计算简单，常被广泛应用于提取遥感影像中的植被。

2.1MBI计算

MBI是通过建立了建筑的光谱和空间特征与形态学算子之间的关系，实现建筑物的提取。

具体的计算步骤如下：

2.1.1获取遥感影像上各像元在可见光波段中光谱值的最大值，作为该像元的亮度值，从而获得亮度影像：

式(1)中，M_k(x)为像元x在波段k中的光谱值，b(x)则代表像元x的亮度值，K为总波段数。

2.1.2利用差分形态学属性的白顶帽值表示建筑的光谱和结构特征，基于多尺度和多方向线性结构元素的差分形态学属性白顶帽值定义为：

TH-DMP(d，s)＝|TH_b(d，s)-TH_b(d，s-Δs)| (2)

其中

代表基于重建的开运算后的亮度影像b，s(s_min≤s≤s_max)和d分别表示线性结构元素SE的尺度和方向，Δs是尺度的间隔。

2.1.3计算形态学房屋指数MBI，即不同尺度和方向的基于差分形态学属性白顶帽值的均值：

其中D和S分别是方向和尺度的数目。方向数D推荐为4，尺度参数是由影像的空间分辨率和建筑的大小决定的。本发明中推荐尺度参数s_min＝2，s_max＝150，四个方向分别为45°、90°、135°和180°。

2.2Harris特征的计算

Harris是一种基于信号的角点检测算法。其主要思想是使用一个检测窗口，当该窗口沿各个方向作微小移动时，检测窗口内的平均能量是否发生变化，当该能量变化值超过设定的阈值时，就将窗口的中心像素点提取为角点。

2.2.1如2.1.1中，计算亮度影像：

2.2.2对亮度影像b分别求水平和垂直方向上的梯度

式中，X和Y分别表示对亮度影像b求水平和垂直方向上的梯度，梯度滤波器f＝(-2，-1，0，1，2)。

2.2.3对梯度积高斯滤波

式中，w表示高斯滤波器，滤波窗口内的元素权重不同，中心位置高，边缘位置低；X²、Y²和XY梯度运算均是基于对应的梯度像素进行的。

梯度积运算可以表示为g(x，y)＝h(x，y)×j(x，y)，x和y分别是梯度影像中的像素位置，则h(x，y)和j(x，y)是梯度积因子，表示求乘积的两个梯度因子对应位置的值，g(x，y)为对应像素位置的梯度积。高斯滤波器

u和v表示窗口大小，这里推荐采用5×5，σ取值为2。

上述各式是先对梯度影像求积，然后再对梯度积进行高斯滤波。

2.2.4建立度量函数

R＝AB-C²-k(A+B)² (10)

式中的k是一个系数，此处取0.04，AB和C²表示求积，指整个影像对应位置像素求积，与梯度积类似，而A+B(两个图像求和)则是表示影像对应位置上的像素求和，可以表示为：t(x，y)＝a(x，y)+b(x，y)，x和y分别是梯度影像中的像素位置，则a(x，y)和b(x，y)是分别表示A和B中(x，y)位置处的像元值，t(x，y)为对应像素位置的像素和。

R是度量角点特征的响应指标，其值大于某个阈值则认为是角点。

2.3计算归一化植被指数NDVI

NDVI可以较好的识别植被，在遥感图像处理中广泛使用，其定义为：

式中的ρ_nir和ρ_red分别表示遥感影像中的近红外波段和红色波段的像元值。

步骤3、建筑特征中去除植被区域

分别对两个时相的NDVI图进行二值化处理得到前后时相植被覆盖区域，然后分别使用得到的植被区域过滤两个时相相应的建筑特征(MBI、Harris)，将植被区域赋值为背景值；

3.1求NDVI阈值

两个时相NDVI图均为单波段影像，分别剔除背景值后，可以对像元值进行统计，得到所有有效像元的均值μ和标准差σ，则自适应阈值T＝μ+σ，经过以下处理则得到二值图：

g₁′(x，y)和g₁(x，y)分别代表了二值化后的像元值和原始NDVI值，即NDVI影像中大于阈值的像元为1，表示植被区域，其余的值为0，代表背景值。

3.2DNVI过滤建筑指数中植被区域

NDVI影像取阈值后，得到只含有0(非植被)和1(植被)的二值图，将二值图中像元值为1的区域提取出来，在建筑指数(MBI和Harris)图中找到相应的位置，并赋为背景值，即完成了建筑指数中植被剔除：

式中g₁′(x，y)表示NDVI二值图的像元值，h₁(x，y)和h₁(x，y)分别表示建筑指数(MBI和Harris)在处理前后的像元值。

步骤4、计算特征热点图：

两个时相对应建筑特征值分别求差，然后对求差后的特征变化图进行半重叠开窗计算求窗口内平均亮度值，最后得到两个特征的热点图：IMBI和IHarris；

4.1两个时相建筑指数加权求差

其中，j∈{MBI，Harris}，w^j表示由两个时相的特征计算得到的加权值，

和

分别表示两个时相的特征，Dj表示两个特征加权求差后的影像值。

4.2半重叠开窗计算得到热点图

首先，需要确定窗口大小，设最小检测单元的面积为Area_min，窗口大小m：

式中，res为影像的分辨率，int是对结果求整。

半重叠滤波：

式中，j∈{MBI，Harris}，I^j为半重叠后的加权差值特征图，记为IMBI和IHarris，

表示使用滤波器对全图进行半重叠滤波运算，f_mean为平均滤波算子，由开窗的大小决定，大小为m×m，每个元素值为1/m。

经过以上的处理后，突出了两个时相建筑信息的差异。

步骤5、特征热点图取阈值：

对两个特征热点取阈值得到其对应二值图，即每个特征确定的新增建筑；

5.1IMBI二值化

先将IMBI影像中最小值和最大值线性拉伸到0-1：

式中，I_max和I_min表示IMBI影像中像素的最大值和最小值，g₂(x，y)和g₂′(x，y)分别表示线性拉伸前后的像素值，对拉伸后的影像取固定的阈值(推荐阈值为0.59)：

NL1(x，y)表示IMBI取阈值后的二值图，该结果就是基于MBI得到的新增建筑。

5.2IHarris二值化

依据本方法的实验经验，阈值按如下的方法确定：

T_top＝sort_0.98(IH) (20)

T_low＝sort_0.94(IH) (21)

式中，IH为IHarris强度图像中的所有有效的像元值，即非零值，sort_0.98(IH)表示将所有像元按照升序的方式排列，取第98％位置处的值，sort_0.94(IH)则表示取第94％位置处的值，T_h则是最后的Harris阈值。IHarris的二值化方法为：

NL2(x，y)为IHarris取阈值后的二值图，该结果就是基于Harris特征得到的新增建筑斑块，h₂(x，y)指IHarris特征。

步骤6：对二值化后的特征图取连通域得到对象，采用对象投票处理得到最终新增建筑斑块；

本过程包含对二值图取连通域得到建筑物对象，然后采用面向对象的方法投票，得到新增建筑对象。

6.1提取连通域

取NL1(由MBI确定的新增建筑)和NL2(由Harris确定的新增建筑)的4方向连通域，则有

其中，i∈[1，N_M]和j∈[1，N_H]，N_M和N_H分别表示NL1和NL2的连通域数目，

和

分别表示利用MBI和Harris特征识别的新增建筑对象斑块，它们是由相邻关系的点构成的对象，Δc和Δr由连通域的大小所决定，c0和r0表示连通域的起点位置。

6.2对象投票

经过连通域的提取后，实现了变化信息由像素层面到对象层面的转化，由于

对象可以较完整的包含新增建筑轮廓信息，但会有伪检测结果，而

在建筑上可以产生较强的反映但是边缘轮廓不完整，故本发明创新性的采用了对象投票方式来完成

对

中伪新增建筑的剔除：

式中，Num(*)表示对象的像素数目，

表示两个按照位置求交集，Tp为对象投票比例参数，按照经验取值，新增建筑提取设置为45％。

对象投票思路就是直接针对于对象的处理，若两个对象的交集像元数占

像元数比例小于某个值，就直接删除该对象。经过对象投票剔除后，得到新增建筑

i∈[1，N1_M]，N1_M是经过对象投票后的对象数。

步骤7：新增建筑后处理，主要包括前期建设用地现状图和最小制图单元滤除。

1)使用前期建设用地现状图过滤结果：

使用前期建设用地现状图过滤掉新增建筑与建设用地重叠的部分，得到更新后的新增建筑斑块。现状图是只含有0(非建设用地)和1(建设用地)的二值图，过滤方法如下：

式中f₁(x，y)表示前期现状图中的像元值，h₃(x，y)和h₃(x，y)分别表示新增建筑斑块

处理前后的像元值。

经过以上步骤的处理，得到的新增建筑结果记为

i∈[1，N2_M]，N2_M是经过前期建设用地过滤后的对象数。

2)最小制图单元过滤：

经过以上步骤的处理，新增建筑斑块正确性得到提高，但是小于检测单元的对象仍然存在，故需要将面积小于最小制图单元的对象剔除掉，方法如下：

此处的过滤仍然是采用面向对象的处理方式：

式中，res为影像的分辨率，Area_min为最小检测单元。

经过上述步骤处理后，得到了最终的新增建筑

i∈[1，N3_M]，N3_M是经过剔除对象后的对象数。

Claims (9)

1.一种高分辨率遥感影像新增建筑自动识别方法，其特征在于，包含以下步骤：

步骤1、对不同时相的高分辨率影像进行几何精纠正和相对辐射校正；

步骤2、获取前后两个时相影像的特征：形态学建筑指数MBI、角点特征Harris、归一化差分植被指数NDVI；

步骤3、建筑特征中去除植被区域：

分别对两个时相的NDVI图进行二值化处理得到前后时相植被覆盖区域，然后分别使用得到的植被区域过滤两个时相相应的建筑特征MBI和Harris，将植被区域赋值为背景值；

步骤4、计算特征热点图：

两个时相对应建筑特征值分别求差，然后对求差后的特征变化图进行半重叠开窗计算求窗口内平均亮度值，最后得到两个特征的热点图：IMBI和IHarris；

步骤5、特征热点图取阈值：

对两个特征热点取阈值得到其对应二值图，即每个特征确定的新增建筑；

步骤6、对二值化后的特征图取连通域得到对象，采用对象投票处理得到新增建筑斑块；

步骤7、新增建筑后处理，包括基于前期建设用地现状图和最小制图单元滤除。

2.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤2中形态学建筑指数MBI是通过建立建筑的光谱和空间特征与形态学算子之间的关系，实现建筑物的提取，具体的计算步骤如下：

2.1.1，获取遥感影像上各像元在可见光波段中光谱值的最大值，作为该像元的亮度值，从而获得亮度影像：

式(1)中，M_k(x)为像元x在波段k中的光谱值，b(x)则代表像元x的亮度值，K为总波段数；

2.1.2，利用差分形态学属性的白顶帽值表示建筑的光谱和结构特征，基于多尺度和多方向线性结构元素的差分形态学属性白顶帽值定义为：

TH-DMP(d，s)＝|TH_b(d，s)-TH_b(d，s-Δs)| (2)

其中

代表基于重建的开运算后的亮度影像b，s和d分别表示线性结构元素SE的尺度和方向，Δs是尺度的间隔；

2.1.3，计算形态学房屋指数MBI，即不同尺度和方向的基于差分形态学属性白顶帽值的均值：

其中D和S分别是方向和尺度的数目。

3.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤2中角点特征Harris是一种基于信号的角点检测算法，其主要思想是使用一个检测窗口，当该窗口沿各个方向作微小移动时，检测窗口内的平均能量是否发生变化，当该能量变化值超过设定的阈值时，就将窗口的中心像素点提取为角点，具体计算过程如下，

2.2.1，计算亮度影像：

M_k(x)为像元x在波段k中的光谱值，b(x)则代表像元x的亮度值，K为总波段数；

2.2.2，对亮度影像b分别求水平和垂直方向上的梯度：

式中，X和Y分别表示对亮度影像b求水平和垂直方向上的梯度，梯度滤波器f＝(-2，-1，0，1，2)；

2.2.3，对梯度积高斯滤波，

式中，w表示高斯滤波器，滤波窗口内的元素权重不同，中心位置高，边缘位置低；X²、Y²和XY梯度运算均是基于对应的梯度像素进行的；

2.2.4，建立度量函数，

R＝AB-C²-k′(A+B)² (10)

式中的k′是一个系数，AB和C²表示求积，指整个影像对应位置像素求积，与梯度积类似，而A+B则是表示影像对应位置上的像素求和；R是度量角点特征的响应指标，其值大于某个阈值则认为是角点。

4.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤2中归一化植被指数NDVI的计算方式如下，

式中的ρ_nir和ρ_red分别表示遥感影像中的近红外波段和红色波段的像元值。

5.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤3的具体实现方式如下，

首先采用一种自适应的阈值进行二值化处理，两个时相NDVI图均为单波段影像，对像元值进行统计，得到所有有效像元的均值μ和标准差σ，则自适应阈值T＝μ+σ，经过以下处理则得到二值图：

g₁′(x，y)和g₁(x，y)分别代表了二值化后的像元值和原始NDVI值，即NDVI影像中大于阈值的像元为1，表示植被区域，其余的值为0，代表背景值。

NDVI影像取阈值后，得到只含有非植被和植被的二值图，将二值图中像元值为1的区域提取出来，在建筑指数(MBI和Harris)中找到相应的位置，并赋为背景值，即完成了建筑指数中植被剔除：

式中g₁′(x，y)表示NDVI二值图的像元值，h₁(x，y)和h₁(x，y)分别表示建筑指数(MBI和Harris)在处理前后的像元值。

6.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤4的具体实现方式如下，

4.1，两个时相建筑指数加权求差，

其中，j∈{MBI，Harris}，w^j表示由两个时相的特征计算得到的加权值，

和

分别表示两个时相的特征，D^j表示两个特征加权求差后的影像值；

4.2，半重叠开窗计算得到热点图，

首先，确定窗口大小，设最小检测单元的面积为Area_min，窗口大小m：

式中，res为影像的分辨率，int是对结果求整；

半重叠滤波：

式中，j∈{MBI，Harris}，I^j为半重叠后的加权差值特征图，记为IMBI和IHarris，

表示使用滤波器对全图进行半重叠滤波运算，f_mean为平均滤波算子，由开窗的大小决定，大小为m×m，每个元素值为1/m；

经过以上的处理后，突出了两个时相建筑信息的差异。

7.如权利要求1所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤5的具体实现方式如下，

5.1，IMBI二值化，

先将IMBI影像中最小值和最大值线性拉伸到0-1：

式中，I_max和I_min表示IMBI影像中像素的最大值和最小值，g₂(x，y)和g₂′(x，y)分别表示线性拉伸前后的像素值，对拉伸后的影像取固定的阈值(推荐阈值为0.59)：

NL1(x，y)表示IMBI取阈值后的二值图，该结果就是基于MBI得到的新增建筑。

5.2，IHarris二值化，

IHarris阈值按如下的方法确定：

T_top＝sort_0.98(IH) (20)

T_low＝sort_0.94(IH) (21)

式中，IH为IHarris强度图像中的所有有效的像元值，即非零值，sort_0.98(IH)表示将所有像元按照升序的方式排列，取第98％位置处的值，sort_0.94(IH)则表示取第94％位置处的值，T_h则是最后的Harris阈值，IHarris的二值化方法为：

NL2(x，y)为IHarris取阈值后的二值图，该结果就是基于Harris特征得到的新增建筑斑块，h₂(x，y)是指IHarris特征。

8.如权利要求7所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤6的具体实现方式如下，

6.1，提取连通域

取NL1和NL2的4方向连通域，其中NL1是由MBI确定的新增建筑，NL2是由Harris确定的新增建筑，则有：

其中，i∈[1，N_M]和j∈[1，N_H]，N_M和N_H分别表示NL1和NL2的连通域数目，

和

分别表示利用MBI和Harris特征别的新增建筑对象斑块，它们是由相邻关系的点构成的对象，Δc和Δr由连通域的大小所决定，c0和r0表示连通域的起点位置；

6.2，对象投票

采用对象投票方式来完成

对

中伪新增建筑的剔除：

式中，Num(*)表示对象的像素数目，

表示两个按照位置求交集，Tp为对象投票比例参数；经过对象投票剔除后，得到新增建筑

i∈[1，N1_M]，N1_M是经过对象投票后的对象数。

9.如权利要求8所述的高分辨率遥感影像新增建筑自动识别方法，其特征在于：步骤7的具体实现方式如下，

1)使用前期建设用地现状图过滤结果：

使用前期建设用地现状图过滤掉新增建筑与建设用地重叠的部分，得到更新后的新增建筑斑块，现状图是只含有0和1的二值图，其中，0表示非建设用地，1表示建设用地，过滤方法如下：

式中f₁(x，y)表示前期建设用地现状图中的像元值，h₃(x，y)和h₃(x，y)分别表示新增建筑斑块

处理前后的像元值；

经过以上步骤的处理，得到的新增建筑结果记为

i∈[1，N2_M]，N2_M是经过前期建设用地过滤后的对象数；

2)最小制图单元过滤：

此处的过滤采用面向对象的处理方式：

式中，res为影像的分辨率，Area_min为最小检测单元；

经过上述步骤处理后，得到最终的新增建筑

i∈[1，N3_M]，N3_M是经过最小制图单元过滤后的对象数。

Images