土地利用变化图斑提取方法

夏 旺，赵 展

（1.武汉大学 测绘学院，湖北 武汉430079）

土地利用变化图斑提取方法

夏 旺1，赵 展1

（1.武汉大学 测绘学院，湖北 武汉430079）

为提高变化图斑提取效率，减少人工操作难度，在多元变化检测（MAD）自动发现变化信息的基础上，提取感兴趣区域生成差异影像，并分割差异影像得到完整变化图斑。利用3种典型变化地物类型进行图斑提取实验，结果表明，该方法可以获取完整﹑连通的变化图斑，在错检率和漏检率方面优于MAD，且Kappa系数分别从0.770﹑0.810﹑0.729提高到了0.916﹑0.894﹑0.934。

变化图斑；分水岭分割；区域合并；变化检测

获取两幅不同时相遥感影像的变化图斑是土地利用遥感监测的关键[1]，但目前通常仍是采用人工勾绘的方式来获取变化图斑，我国每年土地变化信息量十分巨大，手工提取变化图斑耗时耗力，且准确勾绘边界线也对操作人员的技术水平提出了较高的要求，因此减少人工操作难度具有重要意义。

利用变化检测方法可以自动识别不同时相遥感影像的变化区域。根据处理单元可将检测算法分为像素级和对象级两类[2-3]。对象级变化检测是将一组有特定关联的像素集合为一个检测单元来进行变化检测，可对上下文信息建模，有效考虑局部邻接像素之间的关系，能较好地去除由光谱差异或配准误差造成的伪变化区域[4]，但该算法依赖于分割与分类算法，细小的变化区域可能被湮没在背景中而漏检[5]，同时由于光照﹑视角﹑大气条件等影响，前后两个时相影像的分割与分类结果难以保持一致，将导致大量伪变化区域[6]，因此在土地覆盖﹑土地利用调查更新等多个应用领域中，仍广泛使用像素级变化检测。像素级变化检测最大的优点是可以全面检测影像中任意形状﹑大小的变化信息，达到很好的效果；但由于算法没有考虑邻接像素的信息，容易对由光谱差异和配准造成的误差敏感，产生离散﹑破碎﹑不连通的检测结果[7]，且存在大量伪变化区域，不能直接利用变化信息提取变化图斑。

针对变化检测算法结果离散﹑破碎﹑不连通的缺点，本文首先利用MAD结果辅助人工选取感兴趣区域（ROI），再利用ROI切割前后时相影像获得影像切片，最后对切片进行归一化再相减获得差异影像。差异影像每个像素的像素值表征该像素的变化程度。利用分水岭算法分割差异影像，并针对分水岭算法过分割问题，设计了一种二段的区域合并方法。

1 基于MAD的差异影像获取

MAD是基于典型相关分析进行变化区域检测的方法。为了有效集中差异信息，提高检测精度，MAD采用计算一对典型变量并相减的方式，既可在最大程度上消除不同时相﹑不同通道相关信息的影响，又可将所有变化信息分配到互相独立的MAD变量中去[8]。

对MAD变量进行后处理和阈值化[9-10]可得到一 幅二值图像。由于MAD结果中存在大量伪变化区域，不能直接用于获取ROI。本文通过人工目视观测MAD提取的变化结果，确定所提取的变化区域是否为真实变化，并选取包含变化区域的外围矩形为ROI，利用ROI切割前后时相影像，获得影像切片Cut1与Cut2。

虽然用于变化检测的影像一般会先经过辐射畸变校正，但辐射畸变并不能被完全消除，且由于两个时相影像可能来自不同传感器，获取时刻光照强度也不完全相同，两个影像仍存在辐射差异。为了尽可能消除辐射差异的影响，本文先将Cut1与Cut2所有波段分别归一化至0～255，归一化后的切片分别记为NCut1与NCut2，再将NCut1与NCut2对应波段相减并取绝对值获取差异影像Dif（Dif的第j个波段记为Difj）：

2 差异影像的分水岭分割

本文利用分水岭算法对差异影像进行分割，分水岭算法对边缘具有良好的响应，可得到封闭连续边缘。差异影像通道数与原始影像相同，往往是多通道。本文对差异影像的每个波段分别作分水岭分割，再将不同波段的分割结果合并在一个波段上。

对于Difj，先利用Sobel算子求得梯度影像Gj。为了便于分水岭算法操作，需将像素值归一化至0～255。分水岭分割的关键在于种子点的选取，本文将Gj所有的局部最小值作为种子点，即遍历Gj每个像素，判断该像素在其3×3的邻域内是否为唯一极小值，若是，则该点为种子点。获取种子点后就可对Gj作分水岭分割[11]，得到分割结果Segj。Segj是一个标记影像，将Segj像素p的标记值记为Segj(p)。

对Dif各波段都求取相应的分割结果后，需将所有分割结果合并为一个新的分割结果影像，记为Seg，它保留了各波段的分割结果，这是由于变化信息可能分布在不同波谱上，所以需在合并时保留所有变化信息。合并的具体步骤为：

