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

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$p\left(x\right|y)=\frac{1}{\sqrt{2\pi {\Sigma }_k}}\exp [-\frac{1}{2}{(x-{\mu }_k)}^T{\Sigma }_k^{-1}(x-{\mu }_k)]---\left(1\right)$

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

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

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\in N}{\Sigma }{\beta }_1·\delta (y,y_m))---\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\in X}{\min }\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (y,y_m)\}---\left(3\right)$

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

$\arg \underset{x\in X}{\min }\underset{s=1}{\overset{2}{\Sigma }}\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (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\in X}{\min }\underset{s=1}{\overset{2}{\Sigma }}\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (y,y_m)+{\beta }_2{\left({\mathrm{dism}}_p\right)}^{\beta }\}---\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.2 RGB颜色空间到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.1 Gabor滤波器：

x′＝xcosθ+ysinθ                        (9)

y′＝-xsinθ+ycosθ

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

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

$e_{\lambda ,\theta }(x,y)=\sqrt{{\gamma }_{\lambda ,\theta ,0}^2(x,y)+{\gamma }_{\lambda ,\theta ,-\frac{\pi }{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 大小的图像区域，计算图像区域的熵，作为该像素的熵特征。熵的计 算公式如下：

$\mathrm{entry}=\underset{i=1}{\overset{255}{\Sigma }}{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\in 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(x/y)：假定每个时相图像中的每类特征均服 从高斯分布，利用高斯分布函数，计算特征模型。对于每个时相图像 中的任意一个像素p(i，j)，该像素属于第k类的高斯分布函数计算 公式如下：

$p\left(x\right|y)=\frac{1}{\sqrt{2\pi {\Sigma }_k}}\exp [-\frac{1}{2}{(x-{\mu }_k)}^T{\Sigma }_k^{-1}(x-{\mu }_k)]---\left(15\right)$

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

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

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\in N}{\Sigma }{\beta }_1·\delta (y,y_m))---\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\in X}{\min }\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (y,y_m)\}---\left(17\right)$

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

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

$\arg \underset{x\in X}{\min }\underset{s=1}{\overset{2}{\Sigma }}\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (y,y_m)\}---\left(18\right)$

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

3.1、利用分水岭分割方法(J.Roerdink and A.Meijster，“The  Watershed Transform：Definitions，Algorithms，and  Parallellization Strategies，”Fundamenta Informaticae，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\in X}{\min }\underset{s=1}{\overset{2}{\Sigma }}\Sigma \{\frac{1}{2}\ln \left(2\pi {\Sigma }_i\right)+\frac{{(x-{\mu }_k)}^2}{2{\Sigma }_k}+{\beta }_1\underset{y_m\in N}{\Sigma }\delta (y,y_m)+{\beta }_2{\left({\mathrm{dism}}_p\right)}^{\beta }\}---\left(21\right)$

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

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

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

6.1、对能量函数(21)进行求解：采用迭代条件模型(iterative  conditional mode(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.Spectral  characteristics preserving image fusion based on Fourier  Domain filtering，In Proc.SPIE，5574：1-4，Bellingham，2004.)， 融合后图像的大小为1024×1024像素。图像中的主要变化为：部分 荒地上建起了大片的房屋，原有的建筑区域成了绿地。

以人工检测结果作为参考结果，对比三种变化检测方法的准确 性，一种为基于分类后比较方法的变化检测方法，一种是基于变化矢 量分析方法的变化检测方法，另一种是本发明实例方法。其中，基于 分类后比较方法的变化检测方法即为本发明具体实施步骤中第二步 的输出结果；基于变化矢量分析方法的变化检测方法为利用马尔可夫 随机场对变化矢量进行直接分析(L.Bruzzone and D.F.Prieto， “Automatic analysis of the difference image for unsupervised  change detection，”IEEE Trans.Geosci.Remote Sens.vol.38， no.3，pp.1171-1182，May 2000.)

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

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

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