CN103150718B - 基于变化矢量分析与分类后比较的遥感图像变化检测方法 - Google Patents

Abstract

一种基于变化矢量分析与分类后比较的遥感图像变化检测方法，涉及遥感图像技术，包括：第一步，分别对两时相的遥感图像以像素为单位提取颜色和纹理特征；第二步，用马尔可夫随机场理论，将基于分类后比较方法的两时相遥感图像变化检测转化为两时相图像上马尔可夫能量函数比较；第三步，分别对两时相图像进行过分割，根据规则分别对两时相过分割图像区域进行重新调整；第四步，将基于变化矢量分析方法的两时相遥感图像变化检测转换为基于区域特征相似性度量函数项；第五步，将基于分类后比较的马尔可夫能量函数与基于变化矢量分析的相似性度量函数项进行联合，构造联合马尔可夫能量函数；第六步，用优化方法对联合马尔可夫能量函数进行求解，输出变化数据。

Description

基于变化矢量分析与分类后比较的遥感图像变化检测方法

技术领域

本发明涉及遥感图像处理技术领域，尤其是针对变化矢量分析方法和分类后比较方法相结合的两时相遥感图像变化检测方法。

背景技术

遥感图像变化检测是遥感图像处理领域中最重要的应用方向之一。遥感图像变化检测技术在自然灾害监测、国土资源规划管理、军事目标打击评估等众多军民领域都有重要的应用价值。随着遥感图像数据的分辨率增高，变化检测可以获取的变化信息日益丰富，使得变化检测的实际应用范围得到进一步扩大。

分类后比较方法(post-classificationcomparison(PCC))和变化矢量分析方法(changevectoranalysis(CVA))是两类常用的变化检测方法。分类后比较方法通常先对每个时相的图像进行分类，通过直接比较分类后图像的标记图，获取两时相上地面覆盖发生的变化情况。变化矢量分析方法则是直接获取关于两时相图像的变化矢量，通过对变化矢量的分析，得到两时相上地面覆盖发生的变化情况。分类后比较方法的优点是：不需要对两时相的图像进行辐射校正，能够比较方便地获取地面上覆盖类型的变化情况。其缺点是：忽略了两时相图像之间的依存关系，容易导致错误的累积，致使检测精度不高。变化矢量分析方法的优点是：比较充分地利用了两时相图像之间的依存关系，能够比较有效地抑制由于分类错误产生的误差累积。其缺点是：检测受光照、植被等辐射因素影响严重；另外，由于变化矢量获取的不唯一性，导致原始图像信息损失严重，致使检测结果不具备唯一性。

如何综合利用两种方法的优点，实现优势互补，提高检测精度，是一个比较有意义的问题。现有方法(J.Chen,X.Chen,X.Cui,andJ.Chen,"ChangeVectorAnalysisinPosteriorProbabilitySpace:ANewMethodforLandCoverChangeDetection",IEEETrans.Trans.Geosci.RemoteSens.Letters,vol.8,no.2,pp.317-321,Mar.2011.)试图利用变化矢量分析方法降低分类后比较方法的误差积累，但是这类方法的重要缺点就是不能充分利用两时相图像的依存信息，需要手工阈值选取，自动化程度不高。本发明综合利用分类后比较方法和变化矢量分析方法的优点，实现优势互补，第一次将分类后比较方法与变化矢量分析方法相结合的变化检测方法转化为一个联合马尔可夫能量函数表现形式，通过能量函数的优化求解，获取最优的变化检测输出结果。整个检测过程无需人工干预，自动化程度高。

发明内容

本发明的目的是提供一种基于变化矢量分析与分类后比较的遥感图像变化检测方法,以充分有效利用变化矢量分析方法与分类后比较方法的优势，提高变化检测精度。

为实现上述目的，本发明的技术方案是：

一种基于变化矢量分析与分类后比较的遥感图像变化检测方法，其包括以下步骤：

第一步，分别对两时相的遥感图像以像素为单位提取颜色和纹理特征；

第二步，用马尔可夫随机场理论，将基于分类后比较方法的两时相遥感图像变化检测转化为两时相图像上马尔可夫能量函数比较；

第三步，分别对两时相图像进行过分割，根据规则分别对两时相过分割图像区域进行重新调整；

第四步，将基于变化矢量分析方法的两时相遥感图像变化检测转换为基于区域特征的相似性度量函数项；

第五步，将基于分类后比较的马尔可夫能量函数与基于变化矢量分析的相似性度量函数项进行联合，构造一个联合马尔可夫能量函数；

第六步，用优化方法对联合马尔可夫能量函数进行求解，输出变化检测结果。

所述的遥感图像变化检测方法，其所述第一步，像素级特征提取，包括：

