CN104751478A - 一种基于多特征融合的面向对象的建筑物变化检测方法 - Google Patents

Abstract

本发明公开了一种基于多特征融合的面向对象的建筑物变化检测方法，首先求得图像像素点的形态学建筑指数（MBI），纹理特征和慢特征分析图（SFA）；利用MBI指数和纹理特征进行FNEA分割；然后求出每个对象的三个特征值，再作差，并利用K均值聚类算法求阈值，得到特征变化图；再利用AC指数进行后处理；利用熵值法对不同的特征变化图求权重，按照权重设定阈值求得变化图像；最后利用投票法进行后处理。得到变化检测结果。本发明利用MBI，SFA和纹理特征来进行建筑物的变化检测；把MBI和纹理特征值加入到FNEA分割中；提出AC指数和投票法进行后处。本发明可为高分辨率遥感图像应用与土地覆盖和城市扩张提供一种新途径。

Description

一种基于多特征融合的面向对象的建筑物变化检测方法

技术领域

本发明属于遥感图像数据的信息提取技术领域，具体涉及一种基于多特征融合的面向对象的建筑物变化检测方法。

背景技术

随着经济发展和城市化进程的加快，人类的各种生产建设活动正在日益改变着城市及其周边的自然环境和土地覆盖类型。因此快速有效的监测这些变化信息，分析变化原因和影响结果，对我国的可持续发展具有重要意义。遥感图像变化检测技术就是根据同一地区不同时相的两幅遥感图像，来获取地表地物变化的一门技术，该技术的快速发展为更新地理空间数据发挥着重大的作用。遥感变化检测技术是数字图像处理方法、计算机视觉技术和人工智能、模式识别理论的综合应用。

遥感图像变化检测包括三个层次：像元级变化检测、特征级变化检测和目标级变化检测。变化检测的三个层次各有优缺点，在具体的变化检测过程中要检测到哪个层次是根据任务的需要来确定的。像元级变化检测保持了尽可能多的原始信息，具有特征级和目标级层次上所不具备的细节信息，但像元级变化检测仅考虑像素属性的变化，而未考虑其空间等特征属性的变化；特征级变化检测不仅考虑到空间形状的变化，而且还考虑特征属性的变化，但特征级变化检测依赖于特征提取的结果，况且单独使用某一特征进行变化检测有可能造成很大的漏检和错检；目标级变化检测最大的优点是它接近于用户的需求，检测的结果可直接应用，它的不足之处在于目标提取的困难性。

随着卫星技术的快速发展，使得高分辨率遥感图像的变化检测成为可能。分辨率在10m以内的高分辨率遥感图像(Very High Resolution Imagery，VHR)已经广泛应用于社会经济的很多领域，中低分辨率遥感图像中的点目标在高分辨率图像中变成了面，图像包含更多、更丰富的地理和地形信息。高效利用高分辨率遥感图像丰富的地物细节变化信息、像元之间以及像元属性之间的相互关系，能有效地抑制自然地物变化和不同成像条件引起的信息干扰，因此对高分辨率遥感图像的变化检测研究具有重要的理论意义和现实意义。传统的遥感图像变化检测是基于像素的，如图像差值法、图像比值法。图像差值法和比值法对图像的质量和预处理要求相对较高，不可避免的几何配准误差、相对辐射校正精度、阴影等都是基于差异图像的高分辨率遥感图像变化检测中的典型问题。简单的单波段相减由于没有考虑波段之间的统计相关性，使得绝对值不同的数值相减得到同样大小的差值，不同像素灰度值之间相比得到同样大小的比值，忽略了不同地物在不同敏感波段存在的差异，导致潜在可利用信息的丢失。

在这背景下，面向对象的变化检测方法能够将像元-像元之间的差异推广到对象-对象，把传统的以像素为单位的变化检测推广到以对象为单位的变化检测，最后提取变化/未变化信息。

发明内容

为了解决上述的技术问题，本发明提出了一种基于多特征融合的面向对象的建筑物变化检测方法。

本发明所采用的技术方案是：一种基于多特征融合的面向对象的建筑物变化检测方法，其特征在于，包括以下步骤：

步骤1：对所选区域的两幅不同时相的高分辨率遥感图像A和B进行预处理，主要包括几何纠正、辐射纠正、几何配准和图像裁剪；

步骤2：计算A和B两幅图像的每个像素点的形态学建筑指数(MBI)、固定窗口的纹理特征和慢特征分析(SFA)，得到图像A和B的MBI特征图像、纹理特征图像和SFA特征图像；

步骤3：选取步骤1中的一幅图像采用基于改进的FNEA多尺度分割方法进行分割，得到一幅多尺度分割图像，并利用步骤2的计算结果来改进多尺度分割算法进行分割，分割的结果是得到一个个对象；并由分割的结果得到每个对象对应像素点的坐标或索引，称为索引矩阵；然后按照此索引矩阵对步骤1中的另一幅图像进行分割，分割的结果使两幅图像具有同样大小的对象；

步骤4：因步骤2中得到的MBI特征图像、纹理特征图像和SFA特征图像具有尺度不统一的特点，故采用单位标准差归一化方法对不同的特征图像进行特征优化；

步骤5：对步骤4归一化后的结果，求每个对象的特征均值，以得到各个对象的特征图像；

