结合相关系数和特征分析的植被区域自动变化检测研发

潘建平， 徐永杰， 李明明， 胡勇， 王春晓

(1.重庆交通大学土木工程学院，重庆 400074； 2.重庆市规划和自然资源调查监测院，重庆 401123； 3.自然资源部海南基础地理信息中心，海口 570203)

0 引言

遥感变化检测是在不同时间观察判断一个地物对象的状态有无变化、发生怎样变化的过程[1-2]。通过变化检测可以有效地掌握地表信息的变化，从而更好地管理自然资源，推动社会发展。其中耕地、林地和草地等植被类型作为地表覆盖类型的重要组成部分，影响着人类和自然的方方面面，植被覆盖区域变化检测对于保护生态红线和守住耕地红线具有重大意义[3]。

传统的遥感变化检测方法基本可以分为2大类： 分类后比较法和直接比较法[4]。分类后比较法首先对前后期遥感影像进行分类，然后根据分类的结果判断是否发生变化以及发生何种类型的变化，如梅树红等[5]为检测林地变更，在构建归一化植被指数(normalized difference vegetation index，NDVI)和提取纹理特征的基础上进行决策树分类来获得变化区域，该方法极易受分类误差的影响。直接比较法通常利用像元之间的光谱差异构建前后期的差异图像，然后在差异图像上提取变化像元和未变化像元，其原理简单，能够消除误差累积的影响，因此受到很多专家学者的广泛使用。Nielsen等[6-7]先后提出多元变化检测(multivariate alteration detection,MAD)和迭代加权多元变化检测(iteratively reweighted MAD,IR-MAD)方法进行多元遥感变化检测研究； 王晓东等[4]选用交叉相关系数来构建两时相的变化强度图像，通过马尔可夫随机场-最大后验估计(maxium a posteriori estimation of Markov random field，MRF- MAP)方法逐像素的提取变化区域； 黄恺等[8]构建时空自相关指数对植被进行变化检测，并考虑了像元的邻近信息。但以像元作为最基本的检测单元，往往会忽略了高空间分辨率遥感影像丰富的上下文信息，因此针对高分辨率遥感影像，国内外学者多展开面向对象的变化检测研究[9-11]。Walter[12]和Zomeni等[13]利用地理信息系统(geographic information system,GIS)矢量数据获取影像像斑进行变化检测； 王琰等[14]在利用矢量数据获取像斑的基础上进行二次划分，获得了更为精细的子像斑，然后以子像斑为检测单元进行相关性判断； 佃袁勇等[15]通过多尺度分割获取地物对象并进行变化检测。

在分割的基础上提取地物特征进行比较是面向对象的直接比较法的基本检测思路[16]。选取各个特征值组成特征向量，然后计算特征向量的变化强度或相似性程度，比较前后期同一区域的变化强度大小或特征向量的相似性程度来筛选变化的结果。然而，一般的光谱和纹理等特征都具有通用性，针对特定的变化检测目标容易产生伪变化区域。常用的解决办法是先从影像上识别出目标对象，然后针对目标对象展开变化检测[17]； 或选取多种特征进行筛选和降维处理，提取出能有效识别目标地物的特征[18]。这无疑给检测工作带来了较大的工作量，并且对目标识别算法有较高的要求。因此，针对特定目标的变化检测工作，如何提高变化检测精度、减少伪变化信息的产生，是面向对象的直接比较法需要重点解决的问题。

针对高分辨率遥感影像中植被区域的变化检测，本文采用面向对象的变化检测方法。首先对前后期高分影像进行分割得到地物图斑，并提取光谱和纹理特征构建特征向量，计算相关系数值判断向量之间的相关性，然后通过特征分析选取合适的伪变化去除特征，利用伪变化去除算法剔除虚检的变化，从而提高变化检测的精度。同时，为了提高变化检测的效率和自动化程度，文章设计研发了一款图斑级变化检测工具软件。

1 研究思路与关键技术

文章进行变化检测使用的数据包括： 前后2期高分辨率遥感影像和前期地理国情普查成果数据。主要的研究思路分为影像分割、特征提取和变化检测3个步骤[19]，其中变化检测又分为变化信息提取和伪变化去除。图1为变化检测主要技术路线。首先，在后期遥感影像上叠加前期普查成果数据，对其进行矢量约束的多尺度分割，然后将后期分割得到的结果叠加在前期影像上进行棋盘分割，从而得到前后期对应的分割图斑； 紧接着，提取前后期每个分割图斑的光谱和纹理特征； 然后，根据提取的地物特征构建特征向量，计算特征向量之间的相关系数来提取变化的信息； 最后，对相关系数法提取的变化结果进行伪变化信息去除，得到最终的变化检测结果。

