基于多源数据的土地利用变化检测

2018-12-20赵展夏旺闫利

自然资源遥感 2018年4期

赵展，夏旺，闫利

(武汉大学测绘学院，武汉 430079)

0 引言

我国每年都进行年度土地利用变更调查，以保证我国土地利用调查数据成果具有很好的现势性。目前的土地利用变化信息获取还是以人工解译遥感影像为主，任务繁重，效率较低。利用计算机自动分类技术快速获得土地利用变化信息已成为当前的迫切需求。

根据变化检测前是否获得土地利用/覆盖的类别信息，可以将变化检测方法分为分类后比较法和直接比较法。直接比较法是直接对同一区域不同时相遥感影像的光谱和其他特征信息进行比较，利用差异信息直接确定变化的位置与范围，常用方法有影像代数法[1]、变化矢量分析法[2-3]、多元变化检测(multivariate alteration detection，MAD)变换法[4]以及光谱梯度差分法[5]等。在这类方法中，除了光谱特征，距离测度[6]、自相关性[7-8]和建筑物指数[9-11]等也被广泛使用，在直接比较结果的基础上设置阈值或进行二类分类从而获得变化区域[2-3，12-13]。分类后比较法是先获得前后时相的土地利用/覆盖类别，通过对比前后类别发现变化，常用的方法有基于像元的随机森林分类[14]、贝叶斯分类[15]、马尔科夫随机场集成多特征分类[16-18]，以及面向对象分类[19-20]等。分类后比较法的关键问题在于分类本身是一个比较繁琐的过程，且易于将分类误差引入变化检测结果并放大。目前也有学者开始研究利用土地利用数据库辅助变化检测，如将后一时相分类结果与土地利用矢量对比发现变化[19,21]、从土地利用数据中学习规则知识[22]等。

尽管这些方法在国内外取得了较好的变化检测效果，但往往过程较复杂，需要人工设置一定的参数或规则，并存在较多的漏提。土地利用变更调查的原则是尽量不遗漏任何疑似变化区域，同时由于作业范围大，要求方法具有一定的自动化程度和较快的运行速度。针对这些需求，本文提出一种基于多源数据的土地利用变化自动检测方法。该方法利用多源数据实现信息互补，将分类后比较法和直接比较法相结合，以提高变化检测的自动化程度、效率和可靠性。

1 研究区概况与数据源

本文研究区选择了我国南方和北方2个生态经济区的典型区域。一个位于黑龙江省杜尔伯特蒙古族自治县，属于东北平原地区。该县近些年来以“畜牧强县”政策为主发展经济，大力发展畜牧业，并且鼓励牧民进城镇生活，出现了大量新增建设用地。研究区域为该县政府驻地的康泰镇及其周围区域，面积约为150 km2。另一个位于湖南省长沙县，属于南方丘陵地区。该县处于长株潭“两型社会”(资源节约型和环境友好型社会)综合配套改革试验区的核心地带，近年大力发展工业、建立物流基地，同样出现大量土地利用变化。区内包括黄兴镇和干杉镇范围，面积约为160 km2。本文使用的数据为国家年度变更调查标准数据，包括前、后2时相的遥感影像，以及前一时相的土地利用数据。同时采用后一时相国家土地利用变更调查中人工解译提取的变化图斑结果作为对比数据，该结果是作业人员在ArcGIS软件中通过目视对比前后时相影像提取的。其中，杜尔伯特蒙古族自治县研究区数据为2012和2013年度的数据，遥感影像是SPOT5影像(空间分辨率为2.5 m)；长沙县研究区数据为2013和2014年度的数据，遥感影像为国产高分卫星影像(空间分辨率为2 m)。所使用的遥感影像均按照国家年度变更调查标准技术规范进行了几何纠正和辐射校正处理，与土地利用数据精确配准。研究区域遥感影像分别如图1和图2所示。

(a) 2013年影像 (b) 2014年影像叠加人工解译变化图斑

2 基于多源数据的土地利用变化检测

