CN105069811A - 一种多时相遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种多时相遥感图像变化检测方法，它包括以下步骤：步骤1、计算差分图像，步骤2、检测SURF特征点，步骤3、SURF特征点匹配，步骤4、提取非变化类像素样本，步骤5、提取变化类像素样本，步骤6、构造分类模型，步骤7、对差分图像的所有像素分类实现变化检测。本发明具有以下优点：(1)？本发明通过将SURF特征点用于变化检测，能有效地克服噪声对变化检测精度造成的不利影响，有利于对背景较为复杂的遥感图像进行处理，提高了变化检测的精确度；(2)？本发明不需要有关遥感数据分布的任何先验知识，能适用于更多类型的遥感数据；(3)？本发明采用支持向量机分类模型，充分利用了支持向量机具有的良好的泛化性能，使本发明具有良好的鲁棒性。

Description

一种多时相遥感图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，具体涉及一种基于SURF特征点和支持向量机的多时相遥感图像变化检测方法，适用于土地利用、城市规划及军事目标监测等领域中的数字图像处理。

背景技术

遥感图像变化检测是通过分析同一目标或区域的多时相遥感图像，从中获取变化的信息。随着空间遥感技术的发展，遥感变化检测已被广泛应用于土地利用、灾害评估与预测、地理数据库更新、气候监测、军事目标监测等领域，对于推动环境保护、经济发展、社会发展和国防建设起到了重要的作用。

遥感图像变化检测方法，总体上分为监督变化检测及非监督变化检测等两类。非监督变化检测是直接对两个不同时相的同一区域的遥感数据进行比较，不需要有关研究区域的任何先验知识，这一点在遥感变化检测实际应用中显得尤为重要，对非监督变化检测方法的研究，已经成为当前变化检测方法研究的热点。非监督变化检测首先对两时相遥感图像进行预处理(空间配准、辐射校正等)；然后通过求差、求比，或者主成分分析、变化矢量分析变换等方法获得差分图像；最后对差分图像进行分割来确定变化和非变化区域，从而得到最终的变化检测结果。

学者们提出了很多种遥感非监督变化检测方法，常见的一类是基于聚类的遥感图像变化检测方法。“Multiscalechangedetectioninmultitemporalsatelliteimages”，T.Celik，IEEEGeoscienceandRemoteSensingLetters，6(4)，820～824，2009，proposedanunsupervisedchangedetectiontechniquebyconductingk-means(k＝2)clusteringonfeaturevectorswhichareextractedusingthesubbandsoftheUDWTdecompositionofthedifferenceimage(“基于多尺度分析的多时相卫星图像变化检测”，图尔盖.赛因切利克，地球科学与遥感快报，第六卷第4期，第820～824页，2009年提出一种基于多尺度分析的多时相卫星图像非监督变化检测方法，首先构造差分图像并对差分图像实施非抽样离散小波分解从而提取特征向量，然后再利用k(k＝2)均值聚类的方法对特征向量进行分类，最后得到变化及非变化区域)。“Fuzzyclusteringalgorithmsforunsupervisedchangedetectioninremotesensingimages”，A.GhoshandN.S.MishraandS.Ghosh，InformationSciences，181(4)，699～715，2011，proposedanunsupervisedcontext-sensitivetechniquesusingfuzzyclusteringapproachfordetectingchangesinmultitemporal，multispectralremotesensingimages.Thetechniqueisbasedonfuzzyclusteringapproachandtakescareofspatialcorrelationbetweenneighboringpixelsofthedifferenceimageproducedbycomparingtwoimagesacquiredonthesamegeographicalareaatdifferenttimes(“基于模糊聚类的遥感图像非监督变化检测”，阿希什.戈什、尼尔纳德瑞.谢卡尔.米什拉、莎士米塔.戈什，信息科学，第一百八十一卷第4期，第699～715页，2011年提出采用模糊聚类的方法实现一种上下文相关的多时相多光谱遥感图像非监督变化检测技术。该技术首先对同一地域不同时相的遥感图像进行比较以构造差分图像，然后基于差分图像中像素的光谱值并结合相邻像素间的空间相关性、采用模糊聚类的方法对差分图像进行分类，最后得到变化及非变化区域)。

遥感图像成像时由各种原因所产生的噪声对变化检测的性能有着很大的影响。上述利用聚类的方法进行变化检测时，当噪声较强时，其检测性能会受到很大的影响，导致变化检测的准确性差。

发明内容

为克服上述基于聚类的遥感图像变化检测方法存在的问题，本发明提出一种多时相遥感图像变化检测方法，它能有效地、准确地检测出变化和非变化区域，降低噪声对变化检测精度造成的不利影响，提高遥感图像变化检测的精确度。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括有以下步骤：

1、输入两幅已经过预处理的多时相遥感图像I₁及I₂，并通过计算两幅图像间对应位置像素光谱值差值的绝对值构造差分图像I_d，所述预处理包括空间配准和辐射校正；

2、分别从图像I₁及I₂中检测SURF特征点(加速稳健特征点(SpeededUpRobustFeatures，SURF)，以下简称为特征点)，得到特征点集合S₁及S₂；

3、对集合S₁与集合S₂中的特征点进行匹配，得到集合S₁、S₂间相互匹配的特征点集合M₁、M₂，以及不能相互匹配的特征点集合R₁、R₂；集合M₁、M₂中的特征点一一对应、构成相互匹配的特征点对，其中

4、提取非变化类像素样本

对集合M₁、M₂中的每个特征点，计算其在差分图像I_d中对应位置的像素的亮度值并归入非变化类像素样本集合S_u；

5、提取变化类像素样本