a1、以像素为单位，提取每个像素的CIELab颜色特征、Gabor纹理特征、熵特征；其中，Lab颜色特征维数为3；Gabor滤波器尺度参数、方向参数根据实际需要选择维数，Gabor纹理特征的维数为5×8＝40；根据图像分辨率选取窗口大小，计算以当前像素为中心的窗口区域内图像的熵作为该像素的熵特征，熵特征的维数为1；

a2、分别对每类特征进行归一化处理。

所述的遥感图像变化检测方法，其所述第二步，基于马尔可夫能量函数比较的分类后比较变化检测方法：利用马尔可夫随机场模型对每个时相的图像分别建模，以第一步获取的两时相遥感图像特征为基础，分别在每个时相图像上构造马尔可夫能量函数，将每个时相上图像的分类问题转化为马尔可夫能量函数优化问题，通过对两时相图像的马尔可夫能量函数比较，获取基于分类后比较方法的两时相遥感图像变化检测结果。

所述的遥感图像变化检测方法，其所述第二步，具体包括：

b1、根据最大后验概率估计理论和马尔可夫随机场理论，分别对每个时相的图像进行建模：

特征模型：假定每个时相图像中的每类特征均服从高斯分布，利用高斯分布函数，计算特征模型；对于每个时相图像中的任意一个像素p(i,j)，该像素属于第k类的高斯分布计算公式如下：

$p\left(x\right|y)=\frac{1}{\sqrt{2{\mathrm{\π\Σ}}_k}}\exp \[-\frac{1}{2}{(x-{\mathrm{\μ}}_k)}^T{\mathrm{\Σ}}_k^{-1}(x-{\mathrm{\μ}}_k)\]---\left(1\right)$

其中，p(x|y)为条件概率形式，x表示当前像素p(i,j)的特征矢量；y为当前像素p(i,j)对应的分类输出标记，y∈{1,...,k}；μ_k、Σ_k分别表示第k类的均值和方差。

先验模型：先验模型采取一阶Ising模型，即仅认为当前像素p(i,j)只与其一阶邻域内N的像素存在相互作用，而与其他像素之间不存在相互关系；马尔可夫一阶Ising模型表述为：

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\∈N}{\Σ}{\mathrm{\β}}_1\·\mathrm{\δ}\left(y,y_m\right))---\left(2\right)$

其中，p(y)为先验概率形式，δ(y,y_m)定义为卡迪拉克函数：如果y＝y_m，则δ(y,y_m)＝1，否则δ(y,y_m)＝0；β₁为平滑加权系数，控制邻域像素之间的互作用的大小；N为一阶邻域；Z为规整化常数；

b2、马尔可夫能量函数构造：

每个时相的图像上马尔可夫能量即为图像上每个像素p(i,j)的能量之和，其函数表达形式为：

$\arg \underset{x\∈X}{min}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(3\right)$

b3、对两时相图像的马尔可夫能量函数进行比较，获取分类后比较方法的马尔可夫能量函数表述形式：由于两时相的图像是独立不相关的，因此每个时相图像的分类是独立的，对应两时相图像的分类后比较就是两时相图像的马尔可夫能量函数分别最小化；直接将两时相图像上定义的能量函数相加，作为两时相图像分类后比较的能量函数表述形式：

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(4\right)$

所述的遥感图像变化检测方法，其所述第三步，过分割图像区域调整：分别对两时相的遥感图像进行过分割处理，获取区域为单位的同质区域图像，按照规则分别对两个过分割图像进行区域调整，包括：

c1、利用分水岭分割方法分别对每个时相的遥感图像进行过分割处理，每个过分割区域即为同质区域；其中，当某个区域的总像素个数少于200时该区域停止分割；

c2、根据两时相图像的过分割图，对每个时相图像的过分割区域逐一进行比较，并根据比较准则对每个时相上的过分割区域进行调整；其中，假定区域A为时相1图像中的一个过分割区域，区域B为区域A在第二个时相图像上的对应区域，则区域A和区域B比较准则为：

如果A＝B，则区域A和区域B均保持不变；

如果A≠B，且A∩B＝A，则区域A保持不变，将区域B分割为

B＝B1∪B2，其中，B1为区域A在第二个时相图像上的对应区域，B2为B-B1；

如果A∪B>A，且A∪B>B，则将区域A分割为A-(A∩B)和A∩B两部分，区域B分割为B-(A∩B)和A∩B两部分。

所述的遥感图像变化检测方法，其所述第四步，基于区域相似性度量变化矢量分析的变化检测方法：对两时相图像的特征矢量进行比较分析，计算特征矢量之间的相似性，将基于区域为单位的遥感图像变化矢量分析转化为特征矢量的相似性度量函数项，包括：

d1、两时相图像中对应像素p1和p2之间的相似性度量dism_p定义为两个像素的特征欧式距离加权和：像素p1和p2的Lab颜色特征、Gabor纹理特征和熵特征的欧式距离加权和，其函数表示形式为：

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²(5)其中，x1_Lab，x1_Gabor和x1_entory分别表示像素p1的Lab颜色特征矢量、Gabor纹理特征和熵特征；x2_Lab，x2_Gabor和x2_entory分别表示像素p2的Lab颜色特征矢量、Gabor纹理特征和熵特征；ω₁和ω₂为加权系数，为常数，一般ω₁∈[0,1]，ω₂∈[0,1]；

