CN107341795B - 一种知识驱动的高空间分辨率遥感影像自动变化检测方法 - Google Patents

Abstract

本发明针对高空间分辨率影像变化检测的应用需求，公开了一种基于多元特征提取与知识驱动的高空间分辨率遥感影像自动变化检测方法。主要包括：S1知识驱动的地表覆盖区域分离；S2提取多元遥感影像特征；S3基于多元特征与地物分布知识的变化检测；S4基于形态学和连通域分析的变化检测后处理及矢量化。本方法能够在有效减少传统变化检测方法由于影像空间分辨率过高导致的虚警率，并保持较高的感兴趣变化地物检测精度。该方法无需人工干预，计算速度快，可满足海量卫星影像自动化生产的需求。

Description

一种知识驱动的高空间分辨率遥感影像自动变化检测方法

技术领域

本发明属于光学遥感图像处理技术领域，尤其涉及一种知识驱动的高空间分辨率遥感影像自动变化检测方法。

背景技术

遥感对地观测技术具有周期短、成本低、范围广等优点。尤其是随着遥感平台的快速发展，影像空间分辨率逐渐提高，地物细节信息更加丰富，利用遥感影像变化检测技术，能够为国土资源动态监测提供数据支撑。

由于空间分辨率的提升，相比于中低分辨率遥感影像，高空间分辨率影像变化检测技术存在以下难点：(1)现有变化检测方法大多基于影像光谱差异性进行分析，然而高空间分辨率影像光谱分辨率较低，“同物异谱-同谱异物”现象严重，导致差异图像中不同变化类型呈现相似的光谱特性，变化类型判别性能低；(2)影像获取时间不同，难以区分由光照、拍摄角度、季节等因素引起的不感兴趣的变化与由人为活动引起的感兴趣变化，变化检测精度低；(3)空间分辨率的提升导致检测结果存在大量椒盐噪声，影响变化检测精度。

发明内容

本发明的目的是针对高空间分辨率遥感影像变化信息提取，提供一种知识驱动的高空间分辨率遥感影像自动变化检测方法，将地物变化先验知识与遥感影像变化检测技术融合，在市郊等植被覆盖面积较大的区域，变化地物的光谱特性随空间分布一致性较高，光谱特征的差异能够反映地物的变化。而在城市覆盖区域，由于地物分布复杂，需要对感兴趣地物特征进行提取，在特征层分析地物变化，从而分离不感兴趣变化区域。

为达到上述目的，本发明包括如下步骤：

步骤1，基于可见光植被指数分离原始高空间分辨率遥感影像中的城市区域和非城市区域，得到二值的掩模图像；

步骤2，对原始高空间分辨率遥感影像进行多元特征提取，包括以下子步骤：

步骤2.1，利用分水岭方法对原始高空间分辨率遥感影像进行分割，并提取面向对象特征，得到面向对象特征影像；

步骤2.2，对原始高空间分辨率遥感影像提取形态学建筑物指数，得到形态学特征影像；

步骤3，分别对步骤2.1得到的面向对象特征影像和步骤2.2得到的形态学特征影像进行变化检测，得到两个特征下的两幅初始结果，然后利用步骤1的掩模对两幅初始结果进行融合，获取初始变化检测结果；

步骤4，基于形态学运算对初始变化检测结果进行后处理，获得形态学后处理的结果二值图，然后基于连通域分析对结果二值图进行矢量化，得到最终的变化检测矢量结果。

进一步的，所述步骤1的实现方式如下，

步骤1.1，对原始高分辨率影像逐像元计算可见光植被指数，获得VDVI影像，其公式为

其中ρ_g、ρ_r、ρ_b分别表示绿、红、蓝三个可见光波段的像元亮度值；

步骤1.2，计算VDVI影像均值，并利用均值将VDVI图像二值化，得到二值的掩模图像M；其中，二值化规则为，遍历VDVI图像中的每个像素，若该像素灰度值大于均值，则将其标记为1，表示该像素属于非城区；若该像素灰度值小于等于均值，则将其标记为0，表示该像素属于城区。

进一步的，所述步骤2.1的实现方式如下，