步骤6：对每个对象求其在不同时相的特征图像的差值，利用k均值聚类算法对三幅不同的差值图像求其阈值，使其阈值自动化，以此得到三幅不同特征的变化图像；

步骤7：对步骤6中得到的三幅不同特征的变化图像进行AC指数后处理；

步骤8：利用熵值法，对步骤7中得到的三幅不同特征的变化图像加不同的权重，再设定阈值，以得到变化检测结果；

步骤9：利用基于投票法的多尺度融合方法对变化检测的结果进行处理，以得到更高的检测精度；

步骤10：对步骤9得到的结果进行精度评定。

作为优选，步骤1中所述的几何纠正采用基于多项式的遥感图像几何纠正，控制点选取分布均匀，重采样采用双线性内插法，最后得到误差要求标准为RMSE<0.5像素；所述的辐射纠正方法采用的是相对辐射归一化纠正。

作为优选，步骤3中所述的选取步骤1中的一幅图像进行分割，选取方法为：当两幅图像分辨率不一致时，用空间分辨率高的图像来进行分割；当分辨率相同时，按照获取时相的时间顺序，选取后一时期的图像进行分割。

作为优选，步骤3中所述的基于改进的FNEA多尺度分割方法，其具体实现过程为：从一个像元起步，先将单像元合并为较小的对象，然后把具有异质性最小的较小对象合并成较大的对象，这样不断合并，直到判断条件不成立，合并操作就终止，最终分割的结果中所有图像对象的平均异质性最小；

在判断两相邻区域是否能够合并时，用总异质性值和先前设定好的尺度阈值进行比较，如果小于尺度阈值就合并，否则就结束合并操作；

总的异质性h计算公式为：

h＝w_spectral*h_spectral+w_shape*h_shape+w_MBI*h_MBI+w_texture*h_texture；

其中，w_spectral,w_shape,w_MBI,w_texture分别为光谱异质性h_spectral、形状异质性h_shape、MBI异质性h_MBI和纹理异质性h_texture对应的权重；MBI异质性和纹理异质性的计算方法和光谱异质性的计算方法原理相同，即计算每一波段的标准差与该波段权重的乘积，再把个波段的值进行累加。

作为优选，步骤6中所述的利用k均值聚类算法对三幅不同的差值图像求其阈值，其具体实现包括以下子步骤：

步骤6.1：从数据集中随机选取K个图像单元作为初始聚类中心；

步骤6.2：计算各个图像单元到聚类中心的光谱距离，将它们一一归类到最近的那个聚类中心所在的类；

步骤6.3：计算新形成的每一个聚类的图像单元的光谱均值，从而得到新的聚类中心；

步骤6.4：迭代实施步骤6.2和步骤6.3，直至前后两次的聚类中心没有任何变化，说明聚类调整结束，聚类准则函数已经收敛；

当输入图像分别为归一化后的MBI特征图像、SFA特征图像和纹理特征图像时，用K均值二值聚类所得的结果，即将这些特征分为两个分离度最高的图像类别，即低相似性与高相似性类别；因此，由此得到的两个图像类别可分别对应于变化区域与未变化区域，即可得到三个阈值X，Y，Z。

作为优选，步骤7中所述的AC指数后处理，具体的公式如下：

AC＝a*(area)/circle  ；

其中，circle是对每一个检测为变化的对象进行圆形拟合，该圆形包含了该对象的所有被检测为变化建筑物的像素以及部分未变化的像素，area为圆形中该对象所有被检测为变化建筑物的像素，a是用来调整该比值的大小。

作为优选，步骤8的具体实现包括以下子步骤：

步骤8.1：原始数据的收集与处理

设有m个待评项目，n个评价指标，形成原始数据矩阵x＝{x_ij}_mn：

$X=\left[\begin{array}{cccc}{X}_{11}& X_{12}& ...& X_{1n}\\ X_{21}& X_{22}& ...& X_{2n}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ X_{m1}& X_{m2}& ...& X_{\mathrm{mn}}\end{array}\right];$

其中x_ij表示第j个指标下第i个项目的评价指标的数值，m为各个特征变化图像值，m＝3，x_ij为第i个特征图像对应的第j个像素的特征值；

步骤8.2：计算第j个指标下第i个项目的指标的比重P_ij：

$p_{\mathrm{ij}}=\frac{X_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{X}_{\mathrm{ij}}}(i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n);$

步骤8.3：计算第j个指标的熵值e_j：

$e_j=-k\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\&CenterDot;\ln \left(p_{\mathrm{ij}}\right),$ 其中， $k=\frac{1}{\ln \left(m\right)},e_j\&GreaterEqual;0;$

步骤8.4：计算第j个指标的熵权w_j：