图1 技术路线

1.1 影像分割

影像分割是将一幅完整的遥感影像分成若干个相互独立的子区域，是面向对象变化检测方法的首要步骤，也是关键步骤。本文将后时期的遥感影像与前时期的普查成果数据进行叠加套合，利用eCognition软件进行多尺度分割； 将分割的结果导出并叠加在前时期的遥感影像上进行棋盘分割，得到前后期完全一致的分割图斑。

多尺度分割是eCognition软件中最常用的一种分割算法，是根据影像的光谱和形状特征进行分割，并采用自下而上的思路合并相邻的同质区域。软件中通过设置紧致度(compactness)、形状因子(shape)和尺度(scale)参数来控制多尺度分割的结果，针对高分辨率遥感影像，通常形状因子介于[0.1，0.5]之间，紧致度介于[0.5，0.99]之间，分割尺度则依据变化检测要求的最小检测面积来确定[1]。在多尺度分割的过程中添加前期成果数据的约束，可以保证前后未发生变化的区域具有相同的边界[20]，而叠加在后期影像上进行分割，可以提取出后期发生变化的区域边界[21]。棋盘分割是将一个父对象分割为许多正方形的子对象，在eCognition软件中通过设定对象尺寸(Object Size)来决定子对象的大小，如果设定的Object Size值大于父对象的尺寸，则可以得到前后期一致的分割图斑。

1.2 相关系数

相似性度量是一种代价函数，可以反映特征之间的相似性。通过提取前后期影像的光谱、纹理特征构建特征向量，带入代价函数进行相似性分析，最后通过确定变化阈值从而判定影像的变化情况。常用的相似性度量有欧氏距离[22]、KL距离[23]和向量相似性[9]等。

相关系数由统计学家皮尔逊提出，可用于研究变量之间的线性相关程度，常用来衡量2个样本之间相似性[24]。皮尔逊相关系数r的计算公式为：

(1)

在前期图像分割和特征提取的基础上，利用提取出的光谱和纹理特征构建特征向量V，即

(2)

以相关系数作为描述特征向量之间相似性的度量方式。通过计算前后期对应图斑的相关系数并设置变化阈值便可得到相应的变化图斑，相关系数ρ的计算公式为：

(3)

式中：V1i和V2i分别为前、后2期对应图斑的特征向量；i为图斑序号；n为特征向量的维数即特征的个数。

1.3 植被区域伪变化去除

传统的特征向量相关系数法根据选择的特征判断变化区域的相似性，然而由于不同遥感传感器之间的成像差异、不同季节植被覆盖度不同、光照和大气条件等因素的影响，同一块植被区域在不同时期的光谱和纹理存在一定的差异，导致其被认定为发生了变化。如图2所示为部分农田和林地因为季节和成像差异出现的误检图斑。因此，单纯的利用相关系数进行变化检测会产生一些伪变化信息。

(a) 前期无植被耕地(b) 后期有植被耕地(c) 前期无植被林地(d) 后期有植被林地

图2 伪变化区域

为消除这种伪变化信息的影响，本文通过在前后期影像上选取多组地物样本，统计分析前后期影像上不同类型地物在红光、绿光、蓝光3个波段下的光谱均值，并引入标准差来确定每个波段下各类地物像元值的取值范围。表1为4种典型地物类型的像元统计特征。

表1 各地物像元统计特征

从表1中可以看出红光波段下植被的像元值明显偏低，与道路、建筑物和堆掘地等存在较大差异。所以，当前期是植被类型的区域变化为后期是非植被类型时，对应区域的红光波段特征值会发生明显改变； 并且二者灰度值范围少有重叠，可以排除“异物同谱”现象的影响。因此，选用影像的红光波段均值来对变化信息作进一步的筛选，并构建了前后期影像的红光波段均值的比值。如若前后期地物类型发生变化，则比值较大，反之较小。其表达式为：

(4)

式中：RT1和RT2分别为前期和后期影像的红光波段均值；β为界定的阈值，称为伪变化阈值。在前期相关系数法变化检测的基础上，如果变化图斑满足表达式RT2/RT1>β,则为最终的变化图斑。通过在实验区挑选前后期变化图斑样本，发现变化区域前后时相红光波段比值在0.9以上，即β≥0.9。

1.4 阈值确定