d2、两时相图像中任意两个区域之间的相似性度量dism_r即为区域中所有像素相似性度量之和；其计算公式为：

dism_r＝∑dism_p(p1,p2)(6)

所述的遥感图像变化检测方法，其所述第五步，联合马尔可夫能量函数构造：将变化矢量分析的能量项与分类后比较的能量函数进行联合，获取联合马尔可夫能量函数；整个函数一方面兼顾分类后比较方法有最优解，另一方面兼顾两时相图像的相似性最大，使得两时相图像之间的变化得到最优检测，包括：

将基于变化矢量分析方法的能量项与基于分类后比较方法的能量函数相加，得到基于变化矢量分析方法与分类后比较方法的联合能量表达形式：

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\mathrm{\Σ}}}\mathrm{\Σ}\{\frac{1}{2}\ln \left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\mathrm{\Σ}}\mathrm{\δ}(y,y_m)+{\mathrm{\β}}_2{\left({\mathrm{dism}}_p\right)}^{\mathrm{\β}}\}---\left(7\right)$

其中，β₂为加权系数，用于控制变化矢量分析方法的能量项与分类后比较方法的能量函数的比重，当β₂取较小值时，分类后比较方法起主导作用；当β₂取较大值时，变化矢量分析方法起主导作用；β为控制函数，一般为开关函数形式，如果两时相图像中像素对应的类标值相同，则β的取值为1；如果两时相图像中像素对应的类标值不同，则β的取值为-1。

所述的遥感图像变化检测方法，其所述第六步，利用优化方法对联合马尔可夫能量函数进行求解，包括：

f1、对两时相图像中对应坐标的每个像素逐点迭代，计算每个像素的最小联合能量；对图像中所有像素的最小联合能量求和，作为当次迭代的能量；

f2、记录当前最小联合能量对应的变化检测输出类标，即为当次最优输出结果；

f3、迭代：当本次迭代计算所得的最小联合能量小于前次迭代最小联合能量时，将本次联合马尔可夫能量函数对应的格局，即当前两时相图像变化检测输出类标图作为输出结果；否则，将上次迭代对应的变化检测输出类标图作为输出结果；

f4、同时，当迭代次数满足一个事先给定的最大值，或者当连续5次的迭代输出结果保持不变时，停止迭代，输出变化数据。

所述的遥感图像变化检测方法，其所述f4步骤中迭代次数的最大值，即输出变化数据迭代次数上限值，设置为200；迭代次数事先给定的最大值，根据精度与速度要求折中选取。

所述的遥感图像变化检测方法，其所述第六步中迭代结束，能量函数收敛于局部最优值时，对应的图像格局即为最终的变化检测结果。

本发明针对两时相遥感图像提出了一种变化矢量分析与分类后比较相结合的变化检测方法，通过合理构造马尔可夫随机场模型能量函数，将基于分类后比较方法的两时相遥感图像变化检测和基于变化矢量分析方法的两时相遥感图像变化检测分别转化为能量函数表述形式，并将这两个能量函数进行巧妙地联合，将变化矢量分析与分类后比较相结合的两时相遥感图像变化检测转化为一个联合马尔可夫能量函数优化求解问题，利用迭代条件模型对联合马尔可夫能量函数进行优化求解，得到变化检测结果。该方法有效地利用了变化矢量分析方法与分类后比较方法在变化检测中的优势，通过优势互补，提高了变化检测的精度，使得该方法有较好的鲁棒性；同时将变化矢量分析与分类后比较相结合的两时相遥感图像变化检测转化为联合马尔可夫能量函数优化求解问题，提高了变化检测的自动化程度，具有较好实用性。

本发明兼具基于变化矢量分析变化检测方法的优点和基于分类后比较变化检测方法的优点，第一次将分类后比较方法与变化矢量分析方法相结合的变化检测方法转化为一个联合马尔可夫能量函数表现形式，通过能量函数的优化求解，获取最优的变化检测输出结果，检测过程无需人工干预，自动化程度高。

附图说明

图1为本发明的一种基于变化矢量分析与分类后比较的遥感图像变化检测方法流程示意图；

图2为本发明方法第二步基于分类后比较方法的两时相遥感图像变化检测马尔可夫能量函数构造示意图；

图3为本发明方法第三步中过分割图像区域调整准则示意图；

图4为本发明方法第五步联合马尔可夫能量函数构造示意图；

图5为本发明方法应用实例数据；其中：

图5(a)为2002年拍摄图像；

图5(b)为2003年拍摄图像；

图6为本发明的一种基于变化矢量分析与分类后比较的遥感图像变化检测方法检测结果示例；其中：

