基于GDAL的遥感图像变化检测技术

蒋世豪 ，江 洪

1.福州大学 数字中国研究院（福建），福州 350108

2.福州大学 空间数据挖掘与信息共享教育部重点实验室，卫星空间信息技术综合应用国家地方联合工程研究中心，福州 350108

1 引言

图像变化检测是指通过对同一区域不同时间观测的图像进行处理与分析，来确定目标变化情况的过程。遥感图像因其覆盖地域大的特点被广泛应用于地物目标的变化检测[1-4]，在各领域均具有广泛的应用，如环境监测、土地利用、农作物生长状况监测、森林采伐监测、灾情估计等。因此，开展遥感图像变化检测方法的研究具有重要的应用价值[5-8]。

近年来，国内外的学者们针对不同的遥感数据源以及不同的待检测地物特征，提出了许多快捷、有效的变化检测方法。其中，由于较高的计算效率以及较好的适用度，直接比较法是最为常用的一种方法，其原理是：将两幅同位置不同时相的遥感图像进行逐像元点位的分析比对，构造出差异图像后，按照一定的阈值，确定出发生变化的点位[9-11]。但是，构造差异图像及阈值的判定一直是直接比较法的两个难点。构造差异图像需要依据待检测地物的特征，选取敏感波段或进行波段组合来有效反映该种地物的变化状态。利用阈值可以将差异图像的像元分为两类，一类是明显发生变化的区域，另一类是未明显发生变化的区域[6，12]。其难点是如何自动确定阈值，使变化阈值具有更好的普适性。目前，一般采用阈值分割的方法来确定变化阈值，并得到明显发生变化的点位[13-15]。

为了从遥感图像中提取有用信息，传统对遥感图像的实验或单机模式大多是使用专业的图像处理软件（如ENVI、ARCGIS）等。这种方式虽然也能够对遥感图像进行处理和分析，但是更多地是利用软件封装好的方法且离不开人为重复性的输入参数，这样会难以避免过多的人为干预所造成的误差[16-18]。同时，由于遥感图像包含的信息量巨大，这种方式不能有效地扩展遥感图像的处理算法。本文基于GDAL（Geospatial Data Abstraction Library）读取遥感图像的栅格像元数据并提取其对应的波段信息、坐标投影、仿射变换矩阵等元数据，结合遥感图像处理技术和聚类分析技术等，构建了一套遥感图像变化检测的技术流程，实现了遥感图像变化检测的自动化处理，可以为区域尺度的植被长势监测提供决策支持。

2 算法设计

2.1 技术路线

本文算法的技术流程如图1。首先，利用GDAL 读取遥感图像的像元数据及空间元数据；第二步，通过波段组合，计算出能够有效消除地形阴影影响（包括本影与落影）且能够有效反演植被长势信息的植被指数SEVI（Shadow Elimination Vegetation Index）；第三步，对不同时相的SEVI 结果进行标准化，以及利用图像差值法计算出差异图像；第四步，利用OTSU 法自动提取明显发生变化的像元点位，该点位包括植被长势增加的点位及植被长势减少的点位；第五步，利用K均值聚类法提取出变化点位聚集区域的核心，并根据每一核心的所有点位的相对行列号，自动识别出变化区域；第六步，依据前几步获取到的检测结果及第一步提取的空间元数据，输出检测结果。

图1 技术流程图

2.2 GDAL简介

GDAL 是使用C/C++语言编写的一套用于读取空间数据的跨平台开源库，它利用抽象数据模型来表达所支持的各种文件格式，还有一系列命令行工具来进行数据转换与处理。它提供了对多种栅格数据的支持，包括Arc/Info ASCII Grid（asc）、GeoTiff（tiff）、Erdas Imagine-Images（img）、ASCII DEM（dem）等格式[19]。值得一提的是，GDAL 使用抽象数据模型（Abstract Data Model）来解析它所支持的数据格式。抽象数据模型包括数据集（Dataset）、坐标系统、仿射地理坐标转换（Affine Geo-Transform）、大地控制点（GCPs）、元数据（Metadata）、栅格波段（Raster Band）、颜色表（ColorTable）、子数据集域（Subdatasets Domain）、图像结构域（Image_Structure Domain）、XML域（XML Domains）。

基于GDAL 开发遥感图像算法具有明显的优势。这是因为遥感图像是一种包含空间位置特征的地理空间数据，这使它区别于一般的计算机数据。普通的计算机图像处理库不能有效地根据空间数据结构解析出空间元数据，而使用GDAL不仅能够快速读取多种格式遥感图像上的像素值，还能够根据原有数据结构有效地提取遥感图像的投影坐标系、仿射变换等信息，这有助于进行空间分析以及输出包含地理坐标位置的结果图像。基于GDAL的诸多功能及优点，利用它开发遥感图像处理的算法，能极大地提升计算效率，简化人为利用软件处理图像的步骤，也能根据不同的具体应用需求，开发不同的遥感图像处理算法。

