CN106650663B - 建筑物真伪变化的判定方法及含此方法的伪变化去除方法 - Google Patents

Abstract

本发明提供一种建筑物真伪变化的判定方法，具体是：首先，从阴影信息入手基于空间拓扑一致性分析确定某城区的不同时相高分辨率遥感影像中对应建筑物及阴影的位置关系；然后，为具有对应关系的建筑物构造其局部建筑区域对，这些区域对均包含对应的建筑顶部边缘；最后，基于局部建筑区域内建筑顶部边缘结构的特征描述建立真伪变化判定指数，比较对应建筑物的局部建筑区域是否含有相似建筑顶部边缘结构信息，进而判定建筑物是否发生伪变化。应用本发明方法，能对建筑物的真伪变化进行精准判定。本发明还公开一种伪变化去除方法，结合上述建筑物真伪变化的判定方法，去除的伪变化像元可达90％以上，并且无误去除真变化像元现象，取得好的效果。

Description

建筑物真伪变化的判定方法及含此方法的伪变化去除方法

技术领域

本发明属于遥感影像处理技术领域，具体涉及一种建筑物真伪变化的判定方法及含此方法的伪变化去除方法。

背景技术

地球表面的地物(比如土地、建筑物、道路等)，由于自然原因或人为因素，随着时间的推移每天都发生着不同程度的变化。由于遥感手段本身所具备的远距离感知能力，利用搭载在不同平台上(如卫星、飞机、飞艇等)的各类传感器，对全球、国家、地区等不同尺度范围内的地物进行持续的监测，已成为各国行业部门的经常采用的非常有效的技术手段。

随着发展中国家的城市化进程的推进，尤其是我国新型城镇化政策的加深落实，城市区域内将面临着房屋、道路、桥梁等诸多人工地物的变迁。随着传感器技术的发展，遥感影像的空间分辨率能力已经得到跨越式提高，目前我们可以很方便地通过专门途径获取到1米以下甚至厘米级的高分辨率遥感影像，这为各国开展地表变化监测提供了非常有用的数据源。

但是，由于地物光谱特征的多样性、地物空间关系的复杂性、地物形状姿态的不确定性等现象普遍存在于遥感影像中，这严重影响到利用高分辨率遥感影像进行城市地物变化检测的精度。为此，相关领域的专家不断提出新的变化检测方法，来提高利用不同时相(大多是两个时相)的高分辨率遥感影像进行变化检测的精度。但有些因素的客观存在，会一直影响到变化检测精度的提高，其中影响最为严重的就是所谓“伪变化”的存在。“伪变化”指本身未发生变化而检测出的虚假变化。导致这种伪变化的一个主要原因是因传感器的观测角不同而导致相同建筑物在不同时相遥感影像中表现出不同的姿态。这是一般针对地物光谱特征多样性问题的变化检测方法所不能解决的。

为了实现伪变化的消除，相关技术人员做了研究，详情如下：

Wang等人提出的波段融合法，即两时相交叉融合图像中错位区会产生“虚影”，因此，将原图像以及波段融合图像叠加统一分割得到对象后，处于错位区的对象可通过设计的截距、斜率及相关系数三个特征指数分类判别得出。此种方法中的波段融合法需要事先获取实验区域的不同时相的全色影像和多波段影像，对数据要求较高，而且将原图与波段融合图像叠加分割获取统一对象时，有可能因几何配准精度不高以及阴影影响，得到一些“伪变化对象”，这些对象无法通过该文章设计的指数消除；此外，这个方法用来判别变化对象的三个指数是由影像的灰度值计算得到的，若不同时相影像内同一建筑物颜色有差异，该方法也可能将其判断为变化建筑物，出现在最终的变化检测结果中。

Tang等人提出了一种将MBI与Harris角点结合构造建筑物感兴趣点，即通过建筑物感兴趣点的匹配，辨别出在不同传感器拍摄角度得到的不同时相影像中的同一建筑物，从而达到减轻误检现象的目的。此种方法由于使用叠加判断来构造建筑物感兴趣点，易受MBI提取得到的建筑物区域以及Harris角点提取结果精度的影响，可能会出现同一建筑物因某时相影像内部建筑物感兴趣点缺失而无法正确匹配上的问题。

综上所述，急需一种步骤精简、操作方便且伪变化消除效果好的方法以解决现有技术中存在的问题。

发明内容

本发明目的在于提供一种步骤精简、操作方便且针对不同时相建筑物伪变化信息的判定方法，具体技术方案如下：

一种建筑物真伪变化的判定方法，包括以下步骤：

第一步、提取不同时相高分辨率遥感影像中的建筑物阴影信息；

第二步、根据第一步所得建筑物阴影信息获取阴影形心点集，建立不同时相影像内建筑物阴影的空间拓扑匹配关系，确定不同时相影像内具有相似拓扑结构的建筑物阴影；