图6(a)为基于分类后比较方法的变化检测结果；

图6(b)是基于变化矢量分析方法的变化检测结果；

图6(c)是本发明方法中变化矢量分析与分类后比较相结合的变化检测结果；

图6(d)为人工检测结果。

具体实施方式

本发明的一种基于变化矢量分析与分类后比较的遥感图像变化检测方法，将基于分类后比较方法和基于变化矢量分析方法的两时相遥感图像变化检测，分别转化为两时相图像的马尔可夫能函数的比较和相似性度量能量项，通过联合两种方法的能量表述函数，构造一个联合马尔可夫能量函数，利用现有的能量函数优化求解方法，对联合马尔可夫能量函数进行求解，以得到最优的变化检测输出结果。

第一步，分别提取两时相的遥感图像基础特征。方法是：

1.1、以像素为单位，提取每个像素的Lab颜色特征：

1.1.1图像平滑滤波：图像中的每个像素p(i,j)，以坐标(i,j)为中心，取3*3大小(窗口大小可根据需要调整)的图像区域，计算区域中像素的光谱均值作为当前像素p(i,j)的光谱值。

1.1.2RGB颜色空间到CIELab颜色空间的转换计算公式为：

$X=\frac{0.412453*R+0.357580*G+0.180423*B}{0.950456}$

$X=\frac{0.412453*R+0.357580*G+0.180423*B}{0.950456}$

$Y=\frac{0.212671*R+0.715160*G+0.072169*B}{1}$

$Z=\frac{0.019334*R+0.119193*G+0.950227*B}{1.088754}$

XT＝X>th,YT＝Y>th,ZT＝Z>th,th＝0.008856

fX＝XT.*X^1/3+(～XT).*(7.787.*X+16/116)(8)

fY＝YT.*Y^1/3+(～YT).*(7.787.*Y+16/116)

fZ＝ZT.*Z^1/3+(～ZT).*(7.787.*Z+16/116)

L＝YT.*(116*Y^1/3-16.0)+(～YT).*(903.3.*Y)

a＝500*(fX-fY)

b＝200*(fY-fZ)

1.1.3图像的Lab颜色特征维数为3，按图像坐标方式存储Lab颜色特征。

1.2、提取图像的Gabor纹理特征：

1.2.1Gabor滤波器：

x'＝xcosθ+ysinθ(9)

y'＝-xsinθ+ycosθ

其中，λ为波长，θ为Gabor滤波器的条带方向，为Gabor滤波器的相位参数，γ为长宽比，控制Gabor的椭球率。

Gabor能量函数定义如下形式：

$e_{\mathrm{\λ},\mathrm{\θ}}(x,y)=\sqrt{{\mathrm{\γ}}_{\mathrm{\λ},\mathrm{\θ},0}^2(x,y)+{\mathrm{\γ}}_{\mathrm{\λ},\mathrm{\θ},-\frac{\mathrm{\π}}{2}}^2(x,y)}---\left(10\right)$

其中，和分别为g_λ,θ,0(x,y)和与灰度图像卷积的结果。

1.2.2将原始RGB图像转化为灰度图像：直接将图像R、G、B三个通道上每个像素的光谱值相加求均值，作为灰度图像的光谱值。

1.2.3设定Gabor滤波器的尺度参数为5，方向参数为8，可以构造40中Gabor滤波器，使用每个滤波器对灰度图像滤波，即将每个Gabor滤波器与灰度图像进行卷积；Gabor纹理特征的维数为40；

1.2.4按图像坐标方式存储Gabor纹理特征；

1.3、以像素坐标为单位，提取每个像素的熵特征：

1.3.1将原始RGB图像转化为灰度图像：直接将图像R、G、B三个通道上每个像素的光谱值相加求均值，作为灰度图像的光谱值。

1.3.2图像中的每个像素p(i,j)，以坐标(i,j)为中心，取10*10大小的图像区域，计算图像区域的熵，作为该像素的熵特征。熵的计算公式如下：

$entry=\underset{i=1}{\overset{255}{\Σ}}{p}_i\log p_i---\left(11\right)$

其中，p_i表示图像区域中灰度值为i的像素个数占整个图像区域中像素个数的比例。

1.3.3图像的熵特征的维数为1，按图像坐标方式存储熵特征。

1.4、分别对每类特征归一化处理。归一化计算公式为：

f＝(f-min(min(f)))/(max(max(f))-min(min(f)))(12)

第二步，利用马尔可夫随机场模型对每个时相的图像分别建模，以获取的两时相的遥感图像特征矢量为基础，分别在每个时相图像上构造马尔可夫能量函数，对两时相上的马尔可夫能量函数进行比较，获取分类后比较的马尔可夫能量函数表述形式。方法是：