目前常用的阈值确定方法有大津法、构建目标函数法和基于信息熵的方法。这些方法大多基于样本表达的信息量，可以保证检测结果整体的正确率，然而在实际生产中更加注重检测结果的漏检率，即首先要确保不能或较少的出现漏检的图斑，其次再优化检测的虚检率。为此，文章采用基于漏检率的阈值确定方法： 首先设置阈值区间，当漏检率低于10%时，此时的变化阈值为相应对象的阈值。

2 工具软件设计

2.1 开发工具

变化检测工具软件需要对矢量和栅格数据进行交互、编辑等操作，所以需要依托GIS平台。ESRI公司的ArcEngine是基于核心组件库ArcObjects搭建的，可以用来开发嵌入式或者独立的GIS程序，并且支持多种开发环境，如C++，Java和.NET[25]，其中基于.NET框架易于设计可视化界面。文章在Visual Studio 2017平台下引入ArcEngine组件库实现工具软件的开发设计。

2.2 功能设计

图斑级变化检测工具软件主要用于辅助内业进行变化检测。将前后期带有特征属性值的矢量图斑数据导入，通过变化检测工具软件提取变化图斑，叠加前后期遥感影像，辅以人工编辑，提高变化检测区域的发现效率和精度。图3为软件功能设计示意图。图4为工具软件的主界面示意图。

图3 功能设计

图4 软件主界面

输入输出功能模块包括数据加载和结果导出。数据加载包括矢量数据和栅格数据的加载，加载的数据会在软件的主界面进行可视化展示； 结果导出是保存得到的变化检测结果，用户可自行选择保存路径和保存名称。

变化检测功能模块包括相关系数法变化检测、阈值确定和变化提取。文章主要通过构建相关系数这种相似性度量方式来进行变化检测，并设定变化阈值进行变化区域的划分； 变化提取是根据计算得到的相似性度量值和设定的阈值进行变化图斑的提取并在图中高亮显示。

后处理功能模块包括人工编辑和字段删除，主要针对变化检测之后的结果进行相关的处理。人工编辑可以针对变化检测的结果进行人为的干预，例如删除或合并图斑； 字段删除是删除数据属性表中的属性字段。

辅助工具包括视图工具和图斑工具。视图工具包括全图显示、移动视图； 图斑工具包括选择要素、清除要素和属性识别。

3 实验与分析

实验选取重庆市璧山区部分区域，前期数据包括2017年GF-2遥感影像和地理国情普查成果数据，后期为2018年BJ-2遥感影像，空间分辨率都为1 m,影像大小为3 761像素×3 085像素，包含3个可见光波段和1个近红外波段，图5分别为前后期原始影像套合矢量图斑的结果图。研究区包含大量耕地、林地等植被类型，并存在较多植被到建设用地的变化。

(a) 前期栅矢套合(b) 后期栅矢套合

3.1 变化检测实验

首先对2幅影像进行精度较高的配准操作，并利用直方图匹配进行相对辐射校正。然后通过套合前期矢量普查成果数据进行分割，获得前后期对应的地物图斑，并提取近红外波段均值(Mean_NIR)、亮度值(Brightness)、对比度(GLCM_Contrast)以及熵(GLCM_Entropy)4个特征构建特征向量。为了提高变化检测的效率，将上述获取到的特征图斑导入变化检测工具软件进行快速变化检测，从而提取出变化的区域。同时，为验证伪变化去除的有效性，并寻求变化阈值和伪变化去除阈值的最佳组合，分别进行多组对照实验。

3.1.1 实验一

第一组实验设定了9组变化阈值，无伪变化去除操作，如图6所示。灰色图斑表示检测出来的图斑，红色图斑表示标准变化图斑。从图中可以看出，随着变化阈值α的逐渐变大，漏检图斑数量逐渐减少，说明变化阈值增大可以抑制漏检图斑的产生。

(a) α=0.950(b) α=0.960(c) α=0.970(d) α=0.980(e) α=0.985

(f) α=0.990(g) α=0.993(h) α=0.996(i) α=0.999

图6 不同变化阈值检测结果

3.1.2 实验二

第二组实验在第一组实验的基础上添加了伪变化去除，并固定了伪变化阈值β为1.0，图7展示了在9组变化阈值下进行伪变化去除后的结果。对比第一组实验结果可以发现，伪变化去除前，随着变化阈值的增大，检测图斑数量越来越多，虚检图斑数量逐渐增多，漏检图斑数量逐渐减少； 伪变化去除后，在相同变化阈值下，检测图斑数量明显下降，虚检图斑大量减少。但在相同的伪变化阈值下，变化阈值越大，虚检图斑越多。