第三步、对于不同时相内相似拓扑结构的每个节点上的建筑物阴影通过其主方向信息提取出包含部分阴影与部分建筑物边缘的区域，获得一系列局部建筑区域对；

第四步、构建判别指数，通过判别指数判断不同时相对应局部建筑区域内是否包含相同建筑物边缘：若不同时相对应局部建筑区域内不包含相同建筑物边缘，则判定该局部建筑区域对应的建筑物真实发生变化；若不同时相对应局部建筑区域内包含相同建筑物边缘，则判断该局部建筑区域对应的建筑物真实未发生变化。

以上技术方案中优选的，所述第一步中提取不同时相高分辨率遥感影像中的建筑物阴影信息包括以下步骤：

步骤1.1、提取建筑物的高分辨率遥感影像中每个像元的光谱特征和纹理特征，所述光谱特征包括影像像元各波段上的亮度值以及和两个参数，其中：和的计算表达式详见表达式1)和表达式2)：

其中：H、I为将原红、绿、蓝三波段影像转换到HIS颜色空间后的H分量和I分量，MSIIndex为形态学阴影指数；

所述纹理特征采用Gabor算子提取，每个像元可获得一个24维Gabor纹理特征向量；

步骤1.2、采集影像中的部分阴影像元作为训练样本，利用步骤1.1提取得到的光谱特征和纹理特征，对影像进行像元级别的分类，获得阴影掩膜图；

步骤1.3、提取阴影块，具体是：先对步骤1.2得到的阴影掩膜图进行形态学处理；再依据四邻域邻接关系，提取阴影掩模图中的连通分量，视其为一块阴影对象；

步骤1.4、筛选出建筑物产生的阴影，具体是：结合面积和形状指数对步骤1.3所判断得到的阴影对象进行筛选，其中：面积A为每块阴影所占像元个数；形状指数SI通过表达式3)获得，表达式3)如下：

其中：SI为形状指数，A为阴影块面积，L为阴影区域的最小外接矩形的较长边。

以上技术方案中优选的，所述步骤1.4中结合面积和形状指数对步骤1.3所判断出的阴影对象进行筛选详细过程是：当某块阴影的面积A和形状指数SI满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影为建筑物阴影；当某块阴影的面积A和形状指数SI不同时满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影不属于建筑物阴影；其中：A₁、A₂、SI₁和SI₂均为人工设定的阈值，A₁的取值范围为[50，200]，A₂的取值>1000，SI₁的取值范围为[0.05,0.1]，SI₂的取值>0.25。

以上技术方案中优选的，所述第二步中建立不同时相影像内建筑物阴影的匹配关系采用点集匹配算法；

所述第二步中提取两时相影像内具有相似拓扑结构的建筑物阴影的详细过程包括以下步骤：

步骤2.1、获取阴影形心点集，具体是：将每一块阴影块的影像行列号视为横纵坐标，求取该阴影块内部像元的行列号均值，获取该阴影块的形心点；集合所有阴影块的形心点，即得该影像的阴影形心点集；

步骤2.2、从两时相影像阴影形心点集P、Q中各自找到一个含有K个点的子点集P₀和Q₀，其中：K≤min(N，M)，N和M分别为P和Q中的点的个数；对于P中的某一点p_i，对Q中所有点都通过算法计算相似性指标；

步骤2.3、进行判断，具体是：设Q中拥有相似性指标值最大的点位为q_j，则q_j为p_i在点集Q中的最优匹配点；如果p_i亦为q_j在点集P中的最优匹配点，则此时p_i与q_j匹配成功，组成匹配对；每组匹配对中点的序号即为算法找到的点集序号匹配对{i，j}，这些匹配对指示出了同一建筑物在两时相影像中的阴影块。

以上技术方案中优选的，所述第三步具体包括以下步骤：

步骤3.1、阴影主方向提取，具体是：首先，采用最小周长多边形简化阴影轮廓；其次，计算该最小周长多边形的主方向；最后，以阴影的主方向角度的正切值为斜率，过阴影质心画出线段，得到该阴影的主方向；

步骤3.2、提取局部建筑区域，具体是：将扫描线段朝建筑物侧平移固定距离，以包含建筑物与阴影交界边缘部分，即可得到矩形的局部建筑区域。

以上技术方案中优选的，所述步骤3.2中扫描线段的获取方法是：a、先获取阴影块的原始阴影轮廓线，后获取其最小周长多边形及其质心；b、通过原始阴影轮廓线求得阴影区的最小外接矩形；c、通过最小周长多边形质心，以主方向角度的正切值为斜率画直线；d、截取矩形内的直线部分作为最终扫描线段。

以上技术方案中优选的，当阳光照射角度发生变化时，所述步骤3.2中扫描线的移动方向通过以下方法获得：首先、获得阳光照射方向的反方向α；其次、以反方向α为角平分线，得到搜索范围[α-90°，α+90°)；最后、在搜索范围找到与主方向垂直的唯一角度θ，作为平移方向。

