结合邻域信息和结构特征的遥感影像变化检测

叶沅鑫，孙苗苗，王蒙蒙，谭 鑫

1. 西南交通大学地球科学与环境工程学院，四川 成都611756； 2. 高速铁路安全运营空间信息技术国家地方联合工程实验室, 四川 成都611756

随着我国城市化与工业化进程不断加快，土地作为社会经济发展的重要支撑，为人类活动和社会发展提供了基础的生产资料。准确快速地获取土地利用变化信息，对社会发展、环境保护、自然资源管理等具有重要意义[1-3]。而常规的地面调查手段不仅耗时长、成本高，且部分区域难以到达，无法人为判断。近年来，遥感成像技术的快速发展为人类获取地面信息提供了新的手段，如利用遥感影像对土地利用变化信息进行检测，此方法已成为变化检测领域的研究热点。

目前，遥感影像变化检测方法，主要可以分为面向对象法和基于像元法。面向对象法[4-6]是将空间邻近、光谱相似的对象作为基本处理单元，该方法虽然能够利用高分辨率影像丰富的空间信息，但在生成对象时的最优分割尺度却难以确定[7]。基于像元法以单个像元为处理单元，主要包括代数法[8-9]、基于深度学习的方法[10-13]和结合多特征的方法。其中，代数法虽简单易行但检测结果椒盐现象严重。基于深度学习的方法虽然能够通过深层网络结构自动、多层次地提取复杂地物的抽象特征，但由于目前公开的变化检测数据集较少，难以获取大量的训练样本，在现阶段具有一定的局限性[7]。结合多特征的方法[14]可以综合利用影像的光谱、纹理等多种信息，便于构建稳健和适用的变化检测模型，因此得到了学者们的广泛关注。

在结合多特征的变化检测方法中，邻域信息是最常用的特征之一，可以将其与灰度信息、空间信息等影像特征进行结合，优化检测结果[15]，还可以与纹理信息结合，提高检测精度[16]。其中邻域相关影像(neighborhood correlation image，NCI)[15]是一种能够反映影像局部邻域相关性的上下文信息，对影像是否发生变化十分敏锐。但NCI是通过固定窗口分析获得的，因此窗口大小会直接影响检测结果，窗口较大时会削弱影像细节处的变化，而窗口较小会产生严重的椒盐现象。文献[17—19]提出了基于自适应延伸技术的变化检测方法，该类方法依据同类地物像元间的光谱相似性获得以像元为中心的自适应区域，然后基于自适应区域提取变化信息。虽然能够有效削弱检测结果中的椒盐现象，但由于该类方法通常将自适应区域的平均灰度值作为中心像元值，因此会在一定程度上影响变化地物边界的准确性。故仅利用单一的空间邻域信息难以获得理想的检测结果。基于此，本文引入一种反映邻域相关性的测度(匹配误差)作为NCI的补充，其基本思想是：在完全配准的双时相影像中，未变化区域对应的局部邻域具有很高的相关性，且相关性最强的邻域位置几乎没有偏移，或者偏移很小，即匹配误差很小；而对于变化区域，由于地物类型发生了变化，双时相影像中相关性最强的邻域会远离对应邻域位置，即匹配误差较大。

此外，考虑到多时相遥感影像中光谱差异通常较大，而几何结构特征相对于灰度信息更加稳定，受光谱差异影响更小。故本文将几何结构特征属性融入变化检测方法中，提出了一种邻域信息与结构特征相结合的遥感影像变化检测方法。该方法通过将NCI和匹配误差相结合获得更加稳健的邻域信息，利用方向梯度通道信息[20]提取影像结构特征，然后将上述3种属性特征作为决策树的分类属性，进行二值分类，获得初始变化检测结果，最后借鉴马尔可夫随机场(Markov random field，MRF)模型在表征图像像素空间关系方面的成功经验[21]，对初始检测结果进行MRF优化，得到具有空间一致性的变化检测结果。

1 本文方法

本文方法主要包括以下3个部分：①获取邻域信息：通过提取双时相影像间的NCI特征，获得局部邻域的上下文信息，同时在影像间利用邻域像素的互相关性进行模板匹配获得匹配误差。②结构特征提取：利用方向梯度信息构建结构特征描述符，提取双时相影像的结构特征。③生成检测结果：将NCI、匹配误差、双时相影像结构特征作为属性信息进行决策树分类，获得初始变化检测结果，然后以双时相影像波段差异图为特征场，初始变化检测结果为标记场进行MRF优化，获得变化检测二值图。方法的总体流程如图1所示。