将集合R₁、R₂合并为集合T，并计算T中每个像素点在差分图像I_d中对应位置的像素的亮度值并归入集合G；假定集合T中的像素由非变化类、未标记类及变化类像素构成，且每类像素的亮度均服从高斯分布，分别记为 则总体亮度直方图应为三维混合高斯分布：

$Z\left(g\right)=w_{1}N\left(g\right|{\mathrm{\μ}}_1,{\mathrm{\δ}}_1^2)+w_{2}N\left(g\right|{\mathrm{\μ}}_2,{\mathrm{\δ}}_2^2)+w_{3}N\left(g\right|{\mathrm{\μ}}_3,{\mathrm{\δ}}_3^2)$

式中，Z(g)表示三维混合高斯分布，g∈G为像素亮度值；w₁、w₂及w₃分别是非变化类、未标记类及变化类像素高斯分布在混合高斯中的权重；表示非变化类像素的亮度服从均值为μ₁、方差为的高斯分布；表示未标记类像素的亮度服从均值为μ₂、方差为的高斯分布；表示变化类像素的亮度服从均值为μ₃、方差为的高斯分布；

利用EM算法对上述模型的相关参数进行求解，最终将集合G中符合下式的像素亮度值归入变化类像素样本集合S_c：

S_c＝{g|g∈G且g≥μ₃-3δ₃}

6、依据步骤4所得的非变化类像素样本集合S_u和步骤5所得的变化类像素样本集合S_c构造训练集S_train、测试集S_test；再构造相应的训练集标签L_train、测试集标签L_test，其中集合L_train中的元素全为0且同集合S_train的大小相等，集合L_test中的元素全为1且同集合S_test的大小相等；最后在对测试集S_train、训练集S_test进行预处理的基础上训练一个支持向量机分类模型M_svm；

7、依据步骤6所得的支持向量机分类模型M_svm对差分图像I_d的所有像素进行分类，得到变化及非变化区域，最终的变化检测结果将以二值变化检测掩膜的形式体现，其中值为0的像素表示未发生变化的类别，值为1的像素表示发生变化的类别。

由于多时相遥感图像变化检测的全部变化信息由若干局部地物的变化信息综合构成，SURF特征点是图像特征的局部表达，能够很好地反映出图像的局部特殊性，且对图像噪声也具有较好的鲁棒性。本发明采用基于SURF特征点和支持向量机的多时相遥感图像变化检测方法，将SURF局部特征点和支持向量机应用于变化检测，突破传统遥感变化检测方法的约束，与现有的遥感变化检测方法相比，本发明具有以下优点：

(1)本发明通过将SURF特征点用于变化检测，能有效地克服噪声对变化检测精度造成的不利影响，有利于对背景较为复杂的遥感图像进行处理，提高了变化检测的精确度。

(2)本发明不需要有关遥感数据分布的任何先验知识，能适用于更多类型的遥感数据。

(3)本发明采用支持向量机分类模型，充分利用了支持向量机具有的良好的泛化性能，使本发明具有良好的鲁棒性。

附图说明

图1是本发明的流程图；

图2是第一时相遥感图像；

图3是第二时相遥感图像；

图4是参考变化图；

图5是用现有的多尺度分析方法对图2和图3进行变化检测的结果；

图6是用现有的模糊聚类的方法对图2和图3进行变化检测的结果；

图7是用本发明的方法对图2和图3进行变化检测的结果。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：

如图1所示，本发明具体实施方式的步骤如下：

步骤1，计算差分图像

输入两幅经过预处理的多时相遥感图像I₁及I₂，定义差分图像I_d中位于第i行第j列的像素的亮度值为两幅图像间对应位置(i,j)处像素光谱值I₁(i,j)和I₂(i,j)差值的绝对值，即：

I_d(i,j)＝|I₁(i,j)-I₂(i,j)|(1)

式(1)中，1≤i≤H,1≤j≤W，H为遥感图像的行数，W为遥感图像的列数。

步骤2，检测SURF特征点

SURF特征点能很好地反映图像的局部特征，对图像旋转、尺度缩放和亮度变化等具有不变性，对视角变化、光照变化及噪声等也具有较好的鲁棒性，被广泛应用于图像配准和目标识别等领域。本发明将SURF特征点应用于变化检测，利用SURF特征点对噪声具有较好的鲁棒性，降低噪声对变化检测精度的影响。从遥感图像I₁及I₂中检测SURF特征点，具体包括以下步骤：

(2.1)计算积分图像

对于任意一幅图像I，位于该图像中第x行第y列的像素点的积分图像值I_Σ(x,y)是指从图像原点到以该像素点为顶点的矩形区域内所有像素的像素值之和。则遥感图像I₁或I₂中位于第x行第y列的像素点的积分图像值为：

$I_{t\mathrm{\Σ}}=\underset{i=0}{\overset{i\≤x}{\Σ}}\underset{j=0}{\overset{j\≤y}{\Σ}}{I}_t(x,y),t\∈\{1,2\}---\left(2\right)$

式(2)中，t为遥感图像I₁或I₂的下角标代号。

(2.2)构建多尺度空间

现实世界的物体均是由不同尺度的结构所组成，于是在计算机视觉领域，引入一个被视为尺度的参数，通过连续变化尺度参数构造多尺度空间而获得不同尺度下的视觉处理信息，然后综合这些信息以便更为深入地挖掘图像的本质特征。为检测出图像中的SURF特征点，也必须要构造多尺度空间。高斯函数是多尺度分析的最优方式，但在实际应用中，高斯函数经常需要被离散化和剪切，且这种离散和剪切的近似过程是不可避免的。