以上技术方案中优选的，所述第四步中利用方向纹理曲线和方向梯度直方图特征算子构建判别指数；

所述第四步的具体过程是：

步骤4.1、取沿-90°带状区域划分的方向纹理曲线；

步骤4.2、依据-90°带状区域划分方式，以每一行为带状区域，获取局部建筑区域的方向纹理曲线；

步骤4.3、寻找到梯度曲线的最大值点Gra_max，并定义梯度变化截断阈值Gra_cut；

步骤4.4、在坐标轴内画出Gra＝Gra_cut的横线，获取得到截断最大值点所在波峰部分所覆盖的带状区域编号，并提取出其中的区域作为Cell；

步骤4.5、提取出两时相对应的局部建筑区域中Cell后，比较判定该局部建筑区域对应的不同时相建筑物真实是否未发生变化。

以上技术方案中优选的，所述步骤4.5中比较和判定过程如下：

首先，求取Cell的Hog特征向量；

其次，比较两个Cell的Hellinger距离；

最后，判定：若两个Cell的Hog特征向量的Hellinger距离不超过人工设定阈值G，则判定为包含相同的建筑物边缘，继而判断该局部建筑区域对应的建筑物真实未发生变化；若两个Cell的Hog特征向量的Hellinger距离超过阈值G，则判定不包含相同的建筑物边缘，继而判定该局部建筑区域对应不同时相建筑物真实发生变化。

本发明所公开的一种建筑物(尤其是高分辨率遥感影像建筑物)真伪变化的判定方法，具体是：首先，从阴影信息入手基于空间拓扑一致性分析确定某城区的不同时相高分辨率遥感影像中对应建筑物及阴影的位置关系；然后，为具有对应关系的建筑物构造其局部建筑区域对，这些区域对均包含对应的建筑顶部边缘；最后，基于局部建筑区域内建筑顶部边缘结构的特征描述建立真伪变化判定指数，比较对应建筑物的局部建筑区域是否含有相似建筑顶部边缘结构信息，进而判定建筑物是否发生伪变化。应用本发明方法，效果是：

(1)对数据要求不高，无需获取不同时相的全色影像和多波段影像，即使是普通的RGB波段高分辨率遥感影像也能达到一定效果。

(2)本发明基于阴影的建筑物空间拓扑关系，能够很好地体现出某城市建成区内的建筑物空间关系。

(3)本发明局部建筑区域对包含相应的建筑物局部边缘，能够通过判断建筑物边缘结构的变化来反映出建筑物变化检测中出现的伪变化。

(4)本发明方法中提出的伪变化判别指标，不易受少量的杂物、噪声影响，对颜色变化也有一定的容忍能力。

(5)可以作为任意高分辨率遥感影像面向对象变化检测方法的后处理方法，去除其中的伪变化成分，从而提高结果精度。

本发明还公开一种伪变化的去除方法，包括以下步骤：

步骤一、通过变化检测，得到变化区域二值图；通过上述的建筑物真伪变化的判定方法对建筑物的真伪变化进行判定；

步骤二、将步骤一种所得的变化区域二值图和对建筑物的真伪变化进行判定的结果进行二值图求交，消除变化区域二值图中的伪变化。

采用本发明的伪变化去除方法，步骤精简，操作方便，去除的伪变化像元超过90％，并且没有去除掉真变化像元，伪变化去除效果好。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是实施例1中形态学阴影指数MSIIndex的提取流程示意图；

图2是实施例1中点集匹配算法效果示意图；

图3是实施例1中利用主方向信息寻找建筑物边缘区域的示意图；

图4是实施例1中主方向不为0°的阴影构建得到的矩形区域提取方法示意图；

图5是实施例1中参数构建结构示意图；

图6是实施例2影像信息示意图；

图6(a)为1999、2010年两时相影像示意图；

图6(b)为人工标注的两时相影像中的建筑物掩膜示意图；

图6(c)为边缘结构发生变化的1、2号建筑物示意图；

图7是实施例2阴影提取效果以及匹配结果示意图；

图7(a)为1999、2010年影像中提取并经处理过后的阴影示意图；

图7(b)为经匹配流程后剩下的具有对应关系的阴影块示意图；

图8是实施例2中矩形区构建以及局部建筑区域提取效果示意图；

图8(a)为所有匹配对阴影构造得到的矩形区示意图；

图8(b)为图8(a)中的1、2和3号建筑物示意图；

图8(c)为1、2和3号建筑物获取得到的局部建筑区域示意图；

图9是实施例2中参数判别结果示意图；

图9(a)为经参数判别后保留的局部建筑区域匹配对示意图；

图9(b)为图9(a)中对应的一些建筑物示意图；

图10是实施例2中叠加判断示意图；

图11是实施例2中伪变化去除流程示意图；