1.1 获取邻域信息

1.1.1 邻域相关影像

邻域相关影像是利用双时相影像间局部邻域的光谱属性得到3个表示上下文信息的特征影像，即相关系数影像，斜率影像，截距影像。它们用于变化检测的基本原理为：若双时相影像对应邻域范围内的地物没有发生变化，则该区域的光谱值倾向于线性相关，即相关系数较大、趋近于1，斜率趋近于1，截距趋近于0；若对应邻域范围内的地物发生了变化，则该区域的光谱值倾向于线性无关，相关系数较小、偏离1，斜率偏离1，截距偏离0。为获取邻域相关影像，在双时相影像对应像元位置开辟一个3×3的动态窗口，然后采用式(1)—式(4)分别计算对应窗口邻域的相关系数r，斜率a，截距b，并将值赋给窗口中心像元位置，构建相关系数影像，斜率影像，截距影像，统称NCI[15]。

图1 本文方法总体流程Fig.1 Flowchart of proposed method

(1)

(2)

(3)

(4)

式中，cov12为双时相影像对应邻域窗口内所有波段亮度值的协方差；BVi1、BVi2分别为前、后时相对应邻域窗口内所有波段第i个像元的亮度值；μ1、μ2分别为前、后时相对应邻域窗口内所有波段亮度值的平均值；n为邻域窗口内所有波段的像元总数；s1、s2分别为前、后时相对应邻域窗口内所有波段亮度值的标准差。

1.1.2 匹配误差

虽然NCI能度量影像间的变化信息，但检测结果中易存在大量椒盐。为充分利用影像的空间信息，本文在计算双时相影像的相关性时，不仅计算对应邻域位置(记为中心位置)的相关系数，还通过模板匹配的方式，利用相关系数计算局部邻域的匹配误差，以此来获取更多的变化信息。其中，匹配误差用于变化检测的基本原理为：对于已经配准过的双时相影像，若对应邻域未发生变化，则相关性最强的位置几乎没有偏移，匹配误差较小。若对应邻域发生了变化，相关性最强的位置会存在较大的偏移，匹配误差较大。

为获得匹配误差，将前时相影像作为基准影像，在后时相影像上对应像元位置开辟一个9×9大小的搜索窗口，然后在搜索窗口内进行模板滑动，计算模板与前时相邻域的相关系数，获得相关系数矩阵，最后通过计算相关系数最大的点到矩阵中心的欧氏距离d得到匹配误差(匹配误差=d)。图2以黄色点(黄色十字标识符标识的点)和粉色点(粉色十字标识符标识的点)为例对未变化点和变化点的匹配误差进行可视化展示。其中，图2(a)和图2(b)分别展示了两点在前、后时相影像中的位置，图2(c)、图2(d)分别为黄色点和粉色点的相关系数矩阵可视化柱状图，其中采用红色十字标识符标识相关系数矩阵中心位置，红色柱子标识相关系数最大的点所在的位置。由图2可以看出，黄色点为未变化点，矩阵中心的相关系数较小，难以给出正确的变化信息，但该点的匹配误差较小(匹配误差为1)，可以给出正确的变化信息。粉色点为变化点，矩阵中心的相关系数较大，难以给出正确的变化信息，但该点的匹配误差较大(匹配误差约为5.66)，可以给出正确的变化信息。以上初步说明了，相关系数的大小难以精确地反映影像的变化信息，将匹配误差引入到变化检测中，可作为NCI的有效补充。

1.2 结构特征提取

考虑到通常情况下，双时相影像间光谱和灰度信息差异较大，本文基于方向梯度通道特征(channel features of orientated gradients,CFOG)[20]构建结构特征，以抵抗影像间的辐射差异和灰度差异对检测结果的影响。整个构建过程主要包括以下3个步骤：①利用多方向梯度信息构建方向梯度通道；②对方向梯度通道进行特征通道卷积形成CFOG；③对CFOG进行通道叠加形成结构特征。其中，方向梯度通道由m层构成，每层的方向梯度通道gθ(θ为该层通道所在的划分方向)在点(h,l)处的梯度幅值gθ(h,l)是通过点(h,l)处的水平方向梯度幅值gx(h,l)和垂直方向的梯度幅值gy(h,l)根据式(5)[20]计算得到

gθ(h,l)=abs(cosθ×gx(h,l)+sinθ×gy(h,l))

(5)