本文提出利用多源数据进行自动的土地利用变化检测。首先,利用前一时相的土地利用数据，自动获取后一时相影像上的训练样本，进行影像分类，将后一时相分类结果与前一时相土地利用数据对比，自动地发现变化区域；然后，对前、后2时相影像做MAD变换，利用MAD变换影像自动设置阈值去除伪变化；最后，通过目视判别确定提取的疑似变化图斑是否为真实变化，对于真实变化区域采用人工绘制边界方法提取变化图斑。整体技术流程如图3所示。

图3 变化检测与提取整体技术流程

2.1 矢量数据辅助影像自动分类

影像分类是一个比较繁琐的过程，特征选择、样本提取和分类器参数设置等方面均影响分类结果，对操作人员经验要求较高。本文基于第二次土地利用调查中农村土地利用类别体系建立分类类别体系，研究从已有土地利用矢量数据中自动获得样本的影像全自动分类方法，并在此基础上实现自动的变化检测。

2.1.1 自动分类类别体系设计

年度土地利用变更调查重点监测的是非建设用地到建设用地的变化，包括已建成区域和在建区域，其规定的变化类型与土地利用调查土地分类体系中的类别并不完全对应一致。在内业变化图斑提取过程中，只需要发现变化区域，最重要的需求是不遗漏真正变化区域和提取效率尽可能高，最终的变化后地物类别是通过外业核查确定的。针对这一特点，本文将土地利用调查土地分类体系中的地物类别简化为5大地类进行影像分类，将“耕地”、“园地”、“林地”和“草地”简化为植被类；“交通运输用地”中除“机场用地“和“港口码头用地”外，均简化为道路类；“水域及水利设施用地”中除“水工建筑用地”外，均简化为水体类；“其他土地”中除“设施农用地”外，均简化为裸地类；“城镇村及工矿用地”、“机场用地”、“港口码头用地”、“水工建筑用地”和“设施农用地”合并为建筑物类。

2.1.2 样本的自动选择与精化

样本的选择对分类结果有较大的影响。传统方法中分类样本是通过人工选取的，选取的结果严重依赖于作业员的专业知识和经验，使分类结果受作业员的主观性影响较大。在本文方法中，利用已有土地利用数据自动获得训练样本并精化，从而实现全自动的影像分类，在提高效率的同时保证了分类结果的客观性。

在土地利用变更调查中已有前一时相土地利用矢量数据，尽管后一时相部分区域土地利用会发生变化，但在未发生变化的区域，矢量数据的类别和影像是正确对应的。实际情况中，1 a内变化区域相对只是少数，土地利用数据中图斑的类别信息可认为绝大部分都是准确的。因此，将前一时相土地利用数据中的所有图斑作为初始样本，采用迭代过程剔除少数错误样本，精化初始样本。其步骤为：首先，计算每个类别的类别特征中心，计算每个样本与各个类别特征中心的特征距离，选择所有特征距离最小类别与土地利用数据类别一致的样本，组成初始的精化样本集；然后，利用初始的精化样本集训练分类器对每个样本进行分类，对比每一个样本分类类别和在土地利用数据中的类别，将不一致的样本剔除，得到新的精化样本集；最后，利用新的样本集重新训练分类器，如此循环迭代，直至所有样本的分类类别与土地利用数据中的类别一致，得到最终的精化样本集，并采用训练得到的最终分类器对影像进行分类，得到最终分类结果。方法的基本流程如图4所示。

图4 样本精化流程

2.1.3 分类器设计与分类特征计算

支持向量机(support vector machine，SVM)遵循结构风险最小化的原则，具有强大的非线性和高维处理能力及泛化能力，非常适合于遥感影像的分类。大量实验表明，在SVM的各种核函数中，径向基函数(radial basis function,RBF)具有较好的稳定性。因此，本文选择基于RBF的SVM作为影像分类器。通过大量实验，设置RBF的γ参数为2.0，惩罚系数C为1.0。

本文采用的分类特征是影像的光谱特征和光谱特征的归一化四则运算特征。归一化四则运算特征计算方式为

f=(pixelValuei-pixelValuej)/(pixelValuei+pixelValuej)

(1)

式中：i,j表示波段，满足i,j=1,2，...,n(n为影像波段数),且i