本发明利用盒子滤波器近似表示高斯二阶偏导数来构造多尺度空间。使用的初始盒子滤波器的模板大小为9×9，其对应的尺度大小为σ＝1.2，保持原始图像的大小不变，通过不断增加盒子滤波器的大小，并将其在原始图像的x、y及xy等三个方向上做卷积，最终得到多尺度空间。盒子滤波器模板的尺寸以公差大小为6的方式递增，若盒子滤波器模板的尺寸为N×N，则其对应的尺度大小为：

$\mathrm{\σ}=1.2\×\frac{N}{9}---\left(3\right)$

(2.3)计算近似Hessian矩阵

对于遥感图像I₁和I₂中坐标为(x,y)的像素点X，根据步骤(2.1)中计算出的积分图像值和步骤(2.2)中所构建的多尺度空间，按下式计算出该点在尺度σ上所对应的近似Hessian矩阵H₁(X,σ)和H₂(X,σ)：

$H_t(X,\mathrm{\σ})=\left[\begin{array}{cc}{D}_{xxt}(X,\mathrm{\σ})& {\mathrm{\λD}}_{xyt}(X,\mathrm{\σ})\\ {\mathrm{\λD}}_{xyt}(X,\mathrm{\σ})& D_{yyt}(X,\mathrm{\σ})\end{array}\right],t\∈\{1,2\}---\left(4\right)$

式(4)中，t为遥感图像I₁或I₂的下角标代号；D_xxt(X,σ)是由步骤(2.1)得到的积分图像值与步骤(2.2)中所述的x方向上的盒子滤波器相乘求和的结果；D_xyt(X,σ)是由步骤(2.1)得到的积分图像值与步骤(2.2)中所述的xy方向上的盒子滤波器相乘求和的结果；D_yyt(X,σ)是由步骤(2.1)得到的积分图像值与步骤(2.2)中所述的y方向上的盒子滤波器相乘求和的结果；λ为权重系数，在本发明中λ取值为0.9。矩阵H_t(X,σ)对应行列式的值det(H_t)为：

det(H_t)＝D_xxtD_yyt-(λD_xyt)²(5)

对遥感图像I₁及I₂中的每一个像素，均可以根据式(5)计算出该像素在尺度σ上的响应值。

(2.4)确定极值点

首先计算近似Hessian矩阵H_t(X,σ)的特征值，因为H_t(X,σ)为二阶方阵，所以H_t(X,σ)有两个特征值，分别记为θ_t1、θ_t2，显然：

θ_t1×θ_t2＝det(H_t),t∈{1,2}(6)

式(6)中，t为遥感图像I₁或I₂的下角标代号。

然后根据步骤(2.3)中计算出的近似Hessian矩阵对应行列式的值det(H_t)初步选择极值点。如果det(H_t)小于零，则特征值θ_t1、θ_t2有不同的符号，该像素点不是极值点；反之，如果det(H_t)大于零，则特征值θ_t1、θ_t2同为正或负，该像素点是极值点。再将各极值点所对应的行列式的值与预设的阈值进行比较，所有小于该阈值的极值点都将被剔除，仅保留具有强响应值的极值点。阈值的大小需要根据实际的应用而进行设定，若增加阈值，则检测出的特征点会比较少；如果减小阈值，则检测出的特征点会比较多。

(2.5)建立SURF特征点的描述向量

(2.5.1)确定特征点的主方向

检测出特征点之后，为了保证特征点具有旋转和尺度不变性，需要确定特征点的主方向。以特征点为中心，在其周围构建一个半径为6σ的圆形区域，计算该区域内的点在x轴与y轴方向的Haar小波(Haar小波的边长为4σ)响应，其中σ为该特征点对应的尺度。然后在圆形区域中，用圆心角为π/3的扇形以特征点为中心、0.2弧度为步长旋转一圈，计算在每个角度下该扇形区域范围内的所有像素在x、y两个方向的Haar小波响应的和，并将其构成一个矢量，最后选择最大响应值的矢量所在的方向作为该特征点的主方向。

(2.5.2)建立特征点描述向量

为构建特征点描述向量，首先以特征点为原点并将坐标轴旋转到该特征点主方向上，然后按主方向确定一个以特征点为中心的正方形区域(正方形的边长为20σ，σ是该特征点对应的尺度)，最后将该正方形区域分成4×4个大小相等的子区域。每个子区域的边长为5σ，内有25个均匀分布的采样点，对每个子区域用尺度为2σ的Haar滤波器进行处理，得到x、y两个方向的响应，分别用dx、dy表示。计算每个子区域内的dx的响应值之和Σdx，|dx|的响应值之和Σ|dx|，dy的响应值之和Σdy，以及|dy|的响应值之和Σ|dy|，再将Σdx、Σ|dx|、Σdy及Σ|dy|构成每个子区域的特征描述向量：

V_subregion＝[Σd_xΣd_yΣ|d_x|Σ|d_y|](7)

对于每一个特征点，将4×4个子区域的特征描述向量连接起来，就可以得到一个4×4×4维的向量，这个向量就是该特征点的特征描述向量V_subregion。对这些向量进行归一化处理之后，特征点就具有旋转、尺度和光照不变性。

步骤3，SURF特征点匹配

对图像I₁与I₂中的SURF特征点进行匹配，具体步骤为：

(3.1)判断特征点之间特征相似性的高低

设从图像I₁与I₂中检测出的SURF特征点集合分别为S₁＝{P_1k|1≤k≤N₁}及S₂＝{P_2q|1≤q≤N₂}，其中P_1k、P_2q分别为图像I₁和I₂中的特征点，N₁和N₂分别为从图像I₁和I₂中检测出的特征点个数。对集合S₁中的特征点P_1k，按式(8)计算其与集合S₂中所有特征点的特征描述向量间的欧氏距离，从而判断P_1k与图像I₂中每个特征点之间特征相似性的高低。