$w_j=\frac{(1-e_j)}{m-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{e}_j},(0\&le;w_j\&le;1,\underset{1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_j=1);$

即由上式步骤得到的w_j即为各个特征变化图像的权重值。

作为优选，步骤9所述的利用基于投票法的多尺度融合方法对变化检测的结果进行处理，其具体实现过程为：当某像元在不同尺度下的变化结果中，被检测为变化的次数多于用户定义的投票阈值时，则判定该像元在融合结果中为变化，否则，则视其为未变化，具体公式表示如下：

$M_i=\left\{\begin{array}{cc}1& \begin{array}{cc}\mathrm{if}& n_{\mathrm{ic}}>n_v\end{array}\\ 0& \mathrm{otherwise}\end{array}\right.;$

其中，M_i表示像元i的变化检测融合结果，n_ic为该像元在各尺度变化检测结果中被判定为变化次数，n_v表示用户定义的变化次数投票阈值。

基于变化检测的复杂性，现在大多数变化检测都是针对特定目标的变化检测，如建筑物，道路，森林等。本发明属于建筑物变化检测范畴，并且是基于多特征融合的面向对象的变化检测。基于特征级变化检测主要是利用某类地物的特征区别于其他地物的特征来进行变化检测，特征在图像上表现一般较稳定，受辐射差异影响较小，且不易受遥感图像时相变化的影响。如果仅仅利用光谱特征来进行建筑物的变化检测会存在很大的错检，特别是对不同传感器的图像来进行变化检测。因为建筑物的纹理特征比较规则，所以利用纹理特征来进行变化检测会得到较好的效果。所以本发明提出一种特定目标的变化检测，并且是基于多特征的面向对象的建筑物的变化检测，与传统的变化检测方法相比，本发明利用基于改进的FNEA多尺度分割算法对图像进行分割，对分割的结果进行特征提取，提取的特征包括形态学建筑指数(MBI)，纹理特征和慢特征分析(SFA)，SFA主要是使不变化或者变化很小的像素与变化的像素分开的更加明显。由于特征之间存在尺度不统一等问题，本发明对不同的特征图像还进行了特征优化。此外，本发明提出利用k均值聚类算法对得到的不同的特征差值图像分别计算其阈值，使其阈值自动化，以得到不同的二值图像，得到初步的变化检测检测结果。由于建筑物与道路在形状等各方面具有极大地相似性，会对建筑物的检测结果产生极大的错检，所以本发明原创性地提出了AC指数对得到的变化图像进行后处理，以消除道路对建筑物的影响。最后，本发明利用慢特征分析图像，MBI特征图像和纹理特征图像进行多特征融合，根据熵值法对每一副图像加不同的权值，得到更好的变化检测结果。由于面向对象的变化检测的精度在很大程度上是由于分割结果造成的影响，所以本发明在此利用多尺度融合的投票法来减少分割的结果对变化检测精度的影响。

本发明提供的技术方案的有益效果为：原创性地提出加入形态学建筑指数和纹理特征来进行FNEA多尺度分割，相比于传统的分割算法，该方法能够对建筑物得到更好的分割结果；利用形态学建筑指数，慢特征分析以及纹理特征对建筑物进行多特征变化检测能得到更好的检测效果；为了使阈值的自适应性，本发明提出利用k均值聚类算法对不同的特征差值图像求得不同的阈值(X，Y，Z)；由于道路与建筑物在某些方面具有很大的相似性，因此本发明原创性地提出AC指数来减少道路对建筑物的影响；最后本发明将三幅建筑物特征图像按照熵值法赋予不同的权重，通过设定的阈值，来得到变化检测的结果；为了减少分割算法对变化检测精度的影响，本发明提出利用基于投票法的多尺度融合方法来提高变化检测的精度。因此，本发明为建筑物的变化检测提供了一种新途径。

附图说明

附图1：为本发明实施例的流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种基于多特征融合的面向对象的建筑物变化检测方法，包括以下步骤：

步骤1：对所选区域的两幅不同时相的高分辨率遥感图像A和B进行预处理，主要包括几何纠正、辐射纠正、几何配准和图像裁剪；

(1)几何纠正：本发明几何纠正的方法采用基于多项式的遥感图像几何纠正，控制点选取分布均匀，重采样采用双线性内插法，最后得到误差要求标准为RMSE<0.5像素；

(2)辐射纠正：本发明采用的是相对辐射归一化纠正；

(3)几何配准：为了防止因两幅图像对应位置不一致而引起的检测误差，本发明对两幅图像进行了几何配准，同名点的获取用图像自动配准的方法；

(4)图像裁剪：获得的图像不一定都是感兴趣的图像，所以本发明对获得的图像进行裁剪，提取出感兴趣区域。

步骤2：计算A和B两幅图像的每个像素点的形态学建筑指数(MBI)、固定窗口的纹理特征和慢特征分析(SFA)，得到图像A和B的MBI特征图像、纹理特征图像和SFA特征图像；

(一)形态学建筑指数(MBI)是近年来建筑物提取研究方向的重要成果之一，已被证明能有效实现高分辨率遥感图像建筑物提取。根据建筑物屋顶通常在其邻域范围内表现出较亮的光谱特征，该指数利用多尺度、多方向的形态学操作来描述建筑物的光谱、结构特征，并利用一系列线性结构算子(Structure Element，SE)对高帽变换结果进行差分形态学轮廓(DifferentialMorphological Profile，DMP)重建，具体实现方式如下：

(1)计算亮度值：首先，提取高分辨率图像中各像元在不同波段上的最大值，并将其作为亮度图像：

B(x,y)＝max_1≤k≤K(band_k(x,y))  (0)其中，band_k(x,y)表示像元(x,y)在第k波段上的光谱值，K为图像的波段数；

(2)基于重建的白高帽变换：利用线性结构元素(SE)对图像B(x,y)进行形态学开操作(腐蚀和膨胀)，得到的结果再进行重建。基于开重建操作的高帽变换(THR)定义为：

${\mathrm{THR}}^S\left(b\right)=b-{\mathrm{\&gamma;}}_{\mathrm{RE}}^S\left(b\right)---\left(0\right)$

其中：表示对亮度图像B进行开重建操作s代表所选线性SE的尺度；

(3)多方向的THR：单一的线性结构元素不能包含多方向信息，所以提出利用多方向的线性结构元素来计算THR，具体地计算公式如下：

$\stackrel{\&OverBar;}{{\mathrm{THR}}^S\left(b\right)}=\underset{\mathrm{dir}}{\mathrm{mean}}\left({\mathrm{THR}}^{S.\mathrm{dir}}\left(b\right)\right)---\left(0\right)$

其中，为当线性结构算子SE的尺度s唯一时不同方向的THR的均值；

(4)多尺度的THR：由于在高分辨率遥感图像中，建筑物具有不同的形状、大小、高度和面积，因此需要计算多尺度的THR，具体的：

${\mathrm{THR}}_{\mathrm{DMP}}=\{{\mathrm{THR}}_{\mathrm{DMP}}^{S^{\min }},\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{THR}}_{\mathrm{DMP}}^S,\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{THR}}_{\mathrm{DMP}}^{S^{\max }}\}---\left(0\right)$