2.1、以两时相的图像为基础，分别在每个时相的图像上构造一个马尔可夫随机场模型G(V,E)：每个时相图像中每一个像素p(i,j)表示为马尔可夫随机场模型中的一个结点V，每个节点V与其周围坐标欧式距离差值为1的四个像素构成该节点V的一个一阶邻域N，领域中垂直、水平4个方向上的像素之间的特征矢量差作为马尔可夫随机场模型的边E。

2.2、以每个时相图像提取的特征矢量为基础，分别对每个时相的遥感图像进行分类。分类的原则为：图像中相似性较高的像素聚为一类，以保证类内方差小于类间方差，且该比值越小越好。

2.2.1依据最大后验概率估计理论，给定某一个时相图像的特征矢量X，其分类输出Y即为最大后验概率求解问题：

p(Y|X)＝p(X|Y)p(Y)/p(X)(13)

2.2.2对应到图像上的每一个像素p(i,j)，其分类过程的数学描述形式表述为：

$y=\arg \underset{x\∈X}{max}\left\{p\left(x\right|y)p\left(y\right)\right\}---\left(14\right)$

其中，x为当前像素p(i,j)的特征矢量，y为该点对应的分类类标；p(x|y)描述图像特征与分类输出类标之间的分布关系；p(y)描述图像的输出类标之间的相互作用关系。

2.3依据马尔可夫随机场模型-最大后验概率估计理论，在每个时相的图像上构造特征模型p(x|y)和先验模型p(y)：

2.3.1特征模型：假定每个时相图像中的每类特征均服从高斯分布，利用高斯分布函数，计算特征模型。对于每个时相图像中的任意一个像素p(i,j)，该像素属于第k类的高斯分布函数计算公式如下：

$p\left(x\right|y)=\frac{1}{\sqrt{2{\mathrm{\π\Σ}}_k}}\exp \[-\frac{1}{2}{(x-{\mathrm{\μ}}_k)}^T{\mathrm{\Σ}}_k^{-1}(x-{\mathrm{\μ}}_k)\]---\left(15\right)$

其中，x表示当前像素p(i,j)的特征变化矢量；y为当前像素p(i,j)对应的分类输出标记，y∈{1,...,k}；μ_k、Σ_k分别表示第k类特征矢量的均值和方差。

2.3.2先验模型：先验模型采取一阶Ising模型，即仅认为当前像素p(i,j)只与其一阶邻域内N的像素存在相互作用，而与其他像素之间不存在相互作用关系。马尔可夫一阶Ising模型表述为：

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\∈N}{\Σ}{\mathrm{\β}}_1\·\mathrm{\δ}\left(y,y_m\right))---\left(16\right)$

其中，δ(y,y_m)定义为卡迪拉克函数，如果y＝y_m，则δ(y,y_m)＝1，否则δ(y,y_m)＝0；β₁为平滑加权系数，控制邻域像素之间的作用大小；N为一阶邻域；Z为规整化常数。

2.4分别在每个时相图像上构造马尔可夫能量函数：每个时相的图像上马尔可夫能量即为图像上每个像素p(i,j)的能量之和。对后验概率表达式两边分别取log，则最大后验概率等价为如下的函数表达：

$\arg \underset{x\∈X}{min}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(17\right)$

当图像中每一点的类标值最适合时，整个能量函数的值趋于稳定，具有最小值，此时对应的马尔可夫能量格局即为最优的分类输出结果。

2.5、由于两时相的图像是独立不相关的，因此每个时相图像的分类是独立不相关的。当每个时相的图像上的马尔可夫能量函数最小时，即对应每个时相上图像的最优输出类标图。因此，对应两时相图像的分类后比较就等价于两时相图像的马尔可夫能量函数之和的最小化。采取直接将两时相图像上定义的能量函数进行相加的处理，作为两时相图像分类后比较的能量函数表述形式：

$\arg \underset{x\∈X}{min}\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(18\right)$

第三步，分别对两时相的遥感图像进行过分割处理，获取区域为单位的同质区域图像，按照规则对两个过分割图像进行区域比较，并分别对两时相过分割图像区域进行重新调整。方法是：

3.1、利用分水岭分割方法(J.RoerdinkandA.Meijster,“TheWatershedTransform:Definitions,Algorithms,andParallellizationStrategies,”FundamentaInformaticae,vol.41,pp.187-228,2000.)，以每个时相图像的CIELab颜色为基础，分别对每个时相的遥感图像进行过分割处理。每个过分割区域具备一定的区域同质一致性，即为同质区域。在过分割过程中，为了保证区域的最大一致性以及变化检测的实际解译性，选取区域面积作为约束条件，即当某个分割区域的总像素个数少于200时该区域停止分割；

3.2、根据两时相图像的过分割图，对每个时相图像的过分割区域逐一进行比较，并根据比较准则对每个时相上的过分割区域进行调整。其中，假定区域A为时相1图像中的一个过分割区域，区域B为区域A在第二个时相图像上的对应区域，则区域A和区域B比较准则为：

如果A＝B，则区域A和区域B均保持不变；

如果A≠B，且A∩B＝A，则区域A保持不变，将区域B分割为