步骤2.1.1，提取原始影像的梯度影像GI，公式为

其中I_x为水平方向梯度，I_y为垂直方向梯度；

步骤2.1.2，对梯度影像统计灰度累计直方图，首先计算每个灰度值下的像素个数，然后对于灰度i，统计灰度范围在[0,i]之间的像素总数，得到灰度累计直方图，并将灰度累计直方图横轴映射到[0,100]；

步骤2.1.3，根据灰度累计直方图选择阈值G_T，将梯度影像中灰度值小于G_T的像素灰度值设置为G_T；

步骤2.1.4，令M₁,M₂,...,M_R表示图像GI的局部最小值点的坐标集合，C(M_i)为位于局部最小值相联系的汇水盆地内点的坐标集合，min和max代表梯度图像灰度的最小、最大值，T[n]表示灰度小于n的点的坐标集合，令

表示灰度级为n时汇水盆地被水淹没部分的合集，其中C_n(M_i)＝C(M_i)∩T[n]；

步骤2.1.5，开始时设定C[min+1]＝T[min+1]，按照以下规则迭代：令Q代表T[n]中连通分量的集合，对于每个连通分量q∈Q[n]，若q∩C[n-1]为空或q∩C[n-1]包含C[n-1]中的一个连通分量，则将q合并入Q[n-1]构成Q[n]；否则，在q∈Q[n]内构筑水坝；

步骤2.1.6，对于每个独立的分割区域，统计每个区域中的像素灰度均值，作为面向对象均值特征。

进一步的，所述步骤2.2的实现方式如下，

步骤2.2.1，对于每个像元x，计算波段最大值代表其亮度特性，其公式为

其中K为影像波段数，band_k(x)表示像元x在第k个波段的光谱值，b(x)表示像元x的亮度特征；

步骤2.2.2，构建MBI指数，利用基于顶帽变换与重构的差分形态学谱来提取建筑物的光谱-空间结构特性，具体包含以下步骤：

1)基于重构的顶帽变换，其公式为：

其中，

表示对亮度影像做重构开运算，s与d分别表示线性结构元素的长度和方向；

2)基于顶帽变换的形态学谱，其公式为

其中b表示步骤2.2.1中提取的亮度特征；

3)差分形态学谱，其公式为：

其中Δs为形态学谱间隔；

4)计算形态学建筑物指数，其公式为：

其中D和S表示形态学谱方向数和尺度数。

进一步的，所述步骤3的实现方式如下，

步骤3.1，对步骤2.1和步骤2.2得到的不同特征影像，分别利用变化向量分析法提取影像变化特征，公式为

其中m表示第m幅特征影像，

表示两个时相对应特征影像，b为影像第b波段，n是波段总数；

步骤3.2，针对每幅特征影像m，分别利用k-均值聚类方法将影像二值化，未变化像素标记为0，变化像素标记为1，得到不同特征下变化检测结果，其中进一步包括：

步骤3.2.1，从差值影像中CI^m随机取k个像素，作为k个簇的各自的中心；

步骤3.2.2，分别计算剩下的像素特征到k个簇中心的相异度，公式为：

将这些像素分别划归到相异度最低的簇，其中X_b、Y_b分别表示像素特征向量和聚类中心特征向量，n为特征向量维度；

步骤3.2.3，根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有像素各自维度的算术平均数；

步骤3.2.4，将CI^m中全部像素按照新的中心重新聚类；

步骤3.2.5，重复步骤3.2.2～3.2.4，直到聚类结果不再变化；

步骤3.3，利用步骤1中得到的掩模M将不同特征下变化检测结果融合，得到初始变化检测结果，公式为CM_init＝CM¹M+CM²(1-M)，其中CM¹、CM²分别代表步骤2.1、步骤2.2提取特征下由步骤3.2得到的变化检测中间结果。

进一步的，所述步骤4中基于形态学运算对初始变化检测结果进行后处理，获得形态学后处理的结果二值图的实现方式如下，