${\mathrm{THR}}_{\mathrm{DMP}}^S=|\stackrel{\&OverBar;}{{\mathrm{THR}}^{S+\mathrm{\&Delta;S}}\left(b\right)}-\stackrel{\&OverBar;}{{\mathrm{THR}}^S\left(b\right)}|---\left(0\right)$

S^min≤S≤S^max  (0)

其中，△S为SE的尺度增长步长，且满足(6)式；

(5)形态学建筑指数(MBI)：由于建筑物在尺度大小与方向上相对于其他地物类别更具多样性，因此，定义形态学建筑指数(MBI)为在不同尺度不同方向上对白高帽变换结果进行求均值结果：

$\mathrm{MBI}=\underset{s}{\mathrm{mean}}\left({\mathrm{THR}}_{\mathrm{DMP}}\right)---\left(0\right)$

计算得到的MBI值越大，是建筑物的可能性就越大，经实验证明，形态学建筑指数(MBI)能够很好的提取出建筑物。

(二)计算图像中以各个像素点为中心的3*3的小窗口矩阵的纹理特征，以其中的某一个纹理特征值来代替该点的像素灰度值，最后求得整副图像的像素值是用对应的纹理值来代替，即得到纹理特征图像。

(三)慢特征分析是用来提取输入两幅图像中的不变化特征或者慢变化特征的一种算法，是把两幅不同时相的图像转化到新的特征空间中，在此空间中，不变化的像素和变化的像素将被分开的更加明显。实验证明，将慢特征分析用于建筑物的变化检测能得到较好的检测效果，并且SFA算法得到的是全局最优解而不是局部最优解。具体地步骤为:

(1)给定一个双时相图像的光谱向量和i为像素的个数，N为图像的维数。首先我们要对输入的图像进行归一化：

$\begin{array}{ccc}\widehat{x_j^i}=\frac{x_j^i-u_{x_j}}{{\mathrm{\&delta;}}_{x_j}}& \mathrm{and}& \widehat{y_j^i}=\frac{y_j^i-u_{y_j}}{{\mathrm{\&delta;}}_{y_j}}---\left(0\right)\end{array}$

其中u_jδ_j分别为图像X第j波段的均值和方差，这个过程类似于相对辐射校正归一化；

(2)慢特征分析是找到一个函数g(x)使得归一化后的图像满足如下条件：

$\begin{array}{cccc}\frac{1}{p}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{(g_j\left(\widehat{x^i}\right)-g_j\left(\widehat{y^i}\right))}^2& \mathrm{is}& \min & \mathrm{imal}\end{array}---\left(0\right)$

并且在如下的限制条件下：

$\frac{1}{2p}[\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j\left(\widehat{x^i}\right)+\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j\left(\widehat{g^i}\right)]=0---\left(0\right)$

$\frac{1}{2p}[\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j{\left(\widehat{x^i}\right)}^2+\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j{\left(\widehat{g^i}\right)}^2]=1---\left(0\right)$

$\&ForAll;1<:\frac{1}{2p}[\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j\left(\widehat{x^i}\right)g_1\left(\widehat{x^i}\right)+\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{g}_j\left(\widehat{y^i}\right)g_1\left(\widehat{y^i}\right)]=0---\left(0\right)$

(3)第(2)步可转化为求A，B的广义特征向量和特征值AW＝BWΛ，其中

$A=\frac{1}{p}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(\widehat{x^i}-\widehat{y^i}){(\widehat{x^i}-\widehat{y^i})}^T=\mathrm{\&Sigma;\&Delta;}---\left(0\right)$

$B=\frac{1}{2p}[\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left(\widehat{x^i}\right){\left(\widehat{x^i}\right)}^T+\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left(\widehat{y^i}\right){\left(\widehat{y^i}\right)}^T]=\frac{1}{2}({\mathrm{\&Sigma;}}_x+{\mathrm{\&Sigma;}}_y)---\left(0\right)$

(4)对由A，B求得的广义特征值按从小到大的顺序进行排列，并得到对应的新的特征向量W，则g(x)就定义为：