B＝B1∪B2，其中，B1为区域A在第二个时相图像上的对应区域，B2为B-B1；

如果A∪B>A，且A∪B>B，则将区域A分割为A-(A∩B)和A∩B两部分，区域B分割为B-(A∩B)和A∩B两部分；

第四步，对两时相图像的特征矢量进行比较分析，计算特征矢量之间的相似性，将基于区域为单位的遥感图像变化矢量分析转化为特征矢量的相似性度量。方法是：

4.1、像素p1和p2的相似性度量函数dism_p：两时相图像中任意像素p1和p2之间的相似性度量dism_p定义为两个像素特征的欧式距离加权和，其函数表示形式为：

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²，(19)

其中，x1_Lab，x1_Gabor和x1_entory分别表示像素p1的Lab颜色特征矢量、Gabor纹理特征和熵特征；x2_Lab，x2_Gabor和x2_entory分别表示像素p2的Lab颜色特征矢量、Gabor纹理特征和熵特征；ω₁和ω₂为加权系数，为常数，一般ω₁∈[0,1]，ω₂∈[0,1]。通常根据图像自身特点，选取ω₁和ω₂的值。本发明中ω₁和ω₂取值均为1。

4.2、两时相图像中任意两个区域之间的相似性度量dism_r即为区域中所有像素相似性度量之和。其计算公式为：

dism_r＝∑dism_p(x1,x2)(20)

第五步，将变化矢量进行分析的能量项与分类后比较的能量函数进行联合，获取联合马尔可夫能量函数。使得整个函数一方面保证分类后比较方法有最优解，另一方面保证两时相图像的相似性最大，使得两时相图像之间的变化得到最优检测。具体方法是将变化矢量分析方法的相似性度量函数项作为分类后比较方法能量函数的有偏项，构造一个新的能量函数表达式：

5.1、将变化矢量分析方法的能量项与分类后比较的能量函数进行相加，得到变化矢量分析方法与分类后比较方法的联合马尔可夫能量函数表达形式：

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)+{\mathrm{\β}}_2{\left({\mathrm{dism}}_p\right)}^{\mathrm{\β}}\}---\left(21\right)$

其中，β₂为加权系数，用于控制变化矢量分析方法的能量项与分类后比较方法的能量函数的比重，当β₂取较小值时，分类后比较方法起主导作用；当β₂取较大值时，变化矢量分析方法起主导作用；β为控制函数，一般为开关函数形式，如果两时相图像中像素对应的类标值相同，则β的取值为1；如果两时相图像中像素对应的类标值不同，则β的取值为-1。

5.2、本发明中，β₁取常数1，β₂取常数0.5。

第六步，利用优化方法对联合马尔可夫能量函数进行求解。方法是：

6.1、对能量函数(21)进行求解：采用迭代条件模型(iterativeconditionalmode(ICM))对联合马尔可夫能量函数(21)进行优化求解。迭代次数一般在200次左右，最终收敛于局部最优值。

6.1.1假定分类后比较方法中，每个时相的图像分类的类别数为k，k的值一般直接给定。

6.1.2对两时相图像中的每一个像素，分别计算每个时相图像类标值为[1,2,…,k]时马尔可夫能量，即计算类标值分别为[1,2,…,k]时函数(17)的值；

6.1.3计算两时相图像对应的像素之间的相似性度量，即计算函数项(19)的值；

6.1.4根据两时相像素的类标值，选取开关函数β的取值；将两时相图像中像素的分类后比较方法的马尔可夫能量与变化矢量分析方法的马尔可夫能量相加，得到当前像素上联合马尔可夫能量，即函数(21)中每一个像素的能量；

6.1.5计算图像中所有像素的联合马尔可夫能量，进行求和，得到两时相图像的联合马尔可夫能量，即函数(21)对应的函数值；

6.1.6记录当前函数(21)输出能量值对应的两时相图像的类标图；

6.1.7对两时相的类标图逐个像素进行比较：如果两时相图像对应的像素类标值相同，则该像素位置对应的变化检测输出标记为0；否则，该像素位置输出标记为1；

6.1.8迭代：当本次迭代计算所得的联合马尔可夫能量函数(21)的值小于前次迭代值时，将本次联合马尔可夫能量函数对应的格局，即当前两时相图像变化检测输出类标图作为输出结果；否则，将上次迭代对应的变化检测输出类标图作为输出结果；

6.1.9当迭代次数满足一个事先给定的最大值，或者当连续5次的迭代输出结果保持不变时，停止迭代。本发明方法中，迭代次数上限值设置为200(迭代次数可以根据精度与速度要求折中选取)。

6.1.10迭代结束，能量函数收敛于局部最优值时对应的图像格局即为最终的变化检测结果。

以下进一步说明本发明的应用实例：

