顾及空间引力的变权重MRF高分遥感影像变化检测

韦春桃，李倩倩，卢志豪，张冬梅，陈奕州

(重庆交通大学 智慧城市学院，重庆 400074)

0 引言

在城市化高速增长的现代化社会，土地植被退化、城市变迁、自然灾害等频繁发生，极大地影响了社会经济和生态环境，对这些国土资源的变化情况进行实时有效的动态监测是极其重要的。近年来，为了满足人们对更高分辨率遥感影像的需求，卫星、传感器等技术不断革新，高分遥感数据的获得变得更加方便高效。因其空间分辨率高，包含的地物细节特征更加详尽和丰富，采用高分遥感探测技术动态监测地球表面成为获得地物变化信息的有效途径[1]。高分遥感影像变化检测是对同一区域不同时刻的影像进行定量分析从而获取地表变化信息的过程，是获取地表变化信息的关键技术之一，在环境变化动态监测、自然灾害评估、城市规划发展等领域有着重要的经济及应用价值[2-4]。

科研人员提出了多种无监督分类的遥感影像变化检测方法，这些方法一般包括两个部分：差异影像的生成与分析。差异影像的生成方法有差值法[5]、比值法[6]、变化向量分析法(change vector analysis，CVA)[7]等。常用的差异影像分析方法有大律法(OTSU)[8]、K-means聚类法[9]、模糊C均值法(fuzzy c-means clustering algorithm，FCM)[10]、期望最大化算法(expectation-maximization，EM)[11]、水平集方法[12]、主动轮廓模型[13]等，这些传统算法大都只处理影像的像元灰度值信息，邻域像元之间的空间关系被忽视，导致图像细节保留能力差及抗斑点噪声能力不足从而不能获得更高精度的变化检测结果。随机场建模是变化检测领域中常用的表达影像空间信息的方法[14]。Bruzzone等[15]首次对遥感影像的空间信息建立马尔科夫(markov random field，MRF)模型进行变化检测。Bruzzone 等[16]提出了一种针对多光谱图像利用HMRF模型和EM方法进行参数估计的变化检测方法。Kasetkasem等[17]基于MRF模型，提出了利用可以描述像素与其邻域像素统计相关的特性结合最大后验概率(maximum a posterior，MAP)准则的变化检测算法。Tso等[18]提出将模糊融合信息与空间上下文相结合建立MRF的分类方法。Wang等[19]提出基于非参数密度估计和MRF理论的变化检测方法。Wei等[20]提出了一种基于多阈值集成的MRF无监督分类遥感影像变化检测算法。Zhang等[21]在建模过程中提出了设置相对同质性指数来改进空间权重的变化检测方法。Liu等[22]提出利用非局部均值滤波器思想结合MRF重新定义像素间邻域关系进行遥感图像变化检测。Gong等[23]提出新的MRF能量函数并加入多项式修改模糊聚类得到隶属度信息用于SAR图像的变化检测。上述研究表明，传统的马尔科夫随机场虽然考虑了像元间的空间关系，但是对邻域像素的空间关系定义不够准确，无法避免变化检测结果中出现的变化信息缺失及边缘分割清晰度不高的问题。

针对上述问题，本文提出了一种顾及像元空间引力的变权重MRF高分遥感影像变化检测方法。首先，对由变化向量分析法得到的差异影像进行模糊C均值聚类，得到隶属度信息；然后，基于高斯混合模型对差异影像进行特征场建模，根据空间引力模型改进Potts模型并将聚类得到的像元隶属度信息引入到该模型中，对邻域像素间的空间约束关系进行重新定义，增强图像细节保留能力，建立空间标记场模型；最后，设置参数公式，自适应改变组合能量函数中空间标记场能量和光谱特征场能量的权重，提高边缘分割的精度，采用ICM优化迭代，得到最终的分割结果。

1 改进的变化检测算法原理

假定X1和X2分别为t1、t2两个时刻获取的同一区域大小均为M像素×N像素、具有相同空间分辨率的两景高分辨率遥感影像，数据经过严格几何配准。本文方法的流程图如图1所示。

图1 变化检测流程图

1.1 差异影像生成

变化向量分析法是对两时相光谱矢量变化的大小和方向进行描述来检测变化信息的一种方法。

1.2 模糊C均值聚类算法