完成影像分类后，将分类结果与土地利用矢量对比即可获得变化检测结果。土地利用变更调查主要关注非建设用地变化为建设用地。因此，将类别不一致且分类类别为建筑物、道路或裸地的像素标记为变化像素，将相邻的变化像素聚集成疑似变化图斑。尽管在分类过程中不可避免会存在分类误差，但本文方法的主要错误分类是将非建设用地误分为建设用地，因此，变化检测的结果一般不会出现漏提现象，满足变更调查不遗漏变化的要求。另一方面，变化检测结果会存在较多的伪变化，需要进一步去除，本文利用前后时相影像的MAD变换影像来去除伪变化。

2.2 基于MAD变换的伪变化去除

2.2.1 MAD变换原理

MAD变换的基本思想是认为前后时相影像在未变化的区域存在强相关，在变换区域存在弱相关。通过典型相关分析可以分析前后时相影像的整体相关性突出变化区域。

设随机向量X=[X1,X2,...，Xp]T和Y=[Y1,Y2,...，Yp]T分别是拥有p个波段的不同时相影像。将X和Y分别做零均值归一化，然后做线性变换，即

(2)

式中：a=[a1,a2,...,ap]∈Rp,b=[b1,b2,...,bp]∈Rp为线性变换系数矩阵。则MAD定义为

(3)

根据典型相关分析原理，2个新的变量U和V之间的相关性越强，两者之间包含的变差信息越少；反之，相关性越弱，所包含的变差信息越多。所以在归一化约束条件Var{aT·X}=Var{bT·Y}=1下，以Var{aT·X-bT·Y}最大化为准则，找到变换系数矩阵a和b即可完成MAD变换；然后根据以方差倒数为权值的各波段均方根可以得到最终的MAD变换影像。在MAD影像上变化区域像素一般具有较高的灰度值，而非变化区域具有较低的灰度值。

2.2.2 利用MAD变换去除伪变化

尽管现有研究表明，MAD变换相比于其他正交变换(如主成分分析变换)具有更好的突显变化的能力[4],但当影像覆盖范围大、地物复杂时，MAD影像上变化区域和非变化区域灰度值区间仍然存在一定重叠，设置阈值提取变化区域时，要满足变更调查“不遗漏变化”的要求，会造成较多的误提伪变化。因此，本文并不直接利用MAD影像提取变化区域，而是在自动分类变化检测基础上，利用MAD影像剔除伪变化区域。

大范围检测任务下，影像上变化区域所占比例较小。因此本文通过灰度直方图自动确定阈值，即将直方图中一定百分比对应的灰度值作为阈值，小于阈值的像素被认为是一定没有发生变化的像素，如果一个疑似变化区域所有像素的灰度值均小于阈值，则认为是伪变化区域，将其剔除。另外，高空间分辨率影像上地物局部细节的变化可能会造成细小的伪变化，以及前后时相存在几何配准误差也会造成细小的伪变化，这些伪变化可以通过设置面积阈值滤除。

2.3 变化图斑绘制与提取

实际作业中，一般在伪变化去除后仍然会有一定的伪变化存在。这些伪变化除了来源于分类错误之外，更多的是客观地物与土地利用定义不一致造成的。例如农田由于搭建农业大棚等呈现类似建筑物的特征被误分为建筑物，盐碱地等也会呈现类似建筑物的光谱特征；而田间的道路由于一般宽度较窄在土地利用分类中仍被归为农田类(旱地、水浇地和水田等)的一部分。这些伪变化很难通过自动方法滤除。因此本文方法采用目视观察剔除这类伪变化，然后通过目视对比确定变化区域的变化类型，最后人工绘制准确的变化图斑边界。

3 实验与分析

采用本文方法对杜尔伯特蒙古族自治县研究区域数据进行处理，结果如图5所示。

(a) 类别简化的2012年土地利用数据 (b) 2013年影像自动分类结果

(e) 去除伪变化后的变化检测结果 (f) 最终提取变化图斑