$Dis(P_{1k},P_{2q})={\[\underset{h=1}{\overset{64}{\Σ}}{(m_{1kh}-m_{2qh})}^2\]}^{1/2}---\left(8\right)$

式(8)中，m_1kh表示集合S₁中的特征点P_1k的特征描述向量的第h个元素；m_2qh表示集合S₂中的特征点P_2q的特征描述向量的第h个元素；h为特征点的64维特征描述向量的元素序号；Dis(P_1k,P_2q)用于度量特征点P_1k与特征点P_2q间特征相似性的高低，相似性越低，欧氏距离值越大，反之，欧氏距离值越小。

(3.2)判断特征点之间是否满足8-邻接关系

对于任意一幅图像I中的某个像素点P，其坐标为(x,y)，与像素点P在水平方向(左和右)和垂直方向(上和下)相邻的4个像素点的坐标分别为(x,y-1)、(x,y+1)、(x-1,y)及(x+1,y)；与像素点P在对角线方向(45°、135°方向)上相邻的4个像素点的坐标分别为(x+1,y-1)、(x-1,y+1)、(x-1,y-1)及(x+1,y+1)，这8个像素点在位置上与像素点P构成8-邻接关系。

假设在集合S₂中的特征点P_2q与集合S₁中的特征点P_1k最相似，则进一步根据P_1k与P_2q的坐标值，按式(9)来判断二者在差分图像I_d中是否满足8-邻接关系，从而判断二者是否相互匹配。

$Match(P_{1k},P_{2q})=\left\{\begin{array}{cc}true,& \sqrt{{(x_{1k}-x_{2q})}^2+{(y_{1k}-y_{2q})}^2}\≤\sqrt{2}\\ false,& \sqrt{{(x_{1k}-x_{2q})}^2+(y_{1k}-y_{2q})}>\sqrt{2}\end{array}\right.---\left(9\right)$

式(9)中，(x_1k,y_1k)、(x_2q,y_2q)分别为特征点P_1k、P_2q的位置坐标；Match(P_1k,P_2q)用于度量特征点P_1k与特征点P_2q间是否满足8-邻接关系，若值为true，则P_1k与P_2q满足8-邻接关系(即相互匹配)，否则表示P_1k与P_2q不满足8-邻接关系(即不能相互匹配)。

(3.3)将匹配的特征点与不匹配的特征点进行归类

对集合S₁中的每个特征点，重复执行步骤(3.1)至步骤(3.2)，最终得到集合S₁与集合S₂间所有能相互匹配的特征点。在这些相互匹配的特征点中，凡是属于集合S₁的，则将该特征点归入集合M₁，凡是属于集合S₂的，则将该特征点归入集合M₂。集合M₁、M₂中的特征点一一对应、构成相互匹配的特征点对。对于集合S₁中的部分特征点而言，若在集合S₂中不存在与之匹配的特征点，则将这些特征点归入集合R₁；对于集合S₂中的部分特征点而言，若在集合S₁中不存在与之匹配的特征点，则将这些特征点归入集合R₂。集合M₁、M₂、R₁、R₂、S₁及S₂满足如下关系：

步骤4，提取非变化类像素样本

(4.1)对集合M₁、M₂中坐标为(x,y)的特征点，按式(12)来计算其在差分图像I_d中对应位置的像素的亮度值g(x,y)：

当x或y不为整数时，采用双线性插值的方式来计算该特征点在差分图像I_d中对应位置的像素的亮度值，其中，x_u为x的小数部分，x_i为x的整数部分，y_u为y的小数部分，y_i为y的整数部分。

(4.2)对集合M₁、M₂中的每个特征点，重复执行步骤(4.1)，计算其在差分图像I_d中对应位置的像素的亮度值并归入非变化类像素样本集合S_u。

步骤5，提取变化类像素样本

(5.1)将集合R₁、R₂合并为集合T，即T＝R₁∪R₂。

(5.2)对集合T中的每个特征点，按式(12)来计算其在差分图像I_d中对应位置的像素的亮度值并归入集合G。

(5.3)假定集合T中的像素点由非变化类、未标记类及变化类构成，且每类像素点的亮度均服从高斯分布，则总体亮度直方图应为三维混合高斯分布，如式(13)所示：

$Z\left(g\right)=w_{1}N\left(g\right|{\mathrm{\μ}}_1,{\mathrm{\δ}}_1^2)+w_{2}N\left(g\right|{\mathrm{\μ}}_2,{\mathrm{\δ}}_2^2)+w_{3}N\left(g\right|{\mathrm{\μ}}_3,{\mathrm{\δ}}_3^2)---\left(13\right)$

式(13)中，Z(g)表示三维混合高斯分布，g∈G为像素亮度值；w₁、w₂及w₃分别是非变化类、未标记类及变化类像素亮度的高斯分布在混合高斯中的权重；表示非变化类像素的亮度服从均值为μ₁、方差为的高斯分布；表示未标记类像素的亮度服从均值为μ₂、方差为的高斯分布；表示变化类像素的亮度服从均值为μ₃、方差为的高斯分布。显然，上述参数满足如下约束条件：

0＜w₁,w₂,w₃＜1(14)

w₁+w₂+w₃＝1(15)

μ₁＜μ₂＜μ₃(16)

(5.3.1)初始化高斯混合模型

使用EM算法对各高斯成分的参数进行估计求解。首先按式(17)-(22)对式(13)所示高斯混合模型的相关参数进行初始化：

$w_1^0=\frac{card\left(F_1\right)}{H\×W},w_2^0=\frac{card\left(F_2\right)}{H\×W},w_3^0=\frac{card\left(F_3\right)}{H\×W}---\left(20\right)$

${\mathrm{\μ}}_1^0=\frac{\underset{g\∈F_1}{\Σ}g}{card\left(F_1\right)},{\mathrm{\μ}}_2^0=\frac{\underset{g\∈F_2}{\Σ}g}{card\left(F_2\right)},{\mathrm{\μ}}_3^0=\frac{\underset{g\∈F_3}{\Σ}g}{card\left(F_3\right)}---\left(21\right)$

${\left({\mathrm{\δ}}_1^2\right)}^0=\frac{\underset{g\∈F_1}{\Σ}{\[g-{\mathrm{\μ}}_1^0\]}^2}{card\left(F_1\right)},{\left({\mathrm{\δ}}_2^2\right)}^0=\frac{\underset{g\∈F_2}{\Σ}{\[g-{\mathrm{\μ}}_2^0\]}^2}{card\left(F_2\right)},{\left({\mathrm{\δ}}_3^2\right)}^0=\frac{\underset{g\∈F_3}{\Σ}{\[g-{\mathrm{\μ}}_3^0\]}^2}{card\left(F_3\right)}---\left(22\right)$

式17—22中，g为像素亮度值；β∈(0,1)，其大小可视情况自行设定；F₁、F₂、F₃分别表示集合T中初始的非变化类像素、未标记类像素及变化类像素集合；card(·)表示集合中元素的个数；分别为非变化类像素亮度高斯分布的初始权重、初始均值及初始方差；分别为未标记类像素亮度高斯分布的初始权重、初始均值及初始方差；分别为变化类像素亮度高斯分布的初始权重、初始均值及初始方差。

(5.3.2)高斯混合模型参数估计

设n＝0,1,2,...为迭代次数，Ω₁、Ω₂、Ω₃分别表示非变化类像素、未标记类像素及变化类像素类型。

计算非变化类像素亮度高斯分布的参数及比重的迭代步骤为：

$w_1^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_1^n{P}^n\left(g\right|{\mathrm{\Ω}}_1)}{P^n\left(g\right)}}{HW},{\mathrm{\μ}}_1^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_1^n{P}^n\left(g\right|{\mathrm{\Ω}}_1)}{P^n\left(g\right)}g}{\underset{g\∈G}{\Σ}\frac{w_1^n{P}^n\left(g\right|{\mathrm{\Ω}}_1)}{P^n\left(g\right)}},{\left({\mathrm{\δ}}_1^2\right)}^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_1^n{P}^n\left(g\right|{\mathrm{\Ω}}_1)}{P^n\left(g\right)}{(g-{\mathrm{\μ}}_1^n)}^2}{\underset{g\∈G}{\Σ}\frac{w_1^n{P}^n\left(g\right|{\mathrm{\Ω}}_1)}{P^n\left(g\right)}}---\left(23\right)$