图11(a)为图6(a)中两时相影像通过MBI算法检测得到的建筑物变化区域示意图；

图11(b)为“伪变化”区域与变化区域的叠加图；

图11(c)中为实际上没有发生变化的建筑物产生的伪变化示意图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

实施例1：

一种建筑物真伪变化的去除方法，包括以下步骤：

步骤一、通过常规的变化检测方法，得到变化区域二值图；对建筑物的真伪变化进行判定；

步骤二、将步骤一种所得的变化区域二值图和对建筑物的真伪变化进行判定的结果进行二值图求交，去除变化区域二值图中的伪变化对象。

上述步骤一中对建筑物的真伪变化进行判定，包括四大步骤，详细过程如下：

第一步、提取不同时相高分辨率遥感影像中的建筑物阴影信息，以两时相为例，具体包括以下步骤：

步骤1.1、提取建筑物的高分辨率遥感影像中每个像元的光谱特征和纹理特征，所述光谱特征包括影像像元各波段上的亮度值以及和两个参数，其中：和的计算表达式详见表达式1)和表达式2)：

其中：H、I为将原红、绿、蓝三波段影像转换到HIS颜色空间后的H分量和I分量，MSIIndex为形态学阴影指数(其提取流程详见图1)。

纹理特征采用Gabor算子提取，每个像元可获得一个24维的Gabor纹理特征向量。

步骤1.2、采集影像中的部分阴影像元作为训练样本，利用步骤1.1提取得到的光谱、纹理特征，对影像进行像元级别的分类，获得阴影掩膜图，此处采用的监督分类器为ELM(极限学习机)。

步骤1.3、提取阴影块，先对步骤1.2得到的阴影掩膜图进行形态学处理；再依据四邻域邻接关系，提取阴影掩模图中的连通分量，视其为一块阴影对象。

步骤1.4、筛选出建筑物产生的阴影，具体是：结合面积和形状指数对步骤1.3所判断得到的阴影对象进行筛选，其中：面积A为每块阴影所占像元个数；形状指数SI通过表达式3)获得，表达式3)如下：

其中：A为阴影块面积，L为阴影区域的最小外接矩形的较长边。

所述步骤1.4中结合面积和形状指数对步骤1.3所判断出阴影对象进行筛选详细过程是：当某块阴影的面积A和形状指数SI满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影为建筑物阴影；当某块阴影的面积A和形状指数SI不同时满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影不属于建筑物阴影；其中：A₁、A₂、SI₁和SI₂均为人工设定的阈值，A₁一般取[50，200]，A₂的取值>1000；而SI₁的取值范围为[0.05,0.1]，SI₂的取值>0.25。

第二步、根据第一步所得建筑物阴影信息获取阴影形心点集，利用点集匹配算法建立不同时相影像内建筑物阴影的空间拓扑匹配关系，确定不同时相影像内具有相似拓扑结构的建筑物阴影，详细过程包括以下步骤：

步骤2.1、获取阴影形心点集，具体是：将每一块阴影块的影像行列号视为横纵坐标，求取该阴影块内部像元的行列号均值，获取该阴影块的形心点；集合所有阴影块的形心点，即得该影像的阴影形心点集；

步骤2.2、从两时相影像阴影形心点集P、Q中各自找到一个含有K个点的子点集P₀和Q₀，其中：K≤min(N，M)，N和M分别为P和Q中的点的个数；对于P中的某一点p_i，对Q中所有点都通过算法计算相似性指标；

步骤2.3、进行判断，具体是：设Q中拥有相似性指标值最大的点位为q_j，则q_j为p_i在点集Q中的最优匹配点；如果p_i亦为q_j在点集P中的最优匹配点，则此时p_i与q_j匹配成功，组成匹配对；每组匹配对中点的序号即为算法找到的点集序号匹配对{i，j}，这些匹配对指示出了同一建筑物在两时相影像中的阴影块。P₀和Q₀组成的网络的拓扑结构非常相似，如图2所示，若将其中一点集(P₀)整体做简单的平移、旋转操作后，P₀和Q₀的相应节点都处在一个恰好包含来自两个形心点集的各一个点的圆内，圆的半径即由参数δ决定。因此参数δ也控制了匹配对的数目(即节点数目)。δ一般从给定初始值(如取1)，以给定步长(如取0.5)直至某个终止值(如取10)取离散值，记录下寻找得到的匹配对数n，当n获得最大的值时则为δ的最终取值，图2中：小圆圈和红叉分别代表两个不同点集P与Q，大圆圈中圈出的为成功匹配结果。

第三步、对于不同时相内相似拓扑结构的每个节点上的建筑物阴影通过其主方向信息提取出包含部分阴影与部分建筑物边缘的区域，获得一系列局部建筑区域对，具体如下；