g(x)＝w^T*x  (0)

(5)慢特征分析即定义为：

${\mathrm{SFA}}_j=w_j^T*\hat{x}-w_j^T*\hat{y}---\left(0\right)$

步骤3：选取步骤1中的一幅图像采用基于改进的FNEA多尺度分割方法进行分割，得到一幅多尺度分割图像，并利用步骤2的计算结果来改进多尺度分割算法进行分割，分割的结果是得到一个个对象；并由分割的结果得到每个对象对应像素点的坐标或索引，称为索引矩阵；然后按照此索引矩阵对步骤1中的另一幅图像进行分割，分割的结果使两幅图像具有同样大小的对象；

图像分割的方式总共有三种；(1)单独分割方式：对两幅图像都进行分割，该方法最常用于面向对象的变化检测，且只能用于分类后处理的方法中，因此变化检测的精度只能归结于分割算法和分类精度的好坏；(2)多源图像波段复合的分割方式：对不同时相的图像进行多波段复合，然后再进行分割，当多时相图像是由不同高分辨率遥感卫星获取时，该方式将不再适用。其原因在于，不同卫星传感器由于其接收的光谱区间范围不同，所获取的图像波段的光谱位置及带宽都存在差异。因此，在不同光谱区间上的灰度图像不具备可比性，其灰度差异无法反应地物变化，不能直接进行变化检测；(3)映射分割方式：对其中的一副图像进行分割，并由分割的结果得到每个对象对应像素点的坐标或索引，称为索引矩阵，然后按照此索引矩阵对另一幅图像进行分割，分割的结果使两幅图像具有同样大小的对象。当两幅图像分辨率不一致时，用空间分辨率高的图像来进行分割，当分辨率相同时，按照获取时相的时间顺序，选取后一时期的图像进行分割。本发明采用的是基于改进的FNEA多尺度分割算法，在分割的过程中，不仅考虑颜色，形状和尺度，还加入了MBI指数和纹理值，以此得到更好的建筑物分割效果。

(一)多尺度图像分割采用的是异质性最小的区域合并算法(FNEA算法)，此算法是一种自下而上的区域合并算法，从一个像元起步，先将单像元合并为较小的对象，然后把具有异质性最小的较小对象合并成较大的对象，这样不断合并，直到判断条件不成立，合并操作就终止，最终分割的结果中所有图像对象的平均异质性最小。在判断两相邻区域是否能够合并时，用总异质性值和先前设定好的尺度阈值进行比较，如果小于尺度阈值就合并，否则就结束合并操作。图像的异质性(h)由光谱异质性(h_spectral)和形状异质性(h_shape)决定。

h＝w*h_spectral+(1-w)*h_shape  (0)

w为光谱的权重值，w在0到1范围内。

图像对象的形状异质性由紧密度异质性(h_compact)和光滑度异质性(h_smooth)共同决定的。

h_shape＝w_compact*h_compact+(1-w_compact)*h_smooth  (0)

其中w_compact，w_smooth分别为紧密度和光滑度的权重。在判断两个区域对象是否合并时需要分别计算光谱异质性，紧密度异质性和光滑度异质性，最后求出总的异质性。

FNEA算法的步骤：(1)配置所有与图像分割相关的参数，包括图像各波段的权重值，这个主要依据各波段对于分割过程的影响程度；尺度参数，即用于判断是否需要继续合并操作；根据图像的颜色、色调、纹理等特征来决定光谱异质性和形状异质性的权重值；在形状异质性中根据地物的结构特征来决定紧密度异质性和光滑度异质性的权重值。(2)从图像中任选一个像元，以它为中心启动分割操作，这个过程将把这唯一的像元当作一个最小的多边形对象进行异质性值的计算；第一遍结束后，以最小多边形对象为基础做第二遍分割，也计算其异质性值，比较h与设定阈值，如果h小于阈值，分割操作继续下去，否则就结束分割操作，这样就生成了某一尺度值下的图像对象层。

(二)近年来，MBI指数在提取建筑物特征时能够得到较好的效果，所以本发明把MBI指数加入到多尺度分割算法中；相对于建筑物周围的其他要素，建筑物具有更好的纹理特征，所以本发明也把纹理特征加入到分割算法中，以便能更好的分割出建筑物对象。即本发明的分割算法不仅考虑了尺度参数，形状特征，光谱特征，还考虑了MBI指数和纹理特征。

h＝w_spectral*h_spectral+w_shape*h_shape+w_MBI*h_MBI+w_texture*h_texture  (0)

其中，w_spectral,w_shape,w_MBI,w_texture分别为光谱异质性，形状异质性，MBI异质性和纹理异质性对应的权重。MBI异质性和纹理异质性的计算和光谱异质性的权重相似，即计算每一波段的标准差与该波段权重的乘积，再把个波段的值进行累加。通过对不同的特征异质性赋予不同的权重，与原图像进行比较，通过目视观察，能够较好的把建筑物分割出来。在判断两个区域对象是否合并时也要分别计算光谱异质性，形状异质性(紧密度异质性和光滑度异质性)，MBI异质性和纹理异质性，最后求出总的异质性h，再与阈值比较。

步骤4：因步骤2中得到的MBI特征图像、纹理特征图像和SFA特征图像具有尺度不统一的特点，故采用单位标准差归一化方法对不同的特征图像进行特征优化；