计算未标记类像素亮度高斯分布的参数及比重的迭代步骤为：

$w_2^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_2^n{P}^n\left(g\right|{\mathrm{\Ω}}_2)}{P^n\left(g\right)}}{HW},{\mathrm{\μ}}_2^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_2^n{P}^n\left(g\right|{\mathrm{\Ω}}_2)}{P^n\left(g\right)}g}{\underset{g\∈G}{\Σ}\frac{w_2^n{P}^n\left(g\right|{\mathrm{\Ω}}_2)}{P^n\left(g\right)}},{\left({\mathrm{\δ}}_2^2\right)}^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_2^n{P}^n\left(g\right|{\mathrm{\Ω}}_2)}{P^n\left(g\right)}{(g-{\mathrm{\μ}}_2^n)}^2}{\underset{g\∈G}{\Σ}\frac{w_2^n{P}^n\left(g\right|{\mathrm{\Ω}}_2)}{P^n\left(g\right)}}---\left(24\right)$

计算变化类像素亮度高斯分布的参数及比重的迭代步骤为：

$w_3^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_3^n{P}^n\left(g\right|{\mathrm{\Ω}}_3)}{P^n\left(g\right)}}{HW},{\mathrm{\μ}}_3^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_3^n{P}^n\left(g\right|{\mathrm{\Ω}}_3)}{P^n\left(g\right)}g}{\underset{g\∈G}{\Σ}\frac{w_3^n{P}^n\left(g\right|{\mathrm{\Ω}}_3)}{P^n\left(g\right)}},{\left({\mathrm{\δ}}_3^2\right)}^{n+1}=\frac{\underset{g\∈G}{\Σ}\frac{w_3^n{P}^n\left(g\right|{\mathrm{\Ω}}_3)}{P^n\left(g\right)}{(g-{\mathrm{\μ}}_3^n)}^2}{\underset{g\∈G}{\Σ}\frac{w_3^n{P}^n\left(g\right|{\mathrm{\Ω}}_3)}{P^n\left(g\right)}}---\left(25\right)$

式23—25中，g∈G为像素亮度值；分别表示在第n次迭代计算时非变化类像素亮度高斯分布的权重、均值及方差；分别表示在第n次迭代计算时未标记类像素亮度高斯分布的权重、均值及方差；分别表示在第n次迭代计算时变化类像素亮度高斯分布的权重、均值及方差；Pⁿ(g|Ω₁)、Pⁿ(g|Ω₂)、Pⁿ(g|Ω₃)是第n次迭代计算时的条件概率，Pⁿ(g)是第n次迭代计算时的全概率，且满足：

$P^n\left(g\right|{\mathrm{\Ω}}_1)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_1^n}\exp \[-\frac{{(g-{\mathrm{\μ}}_1^n)}^2}{2{\left({\mathrm{\δ}}_1^2\right)}^n}\]---\left(26\right)$