步骤4.1.1，对步骤3的初始变化检测结果进行腐蚀操作：逐像素扫描二值图像，用3×3的结构元素与其覆盖的二值图像做“与”操作，如果都为1，结果图像的该像素为1，否则为0，记经过腐蚀运算的二值图为CM_erode；

步骤4.1.2，对腐蚀后的二值图CM_erode进行开运算操作：先逐像素扫描CM_erode，用3×3的结构元素与其覆盖的CM_erode做“与”操作，如果都为1，结果图像的该像素为1，否则为0，记得到的中间结果为CM_m；再一次逐像素扫描CM_m，用3×3的结构元素与其覆盖的CM_m做“与”操作，如果都为0，结果图像的该像素为0，否则为1，记经过形态学后处理的结果二值图为CM_post。

进一步的，所述步骤4中基于连通域分析对结果二值图进行矢量化，得到最终的变化检测矢量结果的实现方式如下，

步骤4.2.1，基于递归连通域分析方法，搜索并标记CM_post中变化地物的连通域，其中变化地物象元指为1，逐像素遍历二值图CM_post，将未标记白色像素CM_post[i]作为种子，压入堆栈，标记连通域的label＝k，连通域的像素个数初始值N[k]＝0，标记种子像素mask[i]＝k，按顺序检查当前像素的四邻域像素是否存在白色未标记像素，如果存在，将该像素压入栈堆，记录连通域的像素个数N[k]＝N[k]+1；比较并记录当前连通域的最小矩形包围盒坐标(Xmin,Xmax,Ymin,Ymax)，依次从栈堆中取出顶部元素，以其作为新种子，递归搜索整幅图像，直至堆栈中无元素为止，即搜索得到一个连通域k；

步骤4.2.2，重复步骤4.2.1直到CM_post中的所有白色像素都被标记，即可得到所有连通域；

步骤4.2.3，根据应用需求中的最小图斑面积Smin和影像空间分辨率R设定连通域的像元数量下限Nmin＝Smin÷(R²)，遍历所有连通域，筛选并删除变化像元数量小于设定阈值Nmin的连通域；

步骤4.2.4，基于最小矩形包围盒方法，根据各个连通域坐标(Xmin,Xmax,Ymin,Ymax)提取矩形包围盒的4个顶点(Xmin,Ymin),(Xmax,Ymax)，(Xmin,Ymax),(Xmax,Ymin)，即可得到各个连通域的矢量边界，通过对矢量图进行二值化和再次矢量化，通过合并重叠的矢量得到最终的变化检测矢量结果CM_vector。

本发明具有如下优点和有益效果：

(1)能够很好的抑制高分辨率遥感影像复杂背景和成像条件造成的虚警；

(2)适用性广，可应用于不同分辨率高分辨率遥感影像；

(3)无需人工干预，可自动化处理；

(4)结果记录了变化地物的位置和范围，并且便于作业人员对图斑进行编辑和统计；

(5)计算速度快，可处理海量卫星影像，能满足实际生产的需求。

附图说明

图1为本发明中利用可见光植被指数提取掩模分离城区/非城区过程示意。其中，a为VDVI植被指数图像，为灰度图像；b是阈值分割后的掩模图像，其中白色表示非城区，黑色表示城区。

图2为高分影像多元特征提取示意图。其中a为面向对象均值特征；b为形态学建筑物指数特征。

图3为本发明中基于形态学的变化检测结果后处理的效果对比图。图中，a是腐蚀的示意图；b为开运算的示意图；c为原始二值图，白色代表变化区域，黑色代表未变化区域，红色圆圈代表虚警；d为经过腐蚀和开运算后处理的结果图。

图4为本发明中基于连通域分析对变化检测二值图进行矢量化的示意图。图中，a为变化检测结果图，白色为变化区域；b是递归连通域分析得到的各个连通域，并进行筛选合并后的最小矩形包围盒；c为包围盒矢量图。

具体实施方式

为了更好地理解本发明技术方案，下面将结合附图和实施例对本发明做进一步详细说明。

本具体实施方式ENVI/IDL环境下，由IDL结合C++语言进行开发与优化，整个过程可实现自动化处理。