特征归一化的目的就是均衡各特征分量的取值范围，使它们在距离计算中贡献程度大致相同。

单位标准差线性拉伸把分量x拉伸至具有零均值和单位标准差的区间内，其形式为：

$\tilde{x}=\frac{x-u}{\mathrm{\&delta;}}---\left(0\right)$

其中，u为分量样本的均值，δ为分量样本的标准差。

假设特征分量为正态分布，则经归一化后的在区间[-1,1]之间的概率为68％。将分布形式进行平移缩放，有如下形式：

$\tilde{x}=\frac{(x-u)/3\mathrm{\&delta;}+1}{2}---\left(0\right)$

则落入区间[0,1]的概率为99％，超出的部分可以直接赋为0或1。

步骤5：对步骤4归一化后的结果，求每个对象的特征均值，以得到各个对象的特征图像；

步骤6：对每个对象求其在不同时相的特征图像的差值，利用k均值聚类算法对三幅不同的差值图像求其阈值，使其阈值自动化，以此得到三幅不同特征的变化图像；

为了得到变化的结果，需要对不同的图像设定阈值，为了使每一副图像得到的阈值都具有自适应性，本发明提出了利用k均值聚类的方法求得不同特征图像的阈值。K均值聚类算法的工作原理是：

(1)首先，从数据集中随机选取K个图像单元作为初始聚类中心；

(2)然后，计算各个图像单元到聚类中心的光谱距离，将它们一一归类到最近的那个聚类中心所在的类；

(3)计算新形成的每一个聚类的图像单元的光谱均值，从而得到新的聚类中心；

迭代实施步骤(2)、步骤(3)，直至前后两次的聚类中心没有任何变化，说明聚类调整结束，聚类准则函数已经收敛。当输入图像分别为归一化后的MBI特征图像，SFA特征图像和纹理特征图像时，用K均值二值聚类所得的结果，即将这些特征分为两个分离度最高的图像类别，即低相似性与高相似性类别。因此，由此得到的两个图像类别可分别对应于变化区域与未变化区域，即可得到三个阈值X，Y，Z；

步骤7：对步骤6中得到的三幅不同特征的变化图像进行AC指数后处理；

由于道路与建筑物具有很强的相似性，会对检测结果有很大的影响。所以本发明原创性地提出了区别道路与建筑的AC指数来去除道路对建筑物的影响。具体的公式如下：

AC＝a*(area)/circle  (0)

其中，circle是对每一个检测为变化的对象进行圆形拟合，该圆形包含了该对象的所有被检测为变化建筑物的像素以及部分未变化的像素，area为圆形中该对象所有被检测为变化建筑物的像素，a是用来调整该比值的大小。由上述可知，一个道路对象拥有很大的circle值，但其area值相对较小，而一个建筑物对象拥有相对较小的circle值，但其area值相对较大。所有，如果是道路，则其AC值较小，反之为建筑物。通过设定一定的阈值可以把道路和建筑物有效的分开，从而降低错检率；

步骤8：利用熵值法，对步骤7中得到的三幅不同特征的变化图像加不同的权重，再设定阈值，以得到变化检测结果；

权重是评价因子相对于某个评价标准而言的重要性的体现，所以求权重的过程就是分析不同因子对灾害的重要性程度。权重对评价结果起着关键作用，因此，评价指标及权重的合理性直接影响评价结果的科学性和准确性。权重的确定主要有两种方法：一种是主观赋权法，该方法主要是由决策者根据自己主观判断来赋权，会产生一定的主观随意性，准确性不高，如AHP方法、特尔斐法；另一种是客观赋权法，该方法是由各指标数据在评价中经过整理、计算出权重系数，如熵值法、变异系数法。由于主观赋权法主观性太强，因此本发明用客观赋权法-熵值法对不同的特征变化图像进行赋权。

在信息论中，熵是对不确定性的一种度量。信息量越大，不确定性就越小，熵也就越小；信息量越小，不确定性越大，熵也越大。根据熵的特性，可以用熵值来判断某个指标的离散程度，指标的离散程度越大，该指标对综合评价的影响越大。熵值法是一种客观赋权法，在具体使用过程中，其是根据各指标的变异程度，利用信息熵计算出各指标的熵值，再通过熵值对各指标的权重进行修正，从而得出较为客观的指标权重。具体步骤如下：

(1)原始数据的收集与处理

设有m个待评项目，n个评价指标，形成原始数据矩阵x＝{x_ij}_mn：

$X=\left[\begin{array}{cccc}{X}_{11}& X_{12}& ...& X_{1n}\\ X_{21}& X_{22}& ...& X_{2n}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ X_{m1}& X_{m2}& ...& X_{\mathrm{mn}}\end{array}\right]---\left(0\right)$

其中x_ij表示第j个指标下第i个项目的评价指标的数值，在本发明求权重值的方法中，m为3，为各个特征变化图像值，x_ij为第i个特征图像对应的第j个像素的特征值；

(2)计算第j个指标下第i个项目的指标的比重P_ij：

$p_{\mathrm{ij}}=\frac{X_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{X}_{\mathrm{ij}}}(i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n)---\left(0\right)$

(3)计算第j个指标的熵值e_j：