1）将Seg所有像素设为0，标记值记为label，初始label记为1；

2）顺序遍历图像，当前遍历像素为p，若p=0，则标记为label；

3）遍历p的4个邻接像素，令当前遍历的邻接像素为s，若s≠0，且对于任意一个波段j都满足Segj(s)=Segj(p)，则s被标记为label，并将s压入堆栈；

4）从堆栈中取出一个像素作为p返回步骤3）；

5）当堆栈为空时，label加1，返回步骤2）。

3 分割区域合并

分水岭算法对噪声敏感，同时由于种子点的选取非常密集，分割结果会出现过分割现象。针对这两个缺点，本文设计了一种二段区域合并的方法。

对Seg的每个分割区域计算区域面积（区域像素个数）以及各波段的像素平均值。若影像有N个波段，则第i个分割区域各波段的像素平均值为meani1, meani2,…, meaniN。

设区域i与区域j的光谱差异RegionDifij为：

若RegionDifij较大，则区域i与区域j的光谱差异较大。

3.1 根据面积合并

本文选取的种子点较为密集，分割结果将存在严重的过分割现象。分水岭算法对误差十分敏感，由于小面积区域的像素个数较少，噪声对分割区域像素平均值的影响较大，所以需先将小面积区域与邻接区域合并，若该区域存在多个邻接区域，则选择与它光谱差异最小的邻接区域合并。具体算法为：遍历Seg的所有分割区域，令当前遍历的分割区域为Reg，若Reg的像素个数小于某一阈值A，则根据式（2）计算该区域与它所有邻接区域的光谱差异RegionDif，选取RegionDif最小的邻接区域与Reg合并。阈值A的大小与影像大小有关，本文实验所用A的大小为Seg的像素总数除以500。

3.2 根据光谱差异合并

通过面积合并后，误差都淹没在较大面积的分割区域中，所以对分割区域像素平均值影响较小。此时可计算每个分割区域与各自邻接区域的光谱差异性，若光谱差异较大，则代表两个区域不属于同一类别；若较小，则将两个区域合并。可利用RegionDifij大小判断区域i与区域j是否应被合并为一个区域，需寻找合适的阈值B，若RegionDifij＜B，则将区域i与区域j合并为一个区域。

不同的差异影像的变化区域与非变化区域的反差并不相同，所以每个Seg的最佳合并阈值也不相同。本文根据K均值聚类算法的思想，寻找自适应的阈值B。先利用式（2）计算面积合并后剩下的所有区域与其邻接区域的RegionDif，得到所有区域与各自邻接区域RegionDif的集合，再利用K均值聚类算法处理RegionDif的集合。本文设K=2，将RegionDif的集合分为两类：一类中的RegionDif值较小，代表该类中RegionDif值对应的两个分割区域光谱差异较小，应合并为一个区域，将该类中最大的RegionDif值记为Rmax；另一类中的RegionDif值较大，代表该类中RegionDif值对应的两个分割区域光谱差异大，不应该合并，将该类中最小的RegionDif记为Rmin，则阈值B为：

若RegionDifij＜B，则代表相邻的两个区域的光谱值接近，可合并为一个区域。

3.3 变化图斑提取

根据光谱差异合并后，影像会被划分为少量几个区域，每个区域内像素值接近，区域与区域之间像素值差异较大。若区域的像素平均值较大，则该区域为变化区域；反之，则为非变化区域。令区域i所有像素所有波段的均值为：

可以找到合适的阈值C，使得若meani>C，则将区域i判定为变化区域。但是不同的差异影像会有不同的最佳阈值，所以此处仍使用K均值聚类算法来寻找自适应的阈值C。

计算合并后剩下的各区域的mean，对得到的所有集合作K均值聚类，其中K=2。将mean的集合分为两类：一类中的mean值较小，代表该类mean所对应的区域为非变化区域，将该类中最大的mean记为Mmax；另一类中的mean值较大，代表该类mean所对应的区域为变化区域，将该类中最小的mean记为Mmin，则阈值C为：

若meani>C，则区域i为变化区域；反之，则为非变化区域。

4 实验与结果分析

4.1 实验数据

本文实验采用的遥感影像数据为2012年﹑2013年杜尔伯特蒙古族自治县地区的SPOT卫星影像，包含红﹑绿﹑蓝3个波段，分辨率为2.5 m。两期影像的地物类型丰富，主要由水体﹑建筑物﹑道路﹑植被等多类地物构成，且影像的分辨率较高。对实验数据进行MAD变换，得到MAD结果。目视判别MAD结果，选择真实变化区域的外围矩形作为ROI，并利用本文算法获取最后变化图斑。

本文选择了土地利用遥感监测中常见的3种典型变化地物类型，处理结果如图1～3所示。图1a为人工勾勒的变化图斑，作为真实值对算法进行精度评定；图1b为ROI上生成的差异影像；图1c为ROI对应的MAD结果；图1d为利用本文算法得到的结果。对MAD算法和本文算法进行精度评定，结果如表1所示。

4.2 实验分析