2.3 植被长势参数反演

在植被光谱特征中，由于植被的叶绿素在红色波段表现出较强的吸收能力，而植被叶片的结构特征导致在近红色波段的辐射吸收较少，形成较强的反射。因此，以数学组合形式，构建以红波段、近红波段为主的植被指数，可以反映出植被的长势情况[5-6，20]。本文采用阴影消除植被指数（SEVI）进行植被长势信息提取，这是因为在某些复杂地形山区，太阳辐射的传输易受到地形效应的影响。SEVI能够有效消除崎岖地形中本影和落影的影响，获得准确的植被长势信息[21-22]，计算公式如下：

式中，Bnir为近红外波段反射率；Bred为红色波段反射率；Bnir/Bred为比值型植被指数RVI；1/Bred为阴影植被指数SVI；f(Δ)为地形调节因子，通过控制它来调节阴坡、阳坡的植被指数值，调节计算方法如下：

其中，x是SEVI,y1是RVI,y2是SVI,n是子区域图像中的像素数。设定f(Δ)从0变为1，固定步长为0.001。如式（4），当两个相关系数趋近相等时，确定为最优调节因子。

2.4 变化检测

2.4.1 图像差值法

图像差值法是一种简单快速的变化检测方法[9-12]，它的原理是：图像中未发生变化的地类在两个时相的遥感图像上一般具有相近或相等的灰度值，而当地类发生变化时，对应位置的灰度值将有较大差别。因此在差值图像中，没有发生变化区域的像素值接近0，明显发生变化区域的像素值会明显区别于0，从而使变化信息从图像中显现出来。

因为不同时期的SEVI 受到的地形阴影影响不同，采用了不同的地形调节因子，所以需要进行不同时期图像的标准化处理。假设同一自然保护区同一季相的不同年份的遥感图像的植被指数最大值相近，而在岩石具有相近的最小值。这些最大值和最小值可以用于不同年份SEVI数据之间的标准化，公式如下：

式中，SEVI′是标准化的SEVI值；SEVImin是预先判定为某块岩石处SEVI的值；SEVImax是预先判定为某块植被长势最好处SEVI的值。

式中，E(t1-t2)代表从t1到t2时间段内SEVI的变化量；SEVIt1代表t1时相的植被指数；SEVIt2代表t2时相的植被指数。

2.4.2 阈值的选取

本文采用最大类间方差法（OTSU）来自动获取阈值。这是一种被公认的阈值自动选取方法的最优方法之一，具有算法简单、分割速度快等优点[14，23-25]。阈值的自动获取方式较人机交互方式更加智能化，且避免了因人的主观判断所带来的误差。最大类间方差法（OTSU）的基本思想是如果一幅图像由一物体和背景构成，物体与背景有不同的灰度值，以差异图像中的某一灰度为阈值将图像分为目标和背景两组并计算两组间的方差，当被分成的两组之间的方差最大时，就以这个灰度值作为阈值分割图像[25]。

假设图像的大小为M×N，图像中像素的灰度值满足（大于或小于）阈值T的像素个数记作N1，像素值不满足（大于或小于）阈值T的像素个数记作N2，则有：

式中，ω1为当前满足阈值的像素点数占整幅图像的比例，其平均灰度为μ1；ω2为当前不满足阈值的像素点数占整幅图像的比例，其平均灰度为μ2；xi、xj为遥感图像的像素值。

式中，g为类间方差，按照图像的灰度特征，采用遍历的方法得到使类间方差最大的阈值T。

2.5 利用K 均值聚类分割与识别变化区域

大部分的地理现象都具有空间相关特性，即距离越近的两事物越相似。因此，植被长势的变化点位在一定的区域范围内往往符合聚集规律，但是整体仍会呈现出不均匀、不规则的分布状态，甚至在某些区域还会出现异常点。本文根据所有已检测到点位的行列位置，利用K均值聚类算法将这些点位分为不同区域的簇类，然后再根据每个簇类中所有点位的行列位置，绘制出更具直观效果的形状，来辅助地物判断与分析。

K均值聚类的基本思想是：对于给定的聚类数目k，首先随机创建一个初始分类，然后采用迭代方法通过聚类中心的不断移动来尝试着改进划分[26-30]。具体公式如下：

式中，E是所有样本的聚类误差；Ci是第i个簇；x是Ci中的样本点；μi是Ci中的质心。

这种聚类算法简单且高效，但其聚类结果的优劣往往取决于预先设定的聚类数目k。为了提高聚类稳定性并降低时间复杂度，本文选择平均轮廓系数来自动确定聚类数目k，提高自动化程度。具体公式如下：

式中，ai是当前第i个样本与同簇其他样本的平均距离，称为凝聚度；bi是当前第i个样本与最近簇中所有样本的平均距离，称为分离度。Si是当前第i个样本的个体轮廓系数，当求出所有样本的轮廓系数后再求平均值就得到了平均轮廓系数[31]。平均轮廓系数的取值范围为[-1,1]，且簇内样本的距离越近，簇间样本距离越远，平均轮廓系数越大，聚类效果越好[32-34]。