$P^n\left(g\right|{\mathrm{\Ω}}_2)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_2^n}\exp \[-\frac{{(g-{\mathrm{\μ}}_2^n)}^2}{2{\left({\mathrm{\δ}}_2^2\right)}^n}\]---\left(27\right)$

$P^n\left(g\right|{\mathrm{\Ω}}_3)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_3^n}\exp \[-\frac{{(g-{\mathrm{\μ}}_3^n)}^2}{2{\left({\mathrm{\δ}}_3^2\right)}^n}\]---\left(28\right)$

$P^n\left(g\right)=P^n\left(g\right|{\mathrm{\Ω}}_1)w_1^n+P^n\left(g\right|{\mathrm{\Ω}}_2)w_2^n+P^n\left(g\right|{\mathrm{\Ω}}_3)w_3^n---\left(29\right)$

(5.4)将集合G中符合式(30)的像素亮度值归入变化类像素样本集合S_c：

S_c＝{g|g∈G且g≥μ₃-3δ₃}(30)

μ₃、δ₃分别为变化类像素亮度的高斯分布的均值及标准差。

步骤6，构造分类模型

(6.1)选定训练集和测试集

将非变化类像素样本集合S_u近似等分为两组，记为S_u1、S_u2；同样地再将变化类像素样本集合S_c近似等分成两组，记为S_c1、S_c2，且：

再按下述方式构造训练集S_train、测试集S_test：

S_train＝S_u1∪S_c1(33)

S_test＝S_u2∪S_c2(34)

最后将训练集S_train中每个元素的标签值设为0，测试集S_test中每个元素的标签值设为1，得到训练集标签L_train以及测试集标签L_test，集合S_train和L_train的大小相等，集合S_test和L_test的大小相等。

(6.2)数据预处理

对训练集S_train和测试集S_test进行归一化预处理，采用的归一化方式如下：

$S_r=\frac{S_r-min\left(S_r\right)}{max\left(S_r\right)-min\left(S_r\right)},r\∈\{train,test\}---\left(35\right)$

其中，r为训练集S_train和测试集S_test的下角标代号。

(6.3)训练支持向量机分类模型

本发明使用台湾大学林智仁教授等开发设计的Matlab版libsvm-mat-2.89-3软件包作为分类工具，该工具包可在http://www.csie.ntu.edu.tw/～cjlin/libsvm/index.html#download下载。

Libsvm工具包提供了用于训练支持向量机分类模型的接口：

model＝svmtrain(train_label,train_set,['libsvm_options'])(36)

其中，train_label为训练集的标签(即L_train)、train_set为训练集(即S_train)、libsvm_options为一些控制参数、model为训练后得到的支持向量机分类模型。

Libsvm工具包还提供了用于对测试集进行预测的接口：

[predicted_label,accuracy]＝svmpredict(test_label,test_set,model)(37)

其中，test_label为测试集的标签(即L_test)、test_set为测试集(即S_test)、model为支持向量机分类模型、predicted_label为预测后得到的测试集的标签、accuracy为分类准确率。

Libsvm工具包提供了很多默认的参数，本发明仅对惩罚参数Ch和核函数参数Ke进行调节，其它参数均使用默认设置。本发明选用径向基函数作为核函数，并采用交叉验证的方法选择最优的惩罚参数及核函数参数：

①初始化惩罚参数Ch及核函数参数Ke的取值范围，令Ch∈[0,Ch_max]，Ke∈[0,Ke_max]；

②分别从区间[0,Ch_max]及[0,Ke_max]中以步长为1提取出M_Ch个惩罚参数及N_ke个核函数参数。显然，这M_Ch个惩罚参数及N_ke个核函数参数一共可构成M_Ch×N_ke种参数组合，且每个参数组合中仅包含一个惩罚参数及一个核函数参数；

③将每一种参数组合依次作为支持向量机分类模型的惩罚参数及核函数参数，再用训练集S_train对分类器进行训练获得分类模型，最后用分类模型对测试集S_test进行预测。最终选择使得测试集S_test分类准确率最高的参数组合作为支持向量机分类模型的参数，将该最优参数组合记为(Ch_optimal,Ke_optimal)，最优支持向量机分类模型记为M_svm。如果有多个参数组合对应于最高分类准确率，则选取能够达到最高分类准确率中参数Ch最小的那组参数组合作为最优参数；如果对应于最小的Ch有多组Ke，则选取搜索到的第一个参数组合作为最优参数。

步骤7，对差分图像I_d的所有像素进行分类实现变化检测

(7.1)将差分图像I_d中的每个像素对应的标签值设为0，构造标签集合L_predict；

(7.2)利用最优分类模型M_svm对差分图像I_d中的所有像素是否发生变化进行预测，获得变化检测结果：

[CM,accuracy]＝svmpredict(L_predict,I_d,M_svm)(38)

CM＝{cm(i,j)|1≤i≤H,1≤j≤W},cm(i,j)∈{0,1}(39)

其中，CM为预测结果，即二值变化检测结果。CM中值为0的像素表示未发生变化的类别，值为1的像素表示发生变化的类别。

本发明的效果可以通过以下实验结果与分析进一步说明：

1.实验数据

本发明所使用的实验数据为一组真实的QuickBird多时相遥感图片，参见图2和图3，并附有参考变化图，参见图4。该组遥感图像的大小为238×238像素，为北京市石景山区某区域全色多光谱融合数据，共包含4个波段，空间分辨率为2.4米。两时相遥感图像分别拍摄于2008年10月11日和2009年9月13日。从图2和图3可以看出，主要的变化部分是由于人类活动而导致地面新增了一些面积较大的人工建筑。

2.实验内容