步骤1，基于可见光植被指数分离原始影像城市区域/非城市区域。此步骤进一步包括：

步骤1.1，对原始高分辨率影像逐像元计算可见光植被指数，获得VDVI影像，其公式为：

其中ρ_g、ρ_r、ρ_b分别表示绿、红、蓝三个可见光波段的像元亮度值。

步骤1.2，计算VDVI影像均值，并利用均值将VDVI图像二值化，得到二值的掩模图像M。二值化规则为：遍历VDVI图像中的每个像素，若该像素灰度值大于均值，则将其标记为1，表示该像素属于非城区；若该像素灰度值小于等于均值，则将其标记为0，表示该像素属于城区。

步骤2，高空间分辨率遥感影像多元特征提取，本步骤进一步包括：

步骤2.1，利用分水岭方法对原始影像进行分割，并提取面向对象特征，此步骤进一步包括：

步骤2.1.1提取梯度影像，公式为：

其中I_x为水平方向梯度，I_y为垂直方向梯度；

步骤2.1.2对梯度影像统计灰度累计直方图。首先计算每个灰度值下的像素个数，然后对于灰度i，统计灰度范围在[0,i]之间的像素总数，得到灰度累计直方图，并将灰度累计直方图横轴映射到[0,100]；

步骤2.1.3根据灰度累计直方图选择阈值G_T，将梯度影像中灰度值小于G_T的像素灰度值设置为G_T；

步骤2.1.4令M₁,M₂,...,M_R表示图像GI的局部最小值点的坐标集合，C(M_i)为位于局部最小值相联系的汇水盆地内点的坐标集合，min和max代表梯度图像灰度的最小、最大值，T[n]表示灰度小于n的点的坐标集合。令

表示灰度级为n时汇水盆地被水淹没部分的合集，其中C_n(M_i)＝C(M_i)∩T[n]；

步骤2.1.5开始时设定C[min+1]＝T[min+1]，按照以下规则迭代：令Q代表T[n]中连通分量的集合，对于每个连通分量q∈Q[n]，若q∩C[n-1]为空或q∩C[n-1]包含C[n-1]中的一个连通分量，则将q合并入Q[n-1]构成Q[n]；否则，需要在q∈Q[n]内构筑水坝。

步骤2.1.6对于每个独立的分割区域，统计每个区域中的像素灰度均值，作为面向对象均值特征。

步骤2.2，对原始影像提取形态学建筑物指数。此步骤进一步包括：

步骤2.2.1，计算亮度。对于每个像元x，计算波段最大值代表其亮度特性，其公式为：

其中K为影像波段数，在本实施例中K为3，band_k(x)表示像元x在第k个波段的光谱值，b(x)表示像元x的亮度特征。

步骤2.2.2，构建MBI指数。利用基于顶帽变换与重构的差分形态学谱来提取建筑物的光谱-空间结构特性。具体包含以下步骤：

1)基于重构的顶帽变换(W-TH)，其公式为：

其中，

表示对亮度影像做重构开运算，s与d分别表示线性结构元素(SE)的长度和方向。

2)基于顶帽变换的形态学谱，其公式为：

其中b表示步骤2.21中提取的亮度特征；

3)差分形态学谱，其公式为：

其中Δs为形态学谱间隔，参考值为10；

4)计算形态学建筑物指数，其公式为：

其中D和S表示形态学谱方向数和尺度数，其值通过人为设定，参考值分别为4和10。

步骤3，分别对前两个步骤得到的面向对象特征影像和形态学特征影像进行变化检测，得到两个特征下的两幅初始结果，然后利用步骤1的掩模进行融合，获取初始变化检测结果，本步骤进一步包括：

步骤3.1，对上述两个步骤得到的不同特征影像，分别利用变化向量分析法提取影像变化特征，公式为：

其中m表示第m幅特征影像，本实例中m＝2，

表示两个时相对应特征影像，b为影像第b波段，n是波段总数。

步骤3.2，针对每幅特征影像m，分别利用k-均值聚类方法变化将影像二值化，未变化像素标记为0，变化像素标记为1，得到不同特征下变化检测结果，其中进一步包括：