在形成方向梯度通道之后，在X和Y方向上采用二维高斯核、梯度方向上采用[1,2,1]T核进行特征通道卷积和归一化处理形成CFOG，最后将CFOG所有通道对应位置求和取平均构建单通道结构特征，以减少特征维数。图3显示了双时相影像结构特征的构建过程，从中可以清晰地看出，由于CFOG不受大气折光和灰度差异的影响，因此基于CFOG的结构特征能够较好地反映地物的几何结构信息。

1.3 生成检测结果

1.3.1 决策树分类

在结合邻域信息和结构特征的多特征变化检测问题中，涉及多个属性特征。决策树作为一种非参数化以及能够在样本学习中自动估计属性信息重要程度的分类器，能较好地适用于结合多特征的变化检测方法，故本文采用决策树作为分类器获得初始变化检测结果。在进行决策树分类时，首先依据参考变化图中变化像元和未变化像元的比例选取训练样本。然后以训练样本中像元的属性特征(即NCI、匹配误差、前时相结构特征、后时相结构特征)作为输入，其对应的检测结果(变化像元记为255，未变化像元记为0)作为输出，采用自顶向下的递归方法，以基尼系数为标准确定树节点的最优属性及最佳分割点，并在结点上进行属性值的比较，依据像元对应的不同属性值判断从该结点向下的分支，直到某个属性的子集中所有样本都属于同一类别，完成决策树训练。最后依据训练好的决策树分类规则对全部像元进行判断，获得初始变化检测结果。

1.3.2 MRF优化

考虑到决策树分类方法是基于像元的分类方法，其检测结果过于破碎，存在一定的椒盐噪声现象。MRF作为概率论的一个分支理论，可以很好地刻画影像中邻域像元属性间的依赖关系，具有较强的抗噪性，其在图像分割[22]、变化检测领域取得了广泛应用[23-24]，故可将MRF模型用于优化决策树检测结果。即以双时相影像波段差异图为观测场、决策树检测结果为初始标记场，依据MRF与Gibbs随机场的等价性将标记场的全局最优估计问题转化为能量之和的最小化问题，并通过迭代获得最终二值变化图。

具体实现思路为：①参数设置，即初始化标记场，初始化迭代次数为1，总迭代次数mtotal=5(考虑到基于MRF模型的变化检测方法中存在过度平滑的问题[25]，本文采用较小的迭代次数对初始检测结果进行优化)；②依据当前标记场分类情况，计算特征场参数和标准差；③计算特征场能量和标记场能量，并依据能量最小化原则更新当前分割结果；④重复步骤②③④直到迭代轮数达到mtotal，获得最终二值变化图。

2 试验与分析

2.1 试验数据

本文通过选取两组不同传感器的数据集，检验算法的有效性，其中每组试验区域均包含两幅已配准的不同时相的遥感影像和一幅参考变化图，在参考变化图中白色部分表示变化区域，黑色部分表示未变化区域。

第1组数据集为墨西哥数据集(图4)，其中，图4(a)、(b)分别为在2000年4月、2002年5月获得的墨西哥郊外的两幅Landsat 7 ETM+4遥感影像，影像大小均为512×512像素，分辨率为30 m。前时相影像显示了火灾尚未发生时的情形，后时相影像可以清楚地看出火灾发生的位置及范围(图4(b)中颜色较暗的区域)；参考变化图通过目视解译获得。