实例给出了一组典型城市场景QuickBird多光谱图像，图像拍摄地区为北京市的部分区域，拍摄时间分别为2002年和2003年。图像由全色图像(0.6米/像素)和多光谱图像(2.4米/像素，红、绿、兰、近红外四个波段)组成，为了同时具备高的空间分辨率和光谱分辨率，将全色图像和多光谱图像进行了Ehlers融合(M.Ehlers.SpectralcharacteristicspreservingimagefusionbasedonFourierDomainfiltering,InProc.SPIE,5574:1-4,Bellingham,2004.)，融合后图像的大小为1024×1024像素。图像中的主要变化为：部分荒地上建起了大片的房屋，原有的建筑区域成了绿地。

以人工检测结果作为参考结果，对比三种变化检测方法的准确性，一种为基于分类后比较方法的变化检测方法，一种是基于变化矢量分析方法的变化检测方法，另一种是本发明实例方法。其中，基于分类后比较方法的变化检测方法即为本发明具体实施步骤中第二步的输出结果；基于变化矢量分析方法的变化检测方法为利用马尔可夫随机场对变化矢量进行直接分析(L.BruzzoneandD.F.Prieto,“Automaticanalysisofthedifferenceimageforunsupervisedchangedetection,”IEEETrans.Geosci.RemoteSens.vol.38,no.3,pp.1171–1182,May2000.)

从四个方面对变化检测结果进行评价(见表1)：1)误检率；2)漏检率；3)错误率；4)kappa系数。

表1变化检测结果量化比较结果

量化比较结果说明，本发明的基于变化矢量分析与分类后比较的遥感图像变化检测方法检测精度有较大提高。

Claims (8)

1.一种基于变化矢量分析与分类后比较的遥感图像变化检测方法，其特征在于，包括以下步骤：

第一步，分别对两时相的遥感图像以像素为单位提取颜色和纹理特征；

第二步，用马尔可夫随机场理论，将基于分类后比较方法的两时相遥感图像变化检测转化为两时相图像上马尔可夫能量函数比较；

第三步，分别对两时相图像进行过分割，根据规则分别对两时相过分割图像区域进行重新调整；

第四步，将基于变化矢量分析方法的两时相遥感图像变化检测转换为基于区域特征的相似性度量函数项；

第五步，将基于分类后比较的马尔可夫能量函数与基于变化矢量分析的相似性度量函数项进行联合，构造一个联合马尔可夫能量函数；

第六步，用优化方法对联合马尔可夫能量函数进行求解，输出变化检测结果；

其中，第二步中基于马尔可夫能量函数比较的分类后比较变化检测方法：利用马尔可夫随机场模型对每个时相的图像分别建模，以第一步获取的两时相遥感图像特征为基础，分别在每个时相图像上构造马尔可夫能量函数，将每个时相上图像的分类问题转化为马尔可夫能量函数优化问题，通过对两时相图像的马尔可夫能量函数比较，获取基于分类后比较方法的两时相遥感图像变化检测结果；

其中，第三步中，分别对两时相图像进行过分割，根据规则分别对两时相过分割图像区域进行重新调整包括：

c1、利用分水岭分割方法分别对每个时相的遥感图像进行过分割处理，每个过分割区域即为同质区域；其中，当某个区域的总像素个数少于200时该区域停止分割；

c2、根据两时相图像的过分割图，对每个时相图像的过分割区域逐一进行比较，并根据比较准则对每个时相上的过分割区域进行调整；其中，假定区域A为时相1图像中的一个过分割区域，区域B为区域A在第二个时相图像上的对应区域，则区域A和区域B比较准则为：

如果A＝B，则区域A和区域B均保持不变；

如果A≠B，且A∩B＝A，则区域A保持不变，将区域B分割为B＝B1∪B2，其中，B1为区域A在第二个时相图像上的对应区域，B2为B-B1；

如果A∪B>A，且A∪B>B，则将区域A分割为A-(A∩B)和A∩B两部分，区域B分割为B-(A∩B)和A∩B两部分。

2.如权利要求1所述的遥感图像变化检测方法，其特征在于，所述第一步，像素级特征提取，包括：

a1、以像素为单位，提取每个像素的CIELab颜色空间中Lab颜色特征、Gabor纹理特征、熵特征；其中，Lab颜色特征维数为3；Gabor滤波器尺度参数、方向参数根据实际需要选择维数，Gabor纹理特征的维数为5×8＝40；根据图像分辨率选取窗口大小，计算以当前像素为中心的窗口区域内图像的熵作为该像素的熵特征，熵特征的维数为1；

a2、分别对每类特征进行归一化处理。

3.如权利要求1所述的遥感图像变化检测方法，其特征在于，所述第二步，具体包括：

b1、根据最大后验概率估计理论和马尔可夫随机场理论，分别对每个时相的图像进行建模：

特征模型：假定每个时相图像中的每类特征均服从高斯分布，利用高斯分布函数，计算特征模型；对于每个时相图像中的任意一个像素p(i,j)，该像素属于第k类的高斯分布计算公式如下：