(a) α=0.950,β=1.0(b) α=0.960,β=1.0(c) α=0.970,β=1.0(d) α=0.980,β=1.0(e) α=0.985,β=1.0

(f) α=0.990,β=1.0(g) α=0.993,β=1.0(h) α=0.996,β=1.0(i) α=0.999,β=1.0

图7 伪变化去除检测结果

3.1.3 实验三

第三组实验分别选取变化阈值α为0.993，0.996和0.999，并分别设定伪变化阈值β为0.9，1.0和1.1，图8展示了2种阈值的相互组合变化检测结果。在相同变化阈值的情况下，随着伪变化阈值的增大，变化结果的虚检图斑逐渐减少。

(a) α=0.993,β=0.9(b) α=0.993,β=1.0(c) α=0.993,β=1.1(d) α=0.996,β=0.9(e) α=0.996,β=1.0

(f) α=0.996,β=1.1(g) α=0.999,β=0.9(h) α=0.999,β=1.0(i) α=0.999,β=1.1

图8 阈值组合检测结果

但是，随着伪变化阈值的一直增大，图斑的漏检数也会随之增加。图9为在相同变化阈值下，2种伪变化阈值带来的漏检图斑情况。从图9中可以看出，当β=1.4时裸露出来的标准图斑要多于β=0.9，说明当β=1.4时漏检的图斑要多于β=0.9。

(a) α=0.993,β=0.9(b) α=0.993,β=1.4

图9 漏检图斑对比

3.2 实验分析

通过3组对比实验结果可以得出： ①由第一组实验结果可以看出变化阈值和漏检图斑数呈负相关关系； ②对比第一、二组实验，在相同变化阈值情况下，添加伪变化去除可以有效减少虚检图斑数量，但整体变化阈值和虚检图斑数呈正相关关系； ③由第三组实验可以看出伪变化阈值和虚检图斑数呈负相关关系，但在一定程度上会影响结果的漏检率。

继续统计不同阈值组合下的漏检率和虚检率，当变化阈值在[0.950,0.999]、伪变化阈值在[0.9,1.1]之间时，各项指标的统计结果如图10—12所示。从图10—12中可以看出，变化阈值α在(0.993,0.999]区间内分别可以达到10%左右的漏检率，满足一定的变化检测需求。因此，综合考虑变化结果的漏检率和虚检率，变化阈值和伪变化阈值在一定范围内的取值越大可以获得越好的变化检测结果，一般选取变化阈值α∈(0.993,0.999],伪变化阈值β∈[0.9,1.1]。

图10 β=0.9时阈值组合分析

图11 β=1.0时阈值组合分析

图12 β=1.1时阈值组合分析

3.3 精度评价

为验证该工具软件的有效性，以图斑为基本单位对变化检测的结果图斑进行漏检率和虚检率等指标统计。将人工目视判读选取的真实变化图斑与当α=0.999，β=1.1时进行检测得到的变化图斑进行对比统计，构建用于精度评价的混淆矩阵，统计结果如表2所示。

表2 变化检测结果精度评价

由表2可知，实验区一共包含1 899个图斑，其中实际变化的图斑有306个，该工具软件检测出来的变化图斑有363个，漏检了26个变化的图斑； 83个实际未变化图斑错检成变化图斑。整个变化检测的结果的正确率为94.3%，虚检率为22.9%，其中漏检率可以达到8.5%。

4 结论

针对传统的变化检测方法存在的问题，结合高分辨率遥感影像变化检测的思路，在传统相关系数方法的基础上增加了波段比值去除伪变化，并通过ArcEngine二次开发设计出图斑级变化检测工具软件，相较于传统的变化检测手段，其具有如下优势：

1)研究变化检测算法，消除了大量繁琐的人工目视判读工作，能较好地在植被变化检测实际生产工作中完成任务，提高了工作效率和检测质量。

2)基于ArcEngine研发了快速变化检测的工具软件，实现了半自动化的变化检测，实验精度较高，可以为自然资源调查等部门提供技术支撑。

由于本文构建的伪变化去除方法是针对植被区域变化为非植被区域，因此无法适用于多类变化检测； 并且方法选取的特征数量和质量有限，只选取一个光谱特征指标进行伪变化去除，可能对于老、旧住宅小区等建筑物附近植被的伪变化识别存在一定的局限性。因此，今后需要在提高伪变化识别的广度和深度方面继续进行研究。