结合以上原理及式（10）、（11）可以找出最优的聚类簇数Kopt，同时会产生Kopt个聚类簇。此时，为了识别每一个簇所在的区域，可以根据每一个簇中所有点位的行列位置来提取这些簇的外接图形。其中，外接矩形具有较好的目视效果，且相对其他图形比较规整，更容易看出当前簇的形状大小、内部的点位密度等信息。因此，本文选择绘制这些簇的外接矩形，具体方式如下：

3 实验与结果分析

3.1 植被长势参数反演结果

为验证遥感图像变化检测技术的效果，本文选择了Landsat8-OLI 的两幅遥感图像作为数据源（表1），裁剪边界为福建省武夷山自然保护区（图2），该保护区内地形崎岖复杂，最小坡度为0°，平均坡度为27.3°，最大坡度大于60°，标准方差为10.04。地物类型主要是植被，同时包含了一些水体、裸地和建筑物。

依据式（1）～（5）得到了上述两幅遥感图像的SEVI标准化结果（图3），从目视效果来看，SEVI 的标准化结果，整体呈现出平坦的特征，基本没有图2 中山体区域出现的浮雕情况。这与Jiang 等[21]发现的SEVI 呈现的目视效果一致。

3.2 差值检测结果

依据式（6）～（9）得到了两幅遥感图像的变化检测结果（图4 和表2）。已知类别1 是植被长势明显减少的区域，在研究区内占120个像元（0.108 km2）；类别2是植被长势明显增加的区域，在研究区内占1 659 个像元（1.490 km2），是明显减少区域的13 倍左右。两种明显发生变化的像元点位共约占所有像元数的0.10%。

表1 两幅Landsat8-OLI遥感图像的相关参数

图2 Landsat8-OLI遥感图像

图3 SEVI的标准化结果

图4 变化点位

表2 OTSU法阈值检测结果

从图4可以发现，植被长势发生明显变化的地区主要集中于研究区中部偏北区域，呈现小面积聚集状态。

3.3 利用K 均值聚类的识别结果

根据图像差值以及OTSU 阈值提取的明显变化点位进行聚类分析，结果如图5所示。

图5 变化点位的行列号

根据研究区内变化点位的数量以及分布情况（图5），设定聚类簇数K以1 为步长，从1 到60，计算当前聚类簇数K所对应的平均轮廓系数（式（11）），并绘制两者的关系曲线（图6）。可以发现，当聚类簇数为15 的时候，类别1 有最高的平均轮廓系数为0.649 7。而类别2的曲线呈双峰分布且右边的峰值比左边的峰值更大，所以取对应的聚类簇数为56。

图6 平均轮廓系数

依据确定出的Kopt个聚类簇以及每簇中所有点位的行列号，结合式（12）、（13），就可以绘制出所有聚类簇的外接矩形（图7）。通过与图4的对比可以发现，图7具有更明显的目视效果，提高了检测精度。

图7 识别结果

3.4 不同区域的测试结果

为了充分验证本文技术的有效性，选择福建省闽江源国家级自然保护区作为测试研究区（图8）。该保护区内的地形同样崎岖复杂，海拔在250～1 858 m之间，平均坡度为14.2°，标准方差为8.45，包含了植被、裸地、建筑物、阴影这些主要地物类别，相比武夷山自然保护区，它的地物分布情况较为错综复杂。选择的遥感图像详细信息可见表3。

图8 Landsat8-OLI遥感图像（图像是753波段的假彩色合成）

表3 两幅Landsat8-OLI遥感图像的相关参数

图9 结果展示

通过本文技术检测和识别的结果（图9 和表4），可以得出以下结论：（1）植被长势明显增加的点位主要分布在两个区域，分别是闽江源的最北与最南部地区，面积较大，共占2 260个像元，约2.034 0 km2。（2）植被长势明显减少的点位比较少，仅占25 个像元，分布在6 个不同的区域，共约0.022 5 km2。（3）图9（b）的目视效果要优于图9（a），更能直观地看出植被长势增加的有两块区域，且南部区域的密度要整体大于北部区域，而且突出了减少点位区域的目视效果（如需要更加突出，则仅需绘制时加粗线条）。

表4 闽江源国家级自然保护区两幅遥感图像的检测结果

4 结束语

针对遥感图像变化检测的自动化程度较低，现有算法不易扩展等问题，本文提出一套基于GDAL的遥感图像变化检测技术，主要包括：利用阴影消除植被指数（SEVI）反演植被长势；利用图像差值法及最大类间方差法（OTSU）来提取植被长势明显变化点位；利用K均值聚类自动分割并识别变化区域。实验结果表明，本文方法能够有效地检测植被长势变化区域，可为区域尺度的植被长势监测提供技术支持。下一步将研究更多阈值自动化提取算法、分割算法在遥感图像变化检测中的应用。