其中，p(x|y)为条件概率形式，x表示当前像素p(i,j)的特征矢量；y为当前像素p(i,j)对应的分类输出标记，y∈{1,...,k}；μ_k、Σ_k分别表示第k类的均值和方差；

先验模型：先验模型采取一阶Ising模型，即仅认为当前像素p(i,j)只与其一阶邻域内N的像素存在相互作用，而与其他像素之间不存在相互关系；马尔可夫一阶Ising模型表述为：

其中，p(y)为先验概率形式，δ(y,y_m)定义为卡迪拉克函数：如果y＝y_m，则δ(y,y_m)＝1，否则δ(y,y_m)＝0；β₁为平滑加权系数，控制邻域像素之间的互作用的大小；N为一阶邻域；Z为规整化常数；

b2、马尔可夫能量函数构造：

每个时相的图像上马尔可夫能量即为图像上每个像素p(i,j)的能量之和，其函数表达形式为：

其中，X为时相图像的特征矢量；

b3、对两时相图像的马尔可夫能量函数进行比较，获取分类后比较方法的马尔可夫能量函数表述形式：由于两时相的图像是独立不相关的，因此每个时相图像的分类是独立的，对应两时相图像的分类后比较就是两时相图像的马尔可夫能量函数分别最小化；直接将两时相图像上定义的能量函数相加，作为两时相图像分类后比较的能量函数表述形式：

4.如权利要求3所述的遥感图像变化检测方法，其特征在于，所述第四步，基于区域相似性度量变化矢量分析的变化检测方法：对两时相图像的特征矢量进行比较分析，计算特征矢量之间的相似性，将基于区域为单位的遥感图像变化矢量分析转化为特征矢量的相似性度量函数项，包括：

d1、两时相图像中对应像素p1和p2之间的相似性度量dism_p定义为两个像素的特征欧式距离加权和：像素p1和p2的Lab颜色特征、Gabor纹理特征和熵特征的欧式距离加权和，其函数表示形式为：

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²，(5)

其中，x1_Lab，x1_Gabor和x1_entory分别表示像素p1的Lab颜色特征矢量、Gabor纹理特征和熵特征；x2_Lab，x2_Gabor和x2_entory分别表示像素p2的Lab颜色特征矢量、Gabor纹理特征和熵特征；ω₁和ω₂为加权系数，为常数，ω₁∈[0,1]，ω₂∈[0,1]；

d2、两时相图像中任意两个区域之间的相似性度量dism_r即为区域中所有像素相似性度量之和；其计算公式为：

dism_r＝∑dism_p(p1,p2)(6)。

5.如权利要求4所述的遥感图像变化检测方法，其特征在于，所述第五步，联合马尔可夫能量函数构造：将变化矢量分析的能量项与分类后比较的能量函数进行联合，获取联合马尔可夫能量函数；整个函数一方面兼顾分类后比较方法有最优解，另一方面兼顾两时相图像的相似性最大，使得两时相图像之间的变化得到最优检测，包括：

将基于变化矢量分析方法的能量项与基于分类后比较方法的能量函数相加，得到基于变化矢量分析方法与分类后比较方法的联合马尔可夫能量表达形式：

其中，β₂为加权系数，用于控制变化矢量分析方法的能量项与分类后比较方法的能量函数的比重，当β₂取较小值时，分类后比较方法起主导作用；当β₂取较大值时，变化矢量分析方法起主导作用；β为控制函数，为开关函数形式，如果两时相图像中像素对应的类标值相同，则β的取值为1；如果两时相图像中像素对应的类标值不同，则β的取值为-1。

6.如权利要求1所述的遥感图像变化检测方法，其特征在于，所述第六步，利用优化方法对联合马尔可夫能量函数进行求解，包括：

f1、对两时相图像中对应坐标的每个像素逐点迭代，计算每个像素的最小联合能量；对图像中所有像素的最小联合能量求和，作为当次迭代的能量；

f2、记录当前最小联合能量对应的变化检测输出类标，即为当次最优输出结果；

f3、迭代：当本次迭代计算所得的最小联合能量小于前次迭代最小联合能量时，将本次联合马尔可夫能量函数对应的格局，即当前两时相图像变化检测输出类标图作为输出结果；否则，将上次迭代对应的变化检测输出类标图作为输出结果；

f4、同时，当迭代次数满足一个事先给定的最大值，或者当连续5次的迭代输出结果保持不变时，停止迭代，输出变化数据。

7.如权利要求6所述的遥感图像变化检测方法，其特征在于，所述f4步骤中迭代次数的最大值，即输出变化数据迭代次数上限值，设置为200；迭代次数事先给定的最大值，根据精度与速度要求折中选取。

8.如权利要求6所述的遥感图像变化检测方法，其特征在于，所述第六步中迭代结束，能量函数收敛于局部最优值时，对应的图像格局即为最终的变化检测结果。