$e_j=-k\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\&CenterDot;\ln \left(p_{\mathrm{ij}}\right),$ 其中， $k=\frac{1}{\ln \left(m\right)},e_j\&GreaterEqual;0---\left(0\right)$

(4)计算第j个指标的熵权w_j：

$w_j=\frac{(1-e_j)}{m-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{e}_j},(0\&le;w_j\&le;1,\underset{1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_j=1)---\left(0\right)$

即由上式步骤得到的w_j即为各个特征变化图像的权重值；

步骤9：利用基于投票法的多尺度融合方法对变化检测的结果进行处理，以得到更高的检测精度；

由于面向对象的变化检测结果很大程度上归结于分割算法的好坏，所以本发明采用基于投票法的多尺度融合方法来减少分割算法对变化检测精度的影响。

在不同的分割尺度下，利用上述方法流程，能获得不同对象尺度的二值变化检测结果图。为了适应图像区域地物类型的多样性，使变化检测结果尽可能地符合不同地物的尺度大小，需要对各尺度下的变化检测结果进行多尺度融合处理。

在本方法中，采用投票法完成这一处理。当某像元在不同尺度下的变化结果中，被检测为变化的次数多于用户定义的投票阈值时，则判定该像元在融合结果中为变化，否则，则视其为未变化，具体公式表示如下：

$M_i=\left\{\begin{array}{cc}1& \begin{array}{cc}\mathrm{if}& n_{\mathrm{ic}}>n_v\end{array}\\ 0& \mathrm{otherwise}\end{array}\right.---\left(0\right)$

其中，M_i表示像元i的变化检测融合结果，n_ic为该像元在各尺度变化检测结果中被判定为变化次数，n_v表示用户定义的变化次数投票阈值；

步骤10：对步骤9得到的结果进行精度评定。

变化检测结果精度评定：变化检测结果的精度评定对成果的应用、算法的有效性验证都具有重要的意义。面向对象变化检测结果的精度受到几何配准、分割、分类、矢量变化检测等诸多因素的影响，每一个步骤都可以单独的进行精度评定。但作为一个整体过程，最为重要的还是最终检测结果的精度评定。

进行变化检测精度评定首先必须获取标准的变化区域，标准变化区域的获取需要通过人工选取来获得。在数据量很大的情况下，人工选出所有的变化区域显然是不现实的，此时可以通过间隔采样的方法，每隔一定的距离选择一个采样区进行标准变化的采集，再将这些采样区所对应的自动变化检测的结果与人工采集的结果进行对比，以此评定变化检测的精度。

变化检测的评定方法类似于分类的评定方法，分类中的类别对应着变化检测中的是否变化。在评定变化检测精度时，可以把变化和未变化这两种变化情况看作是两个类别，类比地构造出如表1所示的矩阵：

表1：精度评价

变化情况	检测变化	检测未变化	合计
				实际变化	s11	s12	S1j
实际未变化	s21	/	/
				合计	Si1	/	/

其中：s₁₁表示实际发生变化并且正确检出的区域的面积，s₁₂表示实际发生变化但未检测出的区域的面积，s₂₁表示实际未发生变化但检测为变化区域的面积。

变化检测精度指标则不同于分类，包括检全率、检准率、漏检率和错检率四个指标。分别定义如下：

(1)检全率

$P_A=\frac{s_{11}}{S_{1j}}\&times;100\%=\frac{s_{11}}{s_{11}+s_{12}}\&times;100\%---\left(0\right)$

(2)检准率

$P_C=\frac{s_{11}}{S_{1j}}*100\%=\frac{s_{11}}{s_{11}+s_{21}}\&times;100\%---\left(0\right)$

(3)漏检率

P_L＝1-P_A  (0)

(4)错检率

P_E＝1-P_C  (0)

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (8)

1.一种基于多特征融合的面向对象的建筑物变化检测方法，其特征在于，主要包括以下步骤：

步骤1：对所选区域的两幅不同时相的高分辨率遥感图像A和B进行预处理，主要包括几何纠正、辐射纠正、几何配准和图像裁剪；

步骤2：计算A和B两幅图像的每个像素点的形态学建筑指数(MBI)、固定窗口的纹理特征和慢特征分析(SFA)，得到图像A和B的MBI特征图像、纹理特征图像和SFA特征图像；

步骤3：选取步骤1中的一幅图像采用基于改进的FNEA多尺度分割方法进行分割，得到一幅多尺度分割图像，并利用步骤2的计算结果来改进多尺度分割算法进行分割，分割的结果是得到一个个对象；并由分割的结果得到每个对象对应像素点的坐标或索引，称为索引矩阵；然后按照此索引矩阵对步骤1中的另一幅图像进行分割，分割的结果使两幅图像具有同样大小的对象；

步骤4：因步骤2中得到的MBI特征图像、纹理特征图像和SFA特征图像具有尺度不统一的特点，故采用单位标准差归一化方法对不同的特征图像进行特征优化；

步骤5：对步骤4归一化后的结果，求每个对象的特征均值，以得到各个对象的特征图像；

步骤6：对每个对象求其在不同时相的特征图像的差值，利用k均值聚类算法对三幅不同的差值图像求其阈值，使其阈值自动化，以此得到三幅不同特征的变化图像；

步骤7：对步骤6中得到的三幅不同特征的变化图像进行AC指数后处理；

步骤8：利用熵值法，对步骤7中得到的三幅不同特征的变化图像加不同的权重，再设定阈值，以得到变化检测结果；

步骤9：利用基于投票法的多尺度融合方法对变化检测的结果进行处理，以得到更高的检测精度；

步骤10：对步骤9得到的结果进行精度评定。

2.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤1中所述的几何纠正采用基于多项式的遥感图像几何纠正，控制点选取分布均匀，重采样采用双线性内插法，最后得到误差要求标准为RMSE<0.5像素；所述的辐射纠正方法采用的是相对辐射归一化纠正。

3.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤3中所述的选取步骤1中的一幅图像进行分割，选取方法为：当两幅图像分辨率不一致时，用空间分辨率高的图像来进行分割；当分辨率相同时，按照获取时相的时间顺序，选取后一时期的图像进行分割。

4.根据权利要求1或3所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤3中所述的基于改进的FNEA多尺度分割方法，其具体实现过程为：从一个像元起步，先将单像元合并为较小的对象，然后把具有异质性最小的较小对象合并成较大的对象，这样不断合并，直到判断条件不成立，合并操作就终止，最终分割的结果中所有图像对象的平均异质性最小；

在判断两相邻区域是否能够合并时，用总异质性值和先前设定好的尺度阈值进行比较，如果小于尺度阈值就合并，否则就结束合并操作；

总的异质性h计算公式为：

h＝w_spectral*h_spectral+w_shape*h_shape+w_MBI*h_MBI+w_texture*h_texture；

其中，w_spectral,w_shape,w_MBI,w_texture分别为光谱异质性h_spectral、形状异质性h_shape、MBI异质性h_MBI和纹理异质性h_texture对应的权重；MBI异质性和纹理异质性的计算方法和光谱异质性的计算方法原理相同，即计算每一波段的标准差与该波段权重的乘积，再把个波段的值进行累加。

5.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤6中所述的利用k均值聚类算法对三幅不同的差值图像求其阈值，其具体实现包括以下子步骤：

步骤6.1：从数据集中随机选取K个图像单元作为初始聚类中心；

步骤6.2：计算各个图像单元到聚类中心的光谱距离，将它们一一归类到最近的那个聚类中心所在的类；

步骤6.3：计算新形成的每一个聚类的图像单元的光谱均值，从而得到新的聚类中心；

步骤6.4：迭代实施步骤6.2和步骤6.3，直至前后两次的聚类中心没有任何变化，说明聚类调整结束，聚类准则函数已经收敛；

当输入图像分别为归一化后的MBI特征图像、SFA特征图像和纹理特征图像时，用K均值二值聚类所得的结果，即将这些特征分为两个分离度最高的图像类别，即低相似性与高相似性类别；因此，由此得到的两个图像类别可分别对应于变化区域与未变化区域，即可得到三个阈值X，Y，Z。

6.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤7中所述的AC指数后处理，具体的公式如下：

AC＝a*(area)/circle

；

其中，circle是对每一个检测为变化的对象进行圆形拟合，该圆形包含了该对象的所有被检测为变化建筑物的像素以及部分未变化的像素，area为圆形中该对象所有被检测为变化建筑物的像素，a是用来调整该比值的大小。

7.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤8的具体实现包括以下子步骤：

步骤8.1：原始数据的收集与处理

设有m个待评项目，n个评价指标，形成原始数据矩阵x＝{x_ij}_mn：

$x=\left[\begin{array}{cccc}{x}_{11}& x_{12}& ...& x_{1n}\\ x_{21}& x_{22}& ...& x_{2n}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_{m1}& x_{m2}& ...& x_{\mathrm{mn}}\end{array}\right];$

其中x_ij表示第j个指标下第i个项目的评价指标的数值，m为各个特征变化图像值，m＝3，x_ij为第i个特征图像对应的第j个像素的特征值；

步骤8.2：计算第j个指标下第i个项目的指标的比重P_ij：

$p_{i,j}=\frac{x_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ij}}}(i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n);$