利用上述实验数据，分别采用多尺度分析变化检测方法、模糊聚类变化检测方法和本发明的变化检测方法对图2和图3所示的实验数据进行变化检测实验，根据得到的变化检测结果对各变化检测方法的性能进行对比分析。

3.实验结果与分析

用不同的变化检测方法对实验数据进行变化检测，检测结果为：

(1)采用多尺度分析变化检测方法对图2和图3所示的实验数据进行变化检测，其结果如图5所示；

(2)采用模糊聚类变化检测方法对图2和图3所示的实验数据进行变化检测，其结果如图6所示；

(3)采用本发明的变化检测方法对图2和图3所示的实验数据进行变化检测，其结果如图7所示。

可以看出，与现有的采用多尺度分析变化检测方法及采用模糊聚类变化检测方法相比，本发明采用的利用SURF特征点和支持向量机的变化检测方法，有效地抑制了噪声对变化检测精度的影响，增强了变化检测结果对噪声的鲁棒性，从而得到更为精确的变化检测结果。

将各变化检测方法对应的变化检测结果与参考变化图进行对比分析，分别统计错检(检测结果中被错误标记为变化的像素)的像素个数及漏检(检测结果中被错误标记为非变化的像素)的像素个数，最终计算出变化检测总体错误率P_TE，计算公式为：

$P_{TE}=\frac{FA+MA}{H\×W}\×100\%$

其中，FA、MA分别为变化检测结果中错检的像素个数及漏检的像素个数。

表1本发明方法与其它变化检测方法总体错误率的比较

表1列出了各变化检测方法所对应的总体错误率。本发明的方法对应的变化检测总体错误率最低，仅为4.15％；而采用多尺度分析以及采用模糊聚类的变化检测方法的总体错误率分别高达13.19％及27.23％；本发明方法能够非常有效地对多时相遥感图像进行变化检测，且其变化检测性能优于采用多尺度分析以及采用模糊聚类的变化检测方法的性能。

Claims (4)

1.一种多时相遥感图像变化检测方法，包括有以下步骤：

步骤1、输入两幅已经过预处理的多时相遥感图像I₁及I₂，所述预处理包括空间配准和辐射校正；并通过计算两幅图像间对应位置像素光谱值差值的绝对值构造差分图像I_d:

I_d(i,j)＝|I₁(i,j)-I₂(i,j)|

式中，i为图像行的像素，j为图像列的像素，1≤i≤H,1≤j≤W，H为遥感图像的行数，W为遥感图像的列数；

其特征是，还包括有

步骤2、分别从图像I₁及I₂中检测SURF特征点，得到特征点集合S₁及S₂；

步骤3、对集合S₁与集合S₂中的特征点进行匹配，得到集合S₁、S₂间相互匹配的特征点集合M₁、M₂，以及不能相互匹配的特征点集合R₁、R₂；集合M₁、M₂中的特征点一一对应、构成相互匹配的特征点对，其中S₁＝M₁∪R₁,S₂＝M₂∪R₂,

步骤4、提取非变化类像素样本

对集合M₁、M₂中的每个特征点，计算其在差分图像I_d中对应位置的像素的亮度值并归入非变化类像素样本集合S_u；

步骤5、提取变化类像素样本

将集合R₁、R₂合并为集合T，并计算T中每个像素点在差分图像I_d中对应位置的像素的亮度值并归入集合G；假定集合T中的像素由非变化类、未标记类及变化类像素构成，且每类像素的亮度均服从高斯分布，分别记为 则总体亮度直方图应为三维混合高斯分布：

$Z\left(g\right)=w_{1}N\left(g\right|{\mathrm{\μ}}_1,{\mathrm{\δ}}_1^2)+w_{2}N\left(g\right|{\mathrm{\μ}}_2,{\mathrm{\δ}}_2^2)+w_{3}N\left(g\right|{\mathrm{\μ}}_3,{\mathrm{\δ}}_3^2)$

式中，Z(g)表示三维混合高斯分布，g∈G为像素亮度值；w₁、w₂及w₃分别是非变化类、未标记类及变化类像素高斯分布在混合高斯中的权重；表示非变化类像素的亮度服从均值为μ₁、方差为的高斯分布；表示未标记类像素的亮度服从均值为μ₂、方差为的高斯分布；表示变化类像素的亮度服从均值为μ₃、方差为的高斯分布；

利用EM算法对上述模型的相关参数进行求解，最终将集合G中符合下式的像素亮度值归入变化类像素样本集合S_c：

S_c＝{g|g∈G且g≥μ₃-3δ₃}

μ₃、δ₃分别为变化类像素亮度的高斯分布的均值及标准差。

步骤6、依据步骤4所得的非变化类像素样本集合S_u和步骤5所得的变化类像素样本集合S_c构造训练集S_train、测试集S_test；再构造相应的训练集标签L_train、测试集标签L_test，其中集合L_train中的元素全为0且同集合S_train的大小相等，集合L_test中的元素全为1且同集合S_test的大小相等；最后在对测试集S_train、训练集S_test进行预处理的基础上训练一个支持向量机分类模型M_svm；

步骤7、依据步骤6所得的支持向量机分类模型M_svm对差分图像I_d的所有像素进行分类，得到变化及非变化区域，最终的变化检测结果将以二值变化检测掩膜的形式体现，其中值为0的像素表示未发生变化的类别，值为1的像素表示发生变化的类别。

2.根据权利要求1所述的多时相遥感图像变化检测方法，其特征是，所述的检测SURF特征点包括以下步骤：

步骤2.1、计算积分图像

遥感图像I₁或I₂中位于第x行第y列的像素点的积分图像值为：

$I_{t\mathrm{\Σ}}=\underset{i=0}{\overset{i\≤x}{\Σ}}\underset{j=0}{\overset{j\≤y}{\Σ}}{I}_t(x,y),t\∈\{1,2\}$