步骤3.1、阴影主方向提取，具体是：首先，采用最小周长多边形简化阴影轮廓；其次，计算该最小周长多边形的主方向；最后，以阴影的主方向角度的正切值为斜率，过阴影质心画出线段，得到该阴影的主方向；

步骤3.2、提取局部建筑区域，具体是：将扫描线段朝建筑物侧平移固定距离，以包含建筑物与阴影交界边缘部分，即可得到矩形的局部建筑区域。所述扫描线段的获取方法是：a、先获取阴影块的原始阴影轮廓线，后获取其最小周长多边形及其质心；b、通过原始阴影轮廓线求得阴影区的最小外接矩形；c、通过最小周长多边形的质心，以主方向角度的正切值为斜率画出直线；d、截取矩形内的直线部分作为最终扫描线段。

当阳光照射角度发生变化时，扫描线的移动方向通过以下方法获得：首先、获得计算阳光照射方向的反方向α；其次、以反方向α为角平分线，得到搜索范围[α-90°，α+90°)；最后、在搜索范围找到与主方向垂直的唯一角度θ，作为平移方向，详见图3，图3中处于水平方向的实线为通过阴影主方向寻找得到的扫描线，θ为扫描线平移方向，实线外框为将扫描线平移一定距离d后得到保护部分建筑物边缘的矩形区域。旋转提取区域示意图如图4所述，图4中：a图为建筑物示意图，b图为阴影轮廓及生产扫描线，c图为生成矩形区域，d图和e图分别为将原图及矩形区域旋转置主方向为0的情况，f图为区域提取效果。

第四步、利用方向纹理曲线和方向梯度直方图特征算子构建判别指数，通过判别指数判断不同时相对应局部建筑区域内是否包含相同建筑物边缘进行判别：若不同时相对应局部建筑区域内不包含相同建筑物边缘，则判定该局部建筑区域对应的建筑物真实发生变化；若不同时相对应局部建筑区域内包含相同建筑物边缘，则判断该局部建筑区域对应的筑物真实未发生变化，详见图5，具体过程如下：

步骤4.1、取沿-90°带状区域划分的方向纹理曲线；

步骤4.2、依据-90°带状区域划分方式，以每一行为带状区域，获取局部建筑区域的方向纹理曲线；

步骤4.3、寻找到梯度曲线的最大值点Gra_max，并定义梯度变化截断阈值Gra_cut；

步骤4.4、在坐标轴内画出Gra＝Gra_cut的横线，获取得到截断最大值点所在波峰部分所覆盖的带状区域编号(图5中的St1和St2)，并提取出其中的区域作为Cell；

步骤4.5、提取出两时相对应的局部建筑区域中Cell后，进行比较，判定该局部建筑区域对应的不同时相建筑物真实是否未发生变化，所述比较和判定过程如下：

首先，求取Cell的Hog特征向量；

其次，比较两个Cell的Hellinger距离；

最后，判定：若两个Cell的Hog特征向量的Hellinger距离不超过人工设定阈值G，则判定为包含相同的建筑物边缘，继而判断该局部建筑区域对应的建筑物真实未发生变化；若两个Cell的Hog特征向量的Hellinger距离超过阈值G，则判定不包含相同的建筑物边缘，继而判定该局部建筑区域对应不同时相建筑物真实发生变化。其中：人工设定阈值G视HOG的bin数而定，此处：当HOG的bin数取5时，G的合适取值范围为0.1-0.2。

实施例2

本发明建筑物真伪变化的判定和伪变化的去除方法的实际应用。

所采用的影像为Google Earth影像，分辨率为0.12米，拍摄区域为美国华盛顿市区一块区域，影像行列数为1098×2476，摄取时间为1999年4月30日以及2010年4月3日。图6为影像信息，图6中：(a)为1999、2010年两时相影像，(b)为人工标注的两时相影像中的建筑物掩膜，深灰色(如标号为s1处)为1999年建筑物，灰色(如标号为s3处)为2010年建筑物，浅灰色(如标号为s2处)为重叠部分，灰色框(如标号为s4处)中的建筑物为实际发生了变化的建筑物，(c)为边缘结构发生变化的1、2号建筑物。

通过目视解译，1999年共有43栋建筑物，阳光入射角为92.5°；2010年共有42栋建筑物，阳光入射角为97.5°。相对于1999年的建筑物，2010年左上角的1、2号两栋建筑物顶部都进行了增建，并且边缘结构有所改变，如图6(c)所示；此外，上方靠左3号建筑物被拆除掉了，其余的建筑物都没有发生变化，但由于拍摄角度的不同，其建筑物掩膜部分都不完全重叠。

通过原图，便可通过前文提到的方法获取阴影，并提取阴影形心点集以建立起匹配关系。图7为阴影提取效果以及匹配结果，图7中：(a)为1999、2010年影像中提取并经处理过后的阴影；(b)为经匹配流程后剩下的具有对应关系的阴影块。