步骤8.3：计算第j个指标的熵值e_j：

$e_j=-k\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\&CenterDot;\ln \left(p_{\mathrm{ij}}\right),$ 其中， $k=\frac{1}{\ln \left(m\right)},e_j\&GreaterEqual;0;$

步骤8.4：计算第j个指标的熵权w_j：

$w_j=\frac{(1-e_j)}{m-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{e}_j},(0\&le;w_j\&le;1,\underset{1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_j=1);$

即由上式步骤得到的w_j即为各个特征变化图像的权重值。

8.根据权利要求1所述的基于多特征融合的面向对象的建筑物变化检测方法，其特征在于：步骤9所述的利用基于投票法的多尺度融合方法对变化检测的结果进行处理，其具体实现过程为：当某像元在不同尺度下的变化结果中，被检测为变化的次数多于用户定义的投票阈值时，则判定该像元在融合结果中为变化，否则，则视其为未变化，具体公式表示如下：

$M_i=\left\{\begin{array}{cc}1& \mathrm{if}{n}_{\mathrm{ic}}{>n}_v\\ 0& \mathrm{otherwise}\end{array}\right.$

其中，M_i表示像元i的变化检测融合结果，n_ic为该像元在各尺度变化检测结果中被判定为变化次数，n_v表示用户定义的变化次数投票阈值。