式中，t为遥感图像I₁或I₂的下角标代号；

步骤2.2、构建多尺度空间

初始盒子滤波器的模板大小为9×9，其对应的尺度大小为σ＝1.2，保持原始图像的大小不变，通过不断增加盒子滤波器的大小，并将其在原始图像的x、y及xy等三个方向上做卷积，最终得到多尺度空间；

步骤2.3、计算近似Hessian矩阵

对于遥感图像I₁和I₂中坐标为(x,y)的像素点X，按下式计算出该点在尺度σ上所对应的近似Hessian矩阵H₁(X,σ)和H₂(X,σ)：

$H_t(X,\mathrm{\σ})=\left[\begin{array}{cc}{D}_{xxt}(X,\mathrm{\σ})& {\mathrm{\λD}}_{xyt}\left(X\mathrm{\σ}\right)\\ {\mathrm{\λD}}_{xyt}(X,\mathrm{\σ})& D_{yyt}(X,\mathrm{\σ})\end{array}\right],t\∈\{1,2\}$

式中，t为遥感图像I₁或I₂的下角标代号；D_xxt(X,σ)是由步骤2.1得到的积分图像值与步骤2.2中所述的x方向上的盒子滤波器相乘求和的结果；D_xyt(X,σ)是由步骤2.1得到的积分图像值与步骤2.2中所述的xy方向上的盒子滤波器相乘求和的结果；D_yyt(X,σ)是由步骤2.1得到的积分图像值与步骤2.2中所述的y方向上的盒子滤波器相乘求和的结果；λ为权重系数，矩阵H_t(X,σ)对应行列式的值det(H_t)为：

det(H_t)＝D_xxtD_yyt-(λD_xyt)²

det(H_t)为像素点X在尺度σ上的响应值；

步骤2.4、确定极值点

H_t(X,σ)有两个特征值为θ_t1、θ_t2，

θ_t1×θ_t2＝det(H_t),t∈{1,2}

如果步骤2.3的det(H_t)小于零，则特征值θ_t1、θ_t2有不同的符号，该像素点不是极值点；反之，如果det(H_t)大于零，则特征值θ_t1、θ_t2同为正或负，该像素点是极值点；再将各极值点所对应的行列式的值与预设的阈值进行比较，所有小于该阈值的极值点都将被剔除，仅保留具有强响应值的极值点；

步骤2.5、建立SURF特征点的描述向量

按主方向确定一个以特征点为中心的正方形区域，将该正方形区域分成4×4个大小相等的子区域，每个子区域的边长为5σ，内有25个均匀分布的采样点，对每个子区域用尺度为2σ的Haar滤波器进行处理，得到x、y两个方向的响应dx、dy，计算每个子区域内的dx的响应值之和Σdx，|dx|的响应值之和Σ|dx|，dy的响应值之和Σdy，以及|dy|的响应值之和Σ|dy|，再将Σdx、Σ|dx|、Σdy及Σ|dy|构成每个子区域的特征描述向量：

V_subregion＝[Σd_xΣd_yΣ|d_x|Σ|d_y|]

V_subregion为该特征点的特征描述向量。

3.根据权利要求2所述的多时相遥感图像变化检测方法，其特征是，所述的集合S₁与集合S₂中的特征点进行匹配包括以下步骤：

步骤3.1、判断特征点之间特征相似性的高低

对集合S₁中的特征点P_1k，计算其与集合S₂中所有特征点P_2q的特征描述向量间的欧氏距离，

$Dis(P_{1k},P_{2q})={\[\underset{h=1}{\overset{64}{\Σ}}{(m_{1kh}-m_{2qh})}^2\]}^{1/2}$

式中，m_1kh表示集合S₁中的特征点P_1k的特征描述向量的第h个元素；m_2qh表示集合S₂中的特征点P_2q的特征描述向量的第h个元素；h为特征点的64维特征描述向量的元素序号；

欧氏距离值越大，相似度越低，反之，欧氏距离值越小，相似度越高；

步骤3.2、判断特征点之间是否满足8-邻接关系

假设在集合S₂中的特征点P_2q与集合S₁中的特征点P_1k最相似，根据P_1k与P_2q的坐标值，按下式判断二者在差分图像I_d中是否满足8-邻接关系：

$Match(P_{1k},P_{2q})=\left\{\begin{array}{cc}true,& \sqrt{{(x_{1k}-x_{2q})}^2+{(y_k-y_{2q})}^2}\≤\sqrt{2}\\ false,& \sqrt{{(x_{1k}-x_{2q})}^2+{(y_{1k}-y_{2q})}^2}>\sqrt{2}\end{array}\right.$

式中，(x_1k,y_1k)、(x_2q,y_2q)分别为特征点P_1k、P_2q的位置坐标；Match(P_1k,P_2q)用于度量特征点P_1k与特征点P_2q间是否满足8-邻接关系，若值为true，则P_1k与P_2q满足8-邻接关系，即相互匹配，否则，P_1k与P_2q不满足8-邻接关系，即不能相互匹配。

4.根据权利要求3所述的多时相遥感图像变化检测方法，其特征是，所述的对集合M₁、M₂中的每个特征点计算其在差分图像I_d中对应位置的像素的亮度值为：

对集合M₁、M₂中坐标为(x,y)的特征点，按下式计算其在差分图像I_d中对应位置的像素的亮度值g(x,y)：

当x或y不为整数时，采用双线性插值的方式来计算该特征点在差分图像I_d中对应位置的像素的亮度值，其中，x_u为x的小数部分，x_i为x的整数部分，y_u为y的小数部分，y_i为y的整数部分。