步骤3.2.1，从差值影像中CI^m随机取k个像素，作为k个簇的各自的中心；

步骤3.2.2，分别计算剩下的像素特征到k个簇中心的相异度，公式为：

将这些像素分别划归到相异度最低的簇，其中X_b、Y_b分别表示像素特征向量和聚类中心特征向量，n为特征向量维度；

步骤3.2.3，根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有像素各自维度的算术平均数；

步骤3.2.4，将CI^m中全部像素按照新的中心重新聚类；

步骤3.2.5，重复步骤3.2.2、3.2.3、3.2.4，直到聚类结果不再变化。

步骤3.3，利用步骤1.2中得到的掩模M将不同特征下变化检测结果融合，得到初始变化检测结果。公式为：CM_init＝CM¹M+CM²(1-M)，其中CM¹、CM²分别代表步骤2.1、步骤2.2提取特征下由步骤3.2得到的变化检测中间结果。

步骤4，基于形态学和连通域分析的变化检测后处理及矢量化。此步骤进一步包括：

步骤4.1，基于形态学运算对检测结果进行后处理，以消除高分辨率遥感影像中复杂地物细节造成的虚警。此步骤进一步包括：

步骤4.1.1，对步骤3的初始变化检测结果(二值图)进行腐蚀操作。逐像素扫描二值图像，用3×3的结构元素与其覆盖的二值图像做“与”操作，如果都为1，结果图像的该像素为1，否则为0。记经过腐蚀运算的二值图为CM_erode；

步骤4.1.2，对腐蚀后的二值图CM_erode进行开运算操作。形态学的开运算即先腐蚀后膨胀，逐像素扫描CM_erode，用3×3的结构元素与其覆盖的CM_erode做“与”操作，如果都为1，结果图像的该像素为1，否则为0，记得到的中间结果为CM_m。再一次逐像素扫描CM_m，用3×3的结构元素与其覆盖的CM_m做“与”操作，如果都为0，结果图像的该像素为0，否则为1，记经过形态学后处理的结果二值图为CM_post。

下一步骤在形态学后处理得到的变化检测二值图CM_post上进行。

步骤4.2，统计变化像元的个数和分布情况，基于连通域分析对CM_post进行矢量化。此步骤进一步包括：

步骤4.2.1，基于递归连通域分析方法，搜索并标记CM_post中变化地物(像元值为1)的连通域。逐像素遍历二值图CM_post，将未标记白色像素CM_post[i]作为种子，压入堆栈。标记连通域的label＝k，连通域的像素个数初始值N[k]＝0。标记种子像素mask[i]＝k，按顺序检查当前像素的四邻域像素(左、右、上、下)是否存在白色未标记像素，如果存在，将该像素压入栈堆，记录连通域的像素个数N[k]＝N[k]+1。比较并记录当前连通域的最小矩形包围盒坐标(Xmin,Xmax,Ymin,Ymax)。依次从栈堆中取出顶部元素，以其作为新种子，递归搜索整幅图像，直至堆栈中无元素为止，即搜索得到一个连通域k。

步骤4.2.2，重复步骤4.2.1直到的CM_post中的所有白色像素都被标记，即可得到所有连通域；

步骤4.2.3，根据应用需求中的最小图斑面积Smin和影像空间分辨率R设定连通域的像元数量下限Nmin＝Smin÷(R²)，遍历所有连通域，筛选并删除变化像元数量小于设定阈值Nmin的连通域；

步骤4.2.4，基于最小矩形包围盒方法，根据各个连通域坐标(Xmin,Xmax,Ymin,Ymax)提取矩形包围盒的4个顶点(Xmin,Ymin),(Xmax,Ymax)，(Xmin,Ymax),(Xmax,Ymin)，即可得到各个连通域的矢量边界，通过对矢量图进行二值化和再次矢量化，即可合并重叠的矢量得到最终的变化检测矢量结果CM_vector。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于，包括以下步骤：

步骤1，基于可见光植被指数分离原始高空间分辨率遥感影像中的城市区域和非城市区域，得到二值的掩模图像；