模糊C均值聚类是一种灵活的模糊分类，在隶属度的概念中加入模糊思想，通过对样本建立类别的不确定描述，能够更客观地反映实际情况，在聚类分析中应用更加广泛。

假设差异影像X={x1，x2，…，xN}是由N个向量构成的数据集，c为模糊类别数。FCM聚类算法的目的是通过最小化目标函数J来获取差异影像中数据xi(1≤i≤n)对于第k类的隶属度uik(1≤k≤c)，从而实现最优聚类。

1.3 变权重MRF建模

马尔科夫随机场可以表达图像像元之间的空间相关信息，因此在图像分析和遥感影像变化检测等领域中应用广泛。MRF模型对图像中像素与其对应的分类标记的联合条件概率进行建模，将图像分割问题转化为标记场X的概率估计问题，在最大后验概率(maximum posterior probability，MAP)准则下，从场能量的角度，使每个像素的光谱特征场能量与标记场能量之和最小，以获得最优结果。

假设由1.1节变化向量分析法生成的差异影像为X={x1，x2，…，xN}，L={l1，l2，…,lc}是差异影像中像素的分类标记，c为分类数，在最大后验概率准则下，像素的分类标记表示如式(1)所示。

L=argmax{P(L)P(X|L)}

(1)

式中：P(L)是差分影像类标签的先验概率分布；P(X|L)是差分影像中像元灰度值的联合概率密度函数。在MRF变化检测方法中，求取最大后验概率就相当于计算下列能量函数的最小值，即将概率最大问题转化为求取能量最小的优化问题，能量函数如式(2)所示。

UMRF(xi)=Uspectral(xi)+Uspatial(xi)

(2)

式中：Uspectral(xi)为差异影像的光谱能量函数；Uspatial(xi)为通过像素xi的邻域像素Ni来计算的局部空间能量函数。在常量权重作用下，空间标记场分量较大时，无法体现数据的真实分布，降低图像分割结果的精度。当光谱特征场能量占据主导地位时，分割过程中对像素和邻域像素间的空间相互作用关系利用不足，减弱了图像分割结果的边缘一致性，即使人工取得一个合适权重来使这两个能量达到平衡状态，最后估计的结果也不具有全局性。因此，本文提出自适应可变权重的思想，设置一个根据迭代次数变化的权值函数，使总能量函数趋于平衡。在总能量函数中引入变权重值α(t)后，图像的能量函数可以表示为式(3)。

UMRF=α(t)·Uspectral+Uspatial

(3)

在自适应权值函数的选取过程中，设置权值函数为关于迭代次数t的指数函数，在精度上对线性权重和指数型权重进行比较。分别取c为0.01、0.04、0.06、0.08、0.1、0.5、1、5、10，在高分遥感影像上进行的大量实验证明了虽然指数型权重函数运行时间较慢，但与线性权重相比，指数型权重的精度更高。本文对文献[24]中的权值函数做出了参数改进，将权值函数表示如式(4)所示。

α(t)=c1*γt+c2

(4)

式中：c1、γ、c2是常数，实验中c1取0.01，γ取0.95，c2的值设置为1/4；t为迭代次数。光谱特征场能量在前期的迭代过程中会占据较大比重，从而学习到参数的全局值，之后α接近于c2，空间标记场能量占据较大比重，分割结果的区域一致性更好。像素xi的光谱特征场能量函数如式(5)所示。

(5)

(6)

(7)

势函数β大于0且是固定参数，本文中势函数的值设为0.8。β为基团参数，是一个常量。当相邻像素强度值相差越大，处于该点的势函数的值就越接近于0，两个像素被分入一类的可能性就越小，当两个像素强度值相等时，势函数的值达到最小，即-β。Ni为像素xi(i∉Ni)的空间邻域像素。l(xi)和l(xj)(j∉Nj)分别为像素xi的分类标签和相邻标签。

1.4 空间引力模型

选用Potts模型在马尔科夫随机场中进行空间标记场建模，式(7)即为Potts模型。可以看出，传统的Potts模型是按中心像素的类别与其邻域像素的分类类别是否相等定义的，相等为“1”，不相等为“0”，这种直接定义像素间的空间邻域关系的方式过于绝对化，容易造成空间信息的过度利用。为了解决上述问题，本文提出利用空间引力模型将模糊C均值得到的像素隶属度信息引入到改进的Potts模型中，更好地定义像素间的空间相关性。Potts模型可以改写为式(8)。