第2组数据集为印度尼西亚数据集(图5)，该数据集是由Quick Bird卫星于2004年4月和2005年1月拍摄，影像大小为1500×1500像素，分辨率为2.5 m。该数据集反映的是印度尼西亚地区受海啸侵袭影响的地表变化情况，此时间段正直海啸过后，植被、水体等土地覆盖类型发生显著变化，且两幅影像之间存在大气折光差异(前时相影像中存在大气折光，如图5(a)中矩形框所示。参考变化图通过目视解译获得。

注：黄色点为未变化点，粉色点为变化点图2 匹配误差示意Fig.2 Schematic diagram of matching errors

图3 基于CFOG的结构特征构建过程Fig.3 Construction process of structural features based on CFOG

图4 墨西哥数据集Fig.4 Mexico data set

图5 印度尼西亚数据集Fig.5 Indonesia data set

2.2 试验设计

为了验证本文方法的优越性，将其与4种基于像元的变化检测方法进行比较，并在表1中对本文方法和其他4种方法所使用的属性特征进行了说明。其中，方法I为基于灰度信息的变化向量分析(change vector analysis，CVA)法；方法Ⅱ为将影像灰度信息和NCI作为分类属性的NCI法[15]；方法Ⅲ为基于自适应上下文信息(adaptive contextual information)的方法[19]，记为ACI法；方法Ⅳ为CVA与局部自相似纹理差分度量(local-similarity-based texture difference measure，LSTDM)相结合的方法[16]，记为CVA+LSTDM法；方法Ⅴ为结合邻域信息(NCI、匹配误差)和结构特征的本文方法。

2.3 试验结果与分析

采用5种变化检测方法分别对两组数据集进行检测，检测结果如图6、图7所示，其中所有的图(a)为参考变化图，图(b)为CVA法检测结果，图(c)为NCI法检测结果，图(d)为ACI法检测结果，图(e)为CVA+LSTDM法检测结果，图(f)为本文方法检测结果。其中白色表示变化，黑色表示未变化。检测精度分别见表2、表3。

表1 本文对比的变化检测方法

由图6(b)、图7(b)可以看出，CVA法可以识别墨西哥数据集中的主要变化区域，但椒盐现象明显，虚警率较高。对于印度尼西亚数据集，CVA法无法正确识别植被区域内像元间的差异，虚警率也较高。这是因为CVA利用的是像元间的灰度差来探测变化像元，对光谱差异敏感。

由图6(c)可以看出，NCI对存在光谱差异的地区灵敏度较高，但存在大量虚检，无法正确识别变化区域内部的小块未变化区域(如图6(c)中灰色框所示)。由图7(c)可以看出，NCI可以通过局部相关性分析，较好地克服单一像元的孤立性，正确识别植被区域内单个像元间的光谱差异，但难以正确识别存在大气折光和水质差异的未变化区域(如图7(c)中灰色箭头、灰色框所示)，虚警率仍较高。

图6 墨西哥数据集检测结果Fig.6 The detection results of Mexico data set

图7 印度尼西亚数据集检测结果Fig.7 The detection results of Indonesia data set

表2 墨西哥数据集5种方法检测精度

表3 印度尼西亚数据集5种方法检测精度

从ACI法对两个数据集的检测结果(图6(d)和图7(d))可以看出，基于自适应延伸技术的邻域信息可以有效减弱椒盐现象，削弱大气折光和水质差异对印度尼西亚数据集检测结果的影响(如图7(d)中灰色箭头、灰色框所示)，总分类精度和Kappa系数较高、虚警率较低，但由于ACI法将延伸区域的灰度均值作为中心像元值，会影响变化区域边界的准确性(如图6(d)中灰色框所示)，且存在漏检现象(如图6(d)和图7(d)中灰色椭圆框所示)。

图6(e)、图7(e)检测结果表明，CVA和LSTDM相结合的方法受大气折光、水质差异影响较小、椒盐现象较弱，但该方法对变化信息的灵敏性较弱，无法正确识别光谱差异较弱的变化区域(如图6(e)和图7(e)中灰色椭圆框所示)，漏检率较高。

图6(f)、图7(f)的检测结果表明，本文方法受大气折光和水质差异影响较小，椒盐现象不明显，相比于以上4种方法，具有更高的总分类精度和Kappa系数。这是因为匹配误差通过刻画局部邻域对应像素间的偏移距离，提供了一种反映邻域相关性的测度，对NCI进行了有效补充，并通过引入结构特征，增强算法对影像间光谱差异的稳健性，获得较好的检测结果。

综合上述试验结果可以看出，结合邻域信息和结构特征的本文方法可以有效检测出地物的变化情况，降低检测结果中的椒盐现象。获得5种方法中总分类精度最高、Kappa系数最高，漏检率与虚警率都较低的检测结果。

3 结 语

针对利用传统邻域信息进行变化检测时存在椒盐现象严重，漏检率和虚警率较高等问题，本文提出了一种结合邻域信息和结构特征的变化检测方法。该方法构建了一种新的反映邻域相关性的度量测度(即匹配误差)，并将其与NCI结合来获取更加稳健的邻域信息，与此同时，引入结构特征增强算法对影像间光谱差异的稳健性，然后通过决策树分类和MRF优化获得具有空间一致性的检测结果。通过采用两组试验数据，与CVA和NCI等4种方法进行对比，证明了本文方法能够获得精度更高、虚警率和漏检率均较低的变化检测结果。由于在本文方法中匹配误差的搜索窗口是一个重要的参数，针对不同复杂度的双时相影像需要对搜索窗口的大小进行调整，因此搜索窗口大小的自适应性是未来的研究方向之一。