图1c中MAD检测的道路不完整，上部分出现断裂，同时存在离散的噪声点；图2c中MAD检测建筑中存在空洞，同时由于配准误差，上部分未变化的曲线道路被检测为变化区域；图3c中MAD检测结果更加离散﹑破碎，其中细长部分的建筑断开导致区域不连通，同时存在空洞，且左上部分的建筑屋顶由于辐射差异导致部分像素被MAD判定为变化区域。综上可知，MAD可有效发现变化，但检测的结果离散﹑破碎﹑不连通，形状不够完整，同时存在伪变化区域，不能直接用于图斑提取。

图1～3的差异影像充分反映了前后时相影像的变化信息，影像中较亮的部分代表变化区域。观察MAD检测结果和差异影像可以发现：①由配准误差造成的伪变化区域在差异影像上的像素值仍较亮，但分布离散，面积较小；由辐射差异造成的伪变化区域像素值虽比背景区域的像素值高，但比真实变化区域的像素值要低。②MAD确定的变化区域内存在空洞部分，在差异影像上偏暗。

图1 植被-道路类型变化图斑提取结果

图2 植被-建筑类型变化图斑提取结果

图3 植被-推填土类型变化图斑提取结果

表1 变化检测精度评定

从图1～3的分割结果可以看出，本文算法可剔除MAD上由于辐射差异﹑配准误差等因素造成的误差。因为在区域合并时，首先将小面积区域和周围区域合并，而由于配准误差形成的伪变化区域分布离散，面积小，所以会被周围较暗的背景区域合并。辐射差异造成的伪变化区域像素值低于真实变化区域，在利用区域间差异性进行区域合并时，会在K均值聚类中被判定为背景。从分割结果中还可以看出，变化区域内部的空洞被填充，是一个完整的变化图斑，破碎﹑不连通的部分被合并为一个整体。这是由于在第一次区域合并时，空洞对应的较暗区域与周围较亮区域合并，整体呈现变化特征，在最后的K均值聚类中被判定为变化区域。

从表1可以看出本文算法漏提与错提区域明显小于MAD算法，分割结果与人工提取的结果相似，3种检测结果的Kappa系数为分别从0.770﹑0.810﹑0.729提高到了0.916﹑0.894﹑0.934。

5 结 语

常用的变化检测方法检测得到的变化区域离散﹑破碎﹑不连通，且存在大量伪变化区域，不能直接用于图斑提取。本文设计的影像分割与合并算法有效地分割了差异影像，得到了完整的变化图斑，解决了MAD结果离散﹑破碎﹑不连通的问题，在错检率和漏检率方面优于MAD算法。本文算法避免了人工勾勒变化图斑边界，可有效提高变化图斑提取效率，减少人工操作难度。下一步需要研究ROI自动化提取办法，实现全自动化变化图斑提取。

[1] 王素敏, 翟辉琴. 遥感技术在我国土地利用/覆盖变化中的应用[J].地理空间信息,2004,2(2):31-32

[2] 吴芳,刘荣,田维春,等.遥感变化检测技术及其应用综述[J].地理空间信息,2007,5(4):57-60

[3] Tewkesbury A P, Comber A J, Tate N J, et al. A Critical Synthesis of Remotely Sensed Optical Image Change Detection Techniques[J].Remote Sensing of Environment,2015(160):1-14

[4] Volpi M, Tuia D, Bovolo F, et al. Supervised Change Detection in VHR Images Using Contextual Information and Support Vector Machines[J].International Journal of Applied Earth Observation and Geoinformation,2013,20(2):77-85

[5] Listner C, Niemeyer I. Object-based Change Detection[J]. Photogrammetrie Fernerkundung Geoinformation,2011(4):233-245

[6] CHEN X, CHEN J, SHI Y，et al. An Automated Approach for Updating Land Cover Maps Based on Integrated Change Detection and Classification Methods[J].ISPRS Journal of Photogrammetry and Remote Sensing,2012,71(326):86-95

[7] 佃袁勇,方圣辉,姚崇怀.一种面向地理对象的遥感影像变化检测方法[J].武汉大学学报(信息科学版),2014,39(8): 906-912

[8] Bruzzone L, Prieto D F. Automatic Analysis of the Difference Image for Unsupervised Change Detection[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1 171-1 182

[9] Nielsen A A. The Regularized Iteratively Reweighted MAD Method for Change Detection in Multi- and Hyperspectral Data[J]. IEEE Transactions on Image Processing a Publication of the IEEE Signal Processing Society,2007,16(2):463-478

[10] GONG P. Change Detection Using Principal Component Analysis and Fuzzy Set Theory[J].Canadian Journal of Remote Sensing,1993,19(1):22-29

[11] 陈波,张友静,陈亮.标记分水岭算法及区域合并的遥感图像分割[J].国土资源遥感,2007,19(2):35-38

P237

B

1672-4623（2017）09-0093-04

10.3969/j.issn.1672-4623.2017.09.028

2016-06-28。

项目来源：国土资源部公益性行业科研专项资助项目（201511009-01）。

夏旺，硕士研究生，主要研究方向为摄影测量与遥感。