(8)

式中：wij是中心像素xi与其邻域像素xj之间的空间引力。空间引力公式如式(9)所示。

(9)

式中：G为引力常数，用来表示调节空间约束对聚类目标函数的贡献(通过大量实验，本文对引力常数作出改进，设G=0.6时，效果最好)；uki表示邻域窗口中心像元xi属于第i个类别的隶属度；ukj表示邻域窗口内的第j个像元xj属于第k个类别的隶属度(各像素隶属度信息可由模糊C均值聚类算法计算得出)；Rij表示中心像素xi与其相邻像素的欧氏空间距离。

由于空间引力模型引入了聚类隶属度概率和中心像素与相邻像素的距离，在一定程度上解决了细节信息变化导致结果过于平滑的问题。

1.5 算法流程

上述方法具体实现步骤如下。

步骤1：采用变化矢量分析法(CVA)对影像X1、X2求取差异影像D。

步骤2：通过模糊C均值聚类算法对差值影像D进行初始聚类分割，得到隶属度信息，利用空间引力模型将隶属度及像素的空间距离信息引入到Potts模型中建立空间标记场模型。

步骤3：由高斯混合模型以差异影像D的像素为基础建立光谱特征场模型。

步骤4：在贝叶斯框架下，设置组合能量函数的权重公式，采用条件迭代模型(ICM)不断迭代更新，取得优化的MAP解，获得最终的高分遥感影像变化检测结果。

2 实验与分析

2.1 算法运行环境

将IR-MAD、K-means、FCM、OTSU、MRF五种方法与本文算法进行实验和对比评估。所有算法均在MATLAB 2014b上编程实现。运行环境为Intel(R)Core(TM)i5-4200，2.30 GHz主频，处理器4GBRAM。

2.2 实验数据

实验一中所用两时相遥感影像是由QuickBird传感器分别于2004年5月和2009年5月获取的重庆市某地区影像。研究区域影像大小为1 016像素×672像素。QuickBird影像原始数据及人工绘制的参考图如图2所示。

实验二中所用两时相遥感影像是由Ikonos传感器分别于2009年7月和2013年7月获取的重庆市某地区遥感影像。研究区域影像大小为1 272像素×808像素。Ikonos影像原始数据及人工绘制的参考图如图3所示。

图2 实验一数据

图3 实验二数据

2.3 实验结果与分析

表2、表3分别显示了两组实验中IR-MAD、K-means、FCM、OTSU、MRF及所提出算法的定量评估结果。图4、图5为实验结果。IR-MAD算法简单但是因其没有顾及邻近像元之间的空间关系，结果图像中出现椒盐效应，且漏检率较高；K-means算法与FCM算法生成的变化检测结果与真实参考结果接近，K-means算法属于硬分割，对噪声非常敏感，一般来说噪声点与周边区域存在着较大差异，即使噪声点在某个类的内部也不会被分为该类，这两种算法的不足之处都在于没有考虑差异影像中像元之间的空间关系，结果图中含有大量噪声，影响了变化检测的效果；OTSU法漏检的区域较多，漏检率较高，且对噪声的抑制能力较差，从而影响了变化检测的效果；传统 MRF 算法虽然考虑了邻域像素间的空间相关性，但由于权重参数是固定的，造成了过度检测的现象，未变化像素检测为变化像素的数目较多，其检测结果过于平滑；本文算法由于空间引力模型的改进及变权重思想的引入，更加合理地利用了像素间的空间关系，较传统 MRF 方法减少了实际未变化被检测为变化像素的数量，提高了抗噪能力，增强了区域一致性，提高了变化检测的精度。

表1 实验一变化检测结果精度评定

表2 实验二变化检测结果精度评定

图4 实验一变化检测结果

图5 实验二变化检测结果

3 结束语

本文针对传统遥感影像变化检测方法对于空间信息利用不合理的问题，利用改进的空间引力模型对像元间的空间关系进行重新表征，增强了图像细节的保留能力。另外，对于传统MRF变化检测方法权重固定导致边缘分割存在缺陷的问题，设置自适应权重函数，增强了区域一致性。经过对比实验表明，本文方法较经典变化检测算法及传统MRF方法边缘分割能力更强，具有更高的变化检测精度。