本实施例在分类得到阴影像元后，用于进行形态学处理的形态学开闭运算算子为半径为3的圆形。通过连通分量得到阴影块后，用于去除碎部、树木阴影的面积阈值A>＝T_A＝200，形状系数SI>T_SI＝0.1。流程提取得到的阴影如图7中(a)所示；匹配过程中形心点集运算设置的噪声参数δ＝10；图7中(b)为寻找到的匹配对阴影，匹配对数n＝41。

基于影像阴影分类结果，匹配结果能将所有位置不变的建筑物都匹配上了。没匹配上的两队在图中已用标号为1和2的框框出。其中1号框中的阴影没匹配上，一方面是因为框中得到的阴影块数目不同，其中1999年1号框中得到2块阴影，2010年框中得到了3块阴影，形心点数目不一致；且2010年框中左下角的阴影块由于树木遮挡的原因，缺失了一部分阴影，进一步造成了形心点的偏移；因此两时相1号框内的阴影形心点全都没有匹配上；而2号框中的阴影变化属于建筑物拆迁引起的，理应予以排除。

在得到匹配阴影对之后，下一步即为通过阴影构造矩形区，图8为矩形区构造效果，图8中：(a)为所有的阴影块构造得到的矩形区基本上很好的覆盖了建筑物与阴影的交界地带，得到的局部建筑区域块基本都形成了上文提到的上部为阴影，下部为建筑物，边缘基本成水平走势的情况；(b)和(c)为图中框出的一些建筑物提取得到的局部建筑区域块效果。上述在构造矩形区时，扫描线平移距离D_M＝60。

得到局部建筑区域块之后，便可通过参数计算来筛选掉一部分建筑物。图9为参数判别结果，图9中：(a)两时相参数判别而被认为存在有伪变化的矩形区域，剩下的匹配对p＝35，(b)为相对应的实际建筑物。计算判别参数的过程中，采用的梯度截断阈值Gra_cut＝5，获取Cell后，提取Cell的Hog向量的直方图bin数目Hbin＝5，计算得到相应Cell对的HDist之后，用于判别的阈值T₀＝0.15。

在参数判别过程中，丢失了部分实际上没有发生变化的建筑物，一方面是由于部分局部建筑区域中的阴影内存在树枝，另一方面是由于传感器拍摄角度的原因，2010年的建筑物侧面露出部分更多一些，而且建筑物的顶部存在一些存在颜色差异的斑点。在树枝、侧面杂物和顶部斑点的共同作用下，使得一部分匹配对的HDist取得了较大的值，没有通过判定。

得到图7中(a)中的矩形区域后，由于这些区域包含一部分建筑物边缘，而且经过了所设计的参数进行了相似性判定，因此，如果建筑物没有经过部分扩建的情况下，可以认为这些边缘所属的建筑物在两时相内没有发生变化。

由于在变化检测中，检测出来的因拍摄角度引起的“伪变化”同样会是不变建筑物的一部分。考虑到之前的流程已经通过点集匹配来证明这两个建筑物出现在同一拓扑网络的相同节点上，又通过参数判断出这两栋建筑物具有相似性质，从而可以认定这两个建筑物是不变建筑物。因此，这些建筑物产生的“伪变化”对象会有一部分出现在结果矩形区域中，可以将建筑物变化检测结果与上述过程得到的结果矩形区相叠加(见图10)，从而进行变化性质的判断以及去除。这种叠加判断方法适合于面向对象的变化检测结果。

为了直观表述本发明方法得到的结果如何作用于变化检测结果变化区域二值图，利用MBI指数分别获取两时相建筑物对象再进行比较获取建筑物变化区域。MBI指数的构建充分考虑到建筑物的亮度等光谱结构特征，以及局部对比度、形状、大小和方向性等特征，通过一系列形态学操作提取得到。其提取步骤大致分为以下四步：

步骤1：计算亮度值

其中，band_k(x)为第k光谱波段在像素x处的亮度值，K为可见光光谱波段数。将可见光波段的像素最大值作为该像素亮度值。

步骤2：形态学白帽重构

其中，为对亮度图像b的形态学开运算，而d和s分别代表线性结构元素的方向和尺度。由于建筑物与道路光谱类似，该指数建立的一个重点就是如何自动滤除道路，因为道路总是沿一个或者两个方向延伸，而建筑物则是多个方向。因此采用多个方向和尺度的线性结构元素能将二者很好地区分开来。

步骤3：计算微分形态学剖面DMP(Differential Morphological Profiles)

DMP_WTH(d,s)＝|WTH(d,(s+Δs))-WTH(d,s)|

其中，Δs为人工预设的线性结构元素尺度间隔。

步骤4：计算MBI指数