步骤2，对原始高空间分辨率遥感影像进行多元特征提取，包括以下子步骤：

步骤2.1，利用分水岭方法对原始高空间分辨率遥感影像进行分割，并提取面向对象特征，得到面向对象特征影像；所述步骤2.1的实现方式如下，

步骤2.1.1，提取原始影像的梯度影像GI，公式为

其中I_x为水平方向梯度，I_y为垂直方向梯度；

步骤2.1.2，对梯度影像统计灰度累计直方图，首先计算每个灰度值下的像素个数，然后对于灰度i，统计灰度范围在[0,i]之间的像素总数，得到灰度累计直方图，并将灰度累计直方图横轴映射到[0,100]；

步骤2.1.3，根据灰度累计直方图选择阈值G_T，将梯度影像中灰度值小于G_T的像素灰度值设置为G_T；

步骤2.1.4，令M₁,M₂,...,M_R表示图像GI的局部最小值点的坐标集合，C(M_i)为位于局部最小值相联系的汇水盆地内点的坐标集合，min和max代表梯度图像灰度的最小、最大值，T[n]表示灰度小于n的点的坐标集合，令

表示灰度级为n时汇水盆地被水淹没部分的合集，其中C_n(M_i)＝C(M_i)∩T[n]；

步骤2.1.5，开始时设定C[min+1]＝T[min+1]，按照以下规则迭代：令Q代表T[n]中连通分量的集合，对于每个连通分量q∈Q[n]，若q∩C[n-1]为空或q∩C[n-1]包含C[n-1]中的一个连通分量，则将q合并入Q[n-1]构成Q[n]；否则，在q∈Q[n]内构筑水坝；

步骤2.1.6，对于每个独立的分割区域，统计每个区域中的像素灰度均值，作为面向对象均值特征；

步骤2.2，对原始高空间分辨率遥感影像提取形态学建筑物指数，得到形态学特征影像；

步骤3，分别对步骤2.1得到的面向对象特征影像和步骤2.2得到的形态学特征影像进行变化检测，得到两个特征下的两幅初始结果，然后利用步骤1的掩模对两幅初始结果进行融合，获取初始变化检测结果；

步骤4，基于形态学运算对初始变化检测结果进行后处理，获得形态学后处理的结果二值图，然后基于连通域分析对结果二值图进行矢量化，得到最终的变化检测矢量结果。

2.如权利要求1所述的一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于：所述步骤1的实现方式如下，

步骤1.1，对原始高分辨率影像逐像元计算可见光植被指数，获得VDVI影像，其公式为

其中ρ_g、ρ_r、ρ_b分别表示绿、红、蓝三个可见光波段的像元亮度值；

步骤1.2，计算VDVI影像均值，并利用均值将VDVI图像二值化，得到二值的掩模图像M；其中，二值化规则为，遍历VDVI图像中的每个像素，若该像素灰度值大于均值，则将其标记为1，表示该像素属于非城区；若该像素灰度值小于等于均值，则将其标记为0，表示该像素属于城区。

3.如权利要求2所述的一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于：所述步骤2.2的实现方式如下，

步骤2.2.1，对于每个像元x，计算波段最大值代表其亮度特性，其公式为

其中K为影像波段数，band_k(x)表示像元x在第k个波段的光谱值，b(x)表示像元x的亮度特征；

步骤2.2.2，构建MBI指数，利用基于顶帽变换与重构的差分形态学谱来提取建筑物的光谱-空间结构特性，具体包含以下步骤：

1)基于重构的顶帽变换，其公式为：

其中，

表示对亮度影像做重构开运算，s与d分别表示线性结构元素的长度和方向；

2)基于顶帽变换的形态学谱，其公式为

其中b表示步骤2.2.1中提取的亮度特征；

3)差分形态学谱，其公式为：

其中Δs为形态学谱间隔；

4)计算形态学建筑物指数，其公式为：

其中D和S表示形态学谱方向数和尺度数。

4.如权利要求3所述的一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于：所述步骤3的实现方式如下，

步骤3.1，对步骤2.1和步骤2.2得到的不同特征影像，分别利用变化向量分析法提取影像变化特征，公式为