首先将2012年度土地利用数据进行类别简化(图5(a))，然后通过样本精化自动获得训练样本，利用自动训练的分类器对2013年影像进行分类，结果如图5(b)。从图中可以看出，自动分类结果存在一定的误差，主要的误分区域是部分亮水体和道路被误分为建筑物。其中道路和建筑物都属于建设用地，在土地利用变化检测过程中，道路的误分不会造成影响，而水体的误分会造成误提的伪变化，需要通过后续处理去除。另一方面值得注意的是，2013年影像成像时间在雨季，地面积水比较严重，在2012年影像上部分草地和滩涂裸地被水淹没，2013影像的分类结果上表现出较多的水体，这并非真实的变化，而本文变化检测方法中只关注新增建设用地，这一现象不会造成建设用地变化检测的漏提或误提。将分类结果与2012年土地利用结果对比，得到的初步的变化检测结果如图5(c)所示。从图5(c)中可以看出，尽管自动分类存在一定误差，但其错误主要造成伪变化，而不会漏提真实变化区域。初步提取的疑似变化图斑数量较多(23 735个)，需要进一步去除伪变化。首先利用变化区域的面积大小去除过小的伪变化。我国的土地利用规范要求平原地区最小图斑面积为600 m2，影像的空间分辨率为2.5 m，因此将面积阈值设为90个像元。然后计算2012年和2013年影像的MAD变换影像，如图5(d)所示。从图中可以看出发生变化的区域在MAD影像上亮度值都相对较高，可以通过设置阈值去除一定的伪变化。如上文所述采用累积直方图确定阈值，本研究中累积直方图比例设为98%。进行伪变化去除后结果如图5(e)所示，去除伪变化后剩余疑似变化图斑减少到512个。剩下的疑似图斑中仍有比较多的伪变化，这些伪变化除水体外，主要是裸露的农田和农田间的道路，如上文所述，大部分伪变化是客观地物与土地利用定义不一致造成的，只能通过目视观察分辨。这些伪变化往往比较集中，可以通过目视观察快速剔除。经过目视判别后，最终提取并绘制变化图斑37个，结果如图5(f)所示。

采用同样的处理流程对长沙县数据进行处理，取得类似的处理结果，该研究区域地物细碎且变化较多，最终提取变化图斑156个，如图6所示。

(a) 自动变化检测结果 (b) 最终提取变化图斑

将本文方法提取的结果与国家年度土地利用变更调查结果对比。对于杜尔伯特蒙古族自治县，国家土地利用变更调查共提取了24个变化图斑，本文方法提取了37个变化图斑；对于长沙县，国家土地利用变更调查共提取了139个变化图斑，本文方法提取了146个变化图斑。相比之下，本文方法能够提取更多的变化图斑，更不易遗漏变化区域。图7示出本文方法与传统方法提取结果的局部对比。

如图7(a)和(b)所示影像1中，一个新增加的绿色建筑物在年度变更调查中被遗漏，而本文方法能够将其提取出来。除此之外，年度土地利用变更调查方法只对比前后时相影像，对于历史原因造成的土地利用数据与实际情况不一致则无法发现。如图7(c)和(d)所示影像2中，土地利用矢量中类别是旱地，但在2012年影像上已经表现出建设用地特征，二者不一致，这是历史原因造成的土地利用矢量与实际情况不符。由于该区域在2012和2013年影像上均表现为建设用地，仅目视对比影像认为没有变化发生。而本文方法同时利用土地利用矢量和前后时相影像，因此能够发现这一历史原因造成的没有及时更新的变化。

在实际应用效率方面，本文方法在自动变化检测基础上提取变化图斑，相比于现有目视判读整幅影像的方法能够显著减少工作量，提高作业效率。由于无法得知实际变更调查中研究区域数据的实际作业时间，因此本文设计了对比实验衡量本文方法效率，由相同熟练程度的操作人员按照年度土地利用变更技术规程要求目视解译提取研究区域变换图斑，将其效率与本文方法对比。对于杜尔伯特蒙古族自治县研究区，本文方法所需时间约为23 min，而目视解译作业时间约为50 min；对于长沙县研究区，本文方法处理所需时间约为70 min，而目视解译作业时间约为2.5 h。与目视判别相比，本文方法所需时间缩短了一半。