其中，S＝((s_max-s_min)/Δs)+1，s_min和s_max分别为通过量取遥感影像中房屋尺度所确定的最小、最大线性结构元素，D为计算建筑物剖面时的方向数。由于建筑物周围是亮度值较低的阴影,建筑物具有较大的局部对比度，又由于建筑物结构在微分形态学剖面DMP_W-TH大部分方向上具有较大的值，因此MBI值较大的区域是建筑物。

本实施例中参数：D＝6，s_min＝6，s_max＝42，ΔS＝4，由此得到的方向参数向量为d＝[0，30°，60°，90°，120°，150°]，尺度参数向量为s＝[6，10，14，18，22，26，30，34，38，42]。得到影像MBI值之后，通过MBI(x)>0.4获取两时相建筑物像元，经过形态学开闭运算处理后，在两时相影像中会得到一些图斑，再通过外接矩形长宽比<4，面积≥200提取建筑物对象。最后，通过两时相影像对比获取变化区域。图11即为通过MBI得到的两时相建筑物变化区域，以及本发明方法得到的伪变化区域处理对其处理的效果图。图11中(a)为图6(a)中两时相影像通过MBI算法检测得到的建筑物变化区域，(b)为本发明方法得到的“伪变化”区域与变化区域的叠加图，(c)为实际上没有发生变化的建筑物产生的伪变化。其中，浅灰色部分为成功检测出来部分，深灰色为未成功检测出来部分。

矩形区中包含部分建筑物边缘，我们将出现在“伪变化”区域中的建筑物视为未变化建筑物，因此，通过叠加运算后，可以去除掉一部分“伪变化”建筑物对象。图11中(c)中浅灰色对象为本发明方法能够去除的对象。去除效果如表1所示。

表1“伪变化”去除统计

变化检测方法：MBI	变化对象数	变化像元数
			去除前	84	376027
去除量	57	314783
			实际未发生变化量	75	341274
正确去除率	76％	92.24％

注：变化像元、对象包括所有真变化和未变化的像元对象，本发明方法去除的像元及对象都属于“伪变化”，并不包含真实发生的变化。

从表1可知，本发明方法去除的伪变化像元超过了90％，并且没有误去除掉真变化像元，取得了较好的效果。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种建筑物真伪变化的判定方法，其特征在于，包括以下步骤：

第一步、提取不同时相高分辨率遥感影像中的建筑物阴影信息；

第二步、根据第一步所得建筑物阴影信息获取阴影形心点集，建立不同时相影像内建筑物阴影的空间拓扑匹配关系，确定不同时相影像内具有相似拓扑结构的建筑物阴影；

第三步、对于不同时相内相似拓扑结构的每个节点上的建筑物阴影通过其主方向信息提取出包含部分阴影与部分建筑物边缘的区域，获得一系列局部建筑区域对；

第四步、构建判别指数，通过判别指数判断不同时相对应局部建筑区域内是否包含相同建筑物边缘：若不同时相对应局部建筑区域内不包含相同建筑物边缘，则判定该局部建筑区域对应的建筑物真实发生变化；若不同时相对应局部建筑区域内包含相同建筑物边缘，则判断该局部建筑区域对应的建筑物真实未发生变化。

2.根据权利要求1所述的建筑物真伪变化的判定方法，其特征在于，所述第一步中提取不同时相高分辨率遥感影像中的建筑物阴影信息包括以下步骤：

步骤1.1、提取建筑物的高分辨率遥感影像中每个像元的光谱特征和纹理特征，所述光谱特征包括影像像元各波段上的亮度值以及和两个参数，其中：和的计算表达式详见表达式1)和表达式2)：

其中：H、I为将原红、绿、蓝三波段影像转换到HIS颜色空间后的H分量和I分量，MSIIndex为形态学阴影指数；

所述纹理特征采用Gabor算子提取，每个像元可获得一个24维Gabor纹理特征向量；

步骤1.2、采集影像中的部分阴影像元作为训练样本，利用步骤1.1提取得到的光谱特征和纹理特征，对影像进行像元级别的分类，获得阴影掩膜图；

步骤1.3、提取阴影块，具体是：先对步骤1.2得到的阴影掩膜图进行形态学处理；再依据四邻域邻接关系，提取阴影掩模图中的连通分量，视其为一块阴影对象；

步骤1.4、筛选出建筑物产生的阴影，具体是：结合面积和形状指数对步骤1.3所判断得到的阴影对象进行筛选，其中：面积A为每块阴影所占像元个数；形状指数SI通过表达式3)获得，表达式3)如下：

其中：SI为形状指数，A为阴影块面积，L为阴影区域的最小外接矩形的较长边。

3.根据权利要求2所述的建筑物真伪变化的判定方法，其特征在于，所述步骤1.4中结合面积和形状指数对步骤1.3所判断出的阴影对象进行筛选详细过程是：当某块阴影的面积A和形状指数SI满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影为建筑物阴影；当某块阴影的面积A和形状指数SI不同时满足A∈[A₁，A₂]且SI∈[SI₁，SI₂]时，判定某块阴影不属于建筑物阴影；其中：A₁、A₂、SI₁和SI₂均为人工设定的阈值，A₁的取值范围为[50,200]，A₂的取值>1000，SI₁的取值范围为[0.05,0.1]，SI₂的取值>0.25。