其中m表示第m幅特征影像，

表示两个时相对应特征影像，b为影像第b波段，n是波段总数；

步骤3.2，针对每幅特征影像m，分别利用k-均值聚类方法将影像二值化，未变化像素标记为0，变化像素标记为1，得到不同特征下变化检测结果，其中进一步包括：

步骤3.2.1，从差值影像中CI^m随机取k个像素，作为k个簇的各自的中心；

步骤3.2.2，分别计算剩下的像素特征到k个簇中心的相异度，公式为：

将这些像素分别划归到相异度最低的簇，其中X_b、Y_b分别表示像素特征向量和聚类中心特征向量，n为特征向量维度；

步骤3.2.3，根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有像素各自维度的算术平均数；

步骤3.2.4，将CI^m中全部像素按照新的中心重新聚类；

步骤3.2.5，重复步骤3.2.2～3.2.4，直到聚类结果不再变化；

步骤3.3，利用步骤1中得到的掩模M将不同特征下变化检测结果融合，得到初始变化检测结果，公式为CM_init＝CM¹M+CM²(1-M)，其中CM¹、CM²分别代表步骤2.1、步骤2.2提取特征下由步骤3.2得到的变化检测中间结果。

5.如权利要求4所述的一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于：所述步骤4中基于形态学运算对初始变化检测结果进行后处理，获得形态学后处理的结果二值图的实现方式如下，

步骤4.1.1，对步骤3的初始变化检测结果进行腐蚀操作：逐像素扫描二值图像，用3×3的结构元素与其覆盖的二值图像做“与”操作，如果都为1，结果图像的该像素为1，否则为0，记经过腐蚀运算的二值图为CM_erode；

步骤4.1.2，对腐蚀后的二值图CM_erode进行开运算操作：先逐像素扫描CM_erode，用3×3的结构元素与其覆盖的CM_erode做“与”操作，如果都为1，结果图像的该像素为1，否则为0，记得到的中间结果为CM_m；再一次逐像素扫描CM_m，用3×3的结构元素与其覆盖的CM_m做“与”操作，如果都为0，结果图像的该像素为0，否则为1，记经过形态学后处理的结果二值图为CM_post。

6.如权利要求5所述的一种知识驱动的高空间分辨率遥感影像自动变化检测方法，其特征在于：所述步骤4中基于连通域分析对结果二值图进行矢量化，得到最终的变化检测矢量结果的实现方式如下，

步骤4.2.1，基于递归连通域分析方法，搜索并标记CM_post中变化地物的连通域，其中变化地物象元指为1，逐像素遍历二值图CM_post，将未标记白色像素CM_post[i]作为种子，压入堆栈，标记连通域的label＝k，连通域的像素个数初始值N[k]＝0，标记种子像素mask[i]＝k，按顺序检查当前像素的四邻域像素是否存在白色未标记像素，如果存在，将该像素压入栈堆，记录连通域的像素个数N[k]＝N[k]+1；比较并记录当前连通域的最小矩形包围盒坐标(Xmin,Xmax,Ymin,Ymax)，依次从栈堆中取出顶部元素，以其作为新种子，递归搜索整幅图像，直至堆栈中无元素为止，即搜索得到一个连通域k；

步骤4.2.2，重复步骤4.2.1直到CM_post中的所有白色像素都被标记，即可得到所有连通域；

步骤4.2.3，根据应用需求中的最小图斑面积Smin和影像空间分辨率R设定连通域的像元数量下限Nmin＝Smin÷(R2)，遍历所有连通域，筛选并删除变化像元数量小于设定阈值Nmin的连通域；

步骤4.2.4，基于最小矩形包围盒方法，根据各个连通域坐标(Xmin,Xmax,Ymin,Ymax)提取矩形包围盒的4个顶点(Xmin,Ymin),(Xmax,Ymax)，(Xmin,Ymax),(Xmax,Ymin)，即可得到各个连通域的矢量边界，通过对矢量图进行二值化和再次矢量化，通过合并重叠的矢量得到最终的变化检测矢量结果CM_vector。

Images