4.根据权利要求1所述的建筑物真伪变化的判定方法，其特征在于，所述第二步中建立不同时相影像内建筑物阴影的匹配关系采用点集匹配算法；

所述第二步中提取两时相影像内具有相似拓扑结构的建筑物阴影包括以下步骤：

步骤2.1、获取阴影形心点集，具体是：将每一块阴影块的影像行列号视为横纵坐标，求取该阴影块内部像元的行列号均值，获取该阴影块的形心点；集合所有阴影块的形心点，即得该影像的阴影形心点集；

步骤2.2、从两时相影像阴影形心点集P、Q中各自找到一个含有K个点的子点集P₀和Q₀，其中：K≤min(N，M)，N和M分别为P和Q中的点的个数；对于P中的某一点p_i，对Q中所有点都通过算法计算相似性指标；

步骤2.3、进行判断，具体是：设Q中拥有相似性指标值最大的点位为q_j，则q_j为p_i在点集Q中的最优匹配点；如果p_i亦为q_j在点集P中的最优匹配点，则此时p_i与q_j匹配成功，组成匹配对；每组匹配对中点的序号即为算法找到的点集序号匹配对{i，j}，这些匹配对指示出了同一建筑物在两时相影像中的阴影块。

5.根据权利要求1所述的建筑物真伪变化的判定方法，其特征在于，所述第三步具体包括以下步骤：

步骤3.1、阴影主方向提取，具体是：首先，采用最小周长多边形简化阴影轮廓；其次，计算该最小周长多边形的主方向；最后，以阴影的主方向角度的正切值为斜率，过阴影质心画出线段，得到该阴影的主方向；

步骤3.2、提取局部建筑区域，具体是：将扫描线段朝建筑物侧平移固定距离，以包含建筑物与阴影交界边缘部分，即可得到矩形的局部建筑区域。

6.根据权利要求5所述的建筑物真伪变化的判定方法，其特征在于，所述步骤3.2中扫描线段的获取方法是：a、先获取阴影块的原始阴影轮廓线，后获取其最小周长多边形及其质心；b、通过原始阴影轮廓线求得阴影区的最小外接矩形；c、通过最小周长多边形质心，以主方向角度的正切值为斜率画直线；d、截取矩形内的直线部分作为最终扫描线段。

7.根据权利要求6所述的建筑物真伪变化的判定方法，其特征在于，当阳光照射角度发生变化时，所述步骤3.2中扫描线的移动方向通过以下方法获得：首先、获得阳光照射方向的反方向α；其次、以反方向α为角平分线，得到搜索范围[α-90°，α+90°)；最后、在搜索范围找到与主方向垂直的唯一角度θ，作为平移方向。

8.根据权利要求1所述的建筑物真伪变化的判定方法，其特征在于，所述第四步中利用方向纹理曲线和方向梯度直方图特征算子构建判别指数；

所述第四步的具体过程是：

步骤4.1、取沿-90°带状区域划分的方向纹理曲线；

步骤4.2、依据-90°带状区域划分方式，以每一行为带状区域，获取局部建筑区域的方向纹理曲线；

步骤4.3、寻找到梯度曲线的最大值点Gra_max，并定义梯度变化截断阈值Gra_cut；

步骤4.4、在坐标轴内画出Gra＝Gra_cut的横线，获取得到截断最大值点所在波峰部分所覆盖的带状区域编号，并提取出其中的区域作为Cell；

步骤4.5、提取出两时相对应的局部建筑区域中Cell后，比较判定该局部建筑区域对应的不同时相建筑物真实是否未发生变化。

9.根据权利要求8所述的建筑物真伪变化的判定方法，其特征在于，所述步骤4.5中比较和判定过程如下：

首先，求取Cell的Hog特征向量；

其次，比较两个Cell的Hellinger距离；

最后，判定：若两个Cell的Hog特征向量的Hellinger距离不超过人工设定阈值G，则判定为包含相同的建筑物边缘，继而判断该局部建筑区域对应的建筑物真实未发生变化；若两个Cell的Hog特征向量的Hellinger距离超过阈值G，则判定不包含相同的建筑物边缘，继而判定该局部建筑区域对应不同时相建筑物真实发生变化。

10.一种建筑物真伪变化的去除方法，其特征在于，包括以下步骤：

步骤一、通过变化检测得到变化区域二值图；通过如权利要求1-9所述的建筑物真伪变化的判定方法对建筑物的真伪变化进行判定；

步骤二、将步骤一中所得的变化区域二值图和建筑物的真伪变化判定结果进行二值图求交，去除变化区域二值图中的伪变化对象。