基于像元转换的土地覆盖变化监测方法——以北京市区县为例

于 冰, 王继燕, 苏 勇, 肖东升, 刘福臻

(1.西南石油大学土木工程与建筑学院,成都 610500; 2.东华理工大学江西省数字国土重点实验室,南昌 330013)

0 引言

土地覆盖是一个复杂地表要素覆盖的自然属性或人工建造的综合体。检测土地覆盖变化是对同一空间区域的地表覆盖对象、不同时间观测的空间现象来确定变化检测的处理过程，依据不同时期的地理信息数据或者遥感图像来提取和描述随时间产生的对象变化[1-2]。遥感影像因其获取时间短、覆盖范围广等优势而被广泛应用于土地覆盖变化检测。覆盖变化检测方法通常可分为2大类： 分类前变化检测法和分类后比较法[3]。

分类前的变化检测，即影像对影像的比较[4]。该类方法包括分析不同影像的光谱特征来识别感兴趣变量引起的真实光谱差异[5]。相关方法中有的仅得到了变化和未变化的结果，而有的则可得到土地覆盖变化检测矩阵[6-8]。如支持向量机的变化检测[9]，对象变化矢量分析[10]、基于图像像斑统计分析的高空间分辨率影像变化检测[11]和基于像斑光谱向量相似度的变化检测[12]等。对于影像的变化检测也有一些通用软件包，如ENVI，Erdas Imagine和Fragstats等。早在20世纪末，EI-Gamily[13]使用Erdas Imagine 处理软件中的工具来检测土地覆盖变化，结果显示了Kuwait东南部地区环境破坏的变化情况。

分类后的变化检测，即图与图的比较，主要基于空间叠置分析和统计原理分析。该类方法简单易操作，但是方法的有效性以及检测结果的可靠性则依赖于遥感影像的分类精度。李静等[14]早期使用一种叠合像元二值化的变化检测方法进行变化检测，该方式仅能检测变化及未变化区域； 张学等[15]提出一种扩展的基于像元转化的变化检测方法，该方法给出了各种土地覆盖类别相互转化的转换公式，提高了分类后数据变化检测的效率； 唐朴谦等[16]提出了针对遥感分类后数据的一种基于像素比值的面向对象检测方法，该方法可以有效地减少误差传递的影响，尤其是对未发生变化区域分类所产生的误差。尽管如此，分类后变化检测仍然有可探究的新方法和新策略。

本文基于像元转换的思想，提出一种新的像元转换的土地覆盖变化检测方法，旨在优化检测过程中的处理策略。以北京市辖区2000—2010年间土地覆盖数据为例，计算各种土地覆盖类别的变化情况，同时探讨研究区样点土地覆盖变化的驱动力，验证该方法的适用性。

1 变化检测方法原理

利用GlobeLand30土地覆盖图对北京市16个区2000—2010年土地利用/土地覆盖变化进行了统计分析。数据可从http: //globeland30.org/GLC30Download/index.aspx网站免费下载，下载后数据经过拼接、投影变换和裁剪等预处理后得到北京市2期土地覆盖图。原数据类别编码如表1所示，需要进行像元编码转化； 然后，进行叠置分析、分区统计分析等得到北京市辖区土地覆盖变化情况信息； 最后，使用土地覆盖变化量化评价指标和景观测度指标对检测变化结果进行分析。

表1 覆盖类别新旧对照Tab.1 Old and new value for land categories

1.1 转换像元的变化检测方法

先前研究中的土地覆盖变化检测方法，主要是通过叠置分析得到变化区域与未变化区域，即先将分类结果灰度图进行二值化，然后进行变化检测计算，过程存在较大不确定性。由于二值化的阈值是由人为设定，其结果存在较大误差。为克服这一缺陷，本文提出一种基于像元转换思想的变化检测方法，该方法不需要对分类后结果进行二值化，而是在变化检测的过程中将分类结果中的像元值重新编码，然后进行叠置分析达到变化检测的目的。这种方法可以用于多种土地覆盖类别间的变化检测，增强了方法的适用性，简化了处理流程。

原始数据中不同栅格编码值对应不同类别，2期数据同一覆盖类别的栅格值做栅格计算将无法得到类别变化信息。土地覆盖变化结果是由2个不同时相的土地覆盖图叠置运算得到的。将原先不能进行变化叠置计算的图像像元值统一重分类到唯一的像元灰度值编码，重新定义各种类别像元编码值。例如表1中，类别变化完全不同的情况下有可能得到一样的计算差值，例如耕地变为森林和苔原变为人造覆盖的栅格计算后差值都是10，但是该值无法唯一表示从哪类地物变为另外一类覆盖地物。因此，需要对原编码值进行编码，然后计算变化结果。例如重新定义表1中类别新编码值依次为2的N次方(N=0，1，2，…，9)，新编码值可用于栅格作差计算，计算后得到的值可以唯一标识类别变化状态。该唯一标识值表示从A类型变为B类型，8种类型之间的标识一共有57个像元编码值，其中为0的值表示前后没有发生变化。在处理的过程中，可以利用ArcGIS model builder综合数据处理工具，简化土地覆盖类别变化分析流程。变化检测分析流程如图1所示。

图1 数据处理流程Fig.1 Flowchart of data processing

数据关键处理步骤如下：

1)重分类。原始数据中不同栅格编码值对应不同类别，但对2期土地覆盖类别的栅格值做栅格计算将无法得到类别变化情况。对原编码值进行重新编码,可用于重分类(reclassify)，改变像元值，定义表1中类别新编码值依次为2的N次方(N=0，1，2，…，9)。

2)定义对照表。新的编码可以预先定义类别转换状态，以便于后续数据表关联作变化统计。栅格计算器对2期数据求差后得到唯一栅格值指示类别转移状态，其转化对照表如表2所示，第一列表示2000年原始土地覆盖类别状态，第一行代表2010年变化后土地覆盖类别状态，表格中数值128代表人造覆盖变化为荒地，-128代表荒地变为人造覆盖。栅格作差数据与对照表进行字段连接，得到像元点指代的变化情况。

表2 类别变化唯一值对照表Tab.2 Unique value for cover change

3)制表求交分析。通过上述计算得到北京市整体地表覆盖变化分布数据，结合北京市辖区行政矢量图使用制表求交工具计算辖区的变化前后情况。该求交计算可以快速获得行政区域内类别统计分布情况，如类别数、类别面积及比例[17]。

通过上述变化计算流程，得到了变化量及变化空间状况。随后利用土地利用/土地覆盖变化评价指标和景观指标来揭示变化结果，可简要分析变化驱动原因。

1.2 变化评价指标

1.2.1 土地利用/土地覆盖动态变化分析

单一土地利用/土地覆盖动态度指标用于表示某一研究地区一定时间跨度范围内某一类覆盖类别的数量变化情况，通常用该比率值表示变化程度[18]，其表达式为

(1)

式中：LC为单一土地利用动态度；Ua和Ub分别为研究期始末单一土地覆盖类别的数量(如面积大小)；T为研究时间跨度，当T的单位为a时，LC的值代表该研究区某一种土地利用/土地覆盖类别年变化程度。

1.2.2 景观指标

景观指标是反映一个场景要素特征的重要参数。其中斑块数量和大小直接反映出景观格局的聚集度和丰富性等特征[19-20]。通过不同时间段斑块数量及斑块大小的对比分析，可以探寻斑块分布对土地覆盖变化的影响，进而分析变化趋势，预测土地覆盖变化发展方向。斑块描述指标包括斑块个数(NP)、斑块总面积(AT)、斑块面积标准差(SD)、平均斑块面积(AV)、斑块面积变异系数(CV)、最大斑块指数(LPI)、斑块丰富度(PR)等[21]。在此使用Fragstats工具(http: //www.umass.edu)计算NP，LPI，PR和聚集指数(AI)。

1)NP代表整个景观场景中斑块的总个数。斑块是景观格局的基本组成单元，是指不同于周围背景的、相对均质的非线性区域。斑块统计需要遵循邻域规则，规则分为四邻域和八邻域2种。

2)LPI表征某一类型最大斑块在整个景观中所占比例。计算公式为

(2)

式中：ai为景观中第i个斑块的面积；A为整个景观的总面积。

3)PR为景观中不同类型斑块的数量和。

4)AI则基于同类型斑块像元之间公共边界长度来计算。当某类型中所有像元间不存在公共边界时，该类型的聚合程度最低； 当类型中所有像元间存在的公共边界达到最大值时，具有最高聚合程度。计算公式为

(3)

式中gij为景观类型的相似邻接斑块数量。

2 数据处理与结果分析

2.1 研究区概况与数据

本实验中土地覆盖数据来源于全球30 m地表覆盖数据(GlobeLand30)。该数据是目前世界上首套30 m空间分辨率的全球地表覆盖产品，精细程度高出同类数据产品10倍。数据包括10个一级大类，包括湿地、水体、人造覆盖、灌木丛、耕地、森林、草地、荒地、苔原、永久积雪和冰川。这一数据共有2000年和2010年2期，数据分类精度在80%以上[22]。该数据已将全球地表覆盖数据精度从原来的百m级提高到了30 m，不仅可以满足国土资源管理、环境、农业和城市化等科学研究需求，还可为全球变化研究、地球系统模拟和可持续发展研究提供数据支撑[23]。

本文实验选取北京市为土地覆盖变化分析区域。按照北京市行政区划范围，现有16个市辖区。2期土地覆盖数据包含类别一致，共有水体、湿地、人造覆盖、耕地、森林、灌木丛、草地、荒地8个地类，不含有苔原、永久积雪和冰川。研究区概况如图2所示，不同土地覆盖类别面积统计如表3所示。

(a) 2000年 (b) 2010年

图2 研究区范围Fig.2 Study area表3 北京市2期数据土地覆盖数据Tab.3 Two phases land cover data of Beijing City (km2)

2.2 北京市土地覆盖变化

土地覆盖变化是监测地表覆盖状况的重要手段，土地覆盖变化的类型可分为3种： 土地覆盖退化、土地覆盖转换和土地覆盖改良。土地覆盖退化是指某种土地覆盖种类虽然没有发生变化，但是其质量已经降低，例如过度放牧或者垦荒导致草地退化，过度砍伐引起森林植被密度降低； 土地覆盖转换是指某种覆盖类别完全改变成了另一种类别，例如农业耕地被开发建设成住宅、商业用地，耕地变为人造覆盖用地，林地被开垦为耕地，或者建设为旅游商业用地； 土地覆盖改良是指某种覆盖类型获得保护、修复、更新，譬如森林抚育、土壤/草地改良等。这3种土地覆盖变化中，覆盖转换凸显变化最为明显，比较容易引起人们的关注，也比较容易监测，而退化和改良因其记录不完备不易观测。

北京市GlobeLand30数据共包含有8个类别，利用上述变化分析方法对土地覆盖变化做了统计描述。统计分为北京市整体变化与区级详细变化分布，同时给出了变化区域的空间位置。根据表3数据和式(1)计算单一土地利用动态度，北京市单一土地利用动态度如图3所示。

图3 北京市单一土地利用动态度Fig.3 Dynamic degree of single land use of Beijing City

由图3可知，人造地表覆盖面积10 a间增长超过50%，直接说明北京城镇化发展迅猛。其中耕地、草地面积较之前有明显减少，这些地类从农业用地转向为建筑用地(人造覆盖)； 水体面积也明显减少，可能退化或转向为农业耕地。土地利用变化主要表现为农业用地向建设用地扩展，空间格局上表现为外城蔓延和郊区城市化等增长形式[24]。

为进一步探究类别间变化情况，量化计算了各类别状态转移矩阵(表4)，并绘制了类别变化堆积图(图4)。

表4 北京市类别变化转移矩阵Tab.4 Category change transition matrix of Beijing City (km2)

①0.00是指该数值小于0.005 km2。

图4 北京市10 a间土地覆盖类别转移变化Fig.4 Category change of Beijing City for 10 years

从图4可以看出，原先是草地、耕地和森林覆盖类别转换为其他类别面积值最大，其中2000年是草地覆盖的地区到2010年大量变为森林、耕地和人造覆盖； 2000年是耕地覆盖的区域到2010年大量变为人造覆盖、森林和草地类型； 森林类型也主要变为草地、人造覆盖和耕地。原先是人造覆盖的区域也有不少变为耕地。由于草地、森林和人造覆盖转变为耕地的面积远小于耕地转变为其他类别的面积，所以总体上耕地面积是减少的，同理其他类别也因为类别相互变化出现类似情况。

为验证北京市人造覆盖区域的变化发展情况，绘制出2个时期的人造覆盖区域的空间分布变化(图5)。从图5可以看出，城市发展依附于城市本身形态以及向外形成环状或块状城市地区，外围土地利用和开发不断向周围地区蔓延，市区向南北和东方向扩张量较大，表现出明显的轴向扩张。从城市扩张的土地类别来源来看，主要是耕地、草地和森林这些农业用地为主的郊区转变为以非农业用地为主的城市地区。从北京市的郊区发展看，可以推测到由于城市人造覆盖区域面积的扩大，市郊的农业用地以及一些自然草地和水体区域逐渐转变为工厂、商店及住宅等非农业用地。该结论也验证了北京市城市化增长形式： 轴向扩展和郊区城市化。

图5 北京市人造覆盖区域扩张变化分布Fig.5 Artificial land change map of Beijing City

利用上述方法，计算了2期地表覆盖数据的4个景观指数(表5)，NP值减少了2 553个，LPI值有所增加，AI值增大，同时说明北京市土地覆盖类别更加成片化，以前零碎的斑块被整合，这也是用地向集中节约型发展的一个重要线索。PR值没有太大的变化，说明覆盖类别总体比较稳定。

表5 2期数据景观指数Tab.5 Landscape metrics for two phases data

图6显示了变化区域变化前后的覆盖类别。根据土地覆盖变化数量(地类面积增减量、变化率)覆盖变化方向等因素，结合土地利用动态度，可分析出各个区土地利用/土地覆盖变化的强度、速度和方向。其中，门头沟区、大兴区和通州区3个区域人造覆盖区域变化超过了原先面积的50%，昌平区人造覆盖区域变化超过49%，具体数据参见表6(表中用加粗标记了人造覆盖变化在50%左右的数据)。进一步验证城市扩展主要向外围发展。耕地出现锐减的还有西城区和丰台区，其中西城区原有耕地变为林地。朝阳区耕地和草地有大量的减少，主要转变为人造覆盖和林地，说明该区域建筑用地和绿化面积有很大增长。

(a) 变化前 (b) 变化后

图6 各辖区2期数据土地覆盖类别空间分布Fig.6 Category change of spatial distribution during two periods表6 北京市各辖区类别变化统计Tab.6 Land cover change of zone and county in Beijing City (%)

3 结论

针对同一地区不同时相获得的地表覆盖数据的变化检测方法简单，但是需要通过多个步骤计算获得覆盖变化量、变化方向和格局变化特征。因此，提出一种基于像元转换的变化检测方法。此方法主要是依据土地覆盖数据中不同土地覆盖类别具有不同的编码值的特点，通过重新编辑编码值叠合计算得到可以唯一标识变化情况的值，进而应用分区统计分析重叠计算得到土地覆盖变化图像。

选取北京市的Globeland30数据，利用该方法和处理策略探究了北京市不同市辖区2000—2010年间土地覆盖类别变化情况。土地覆盖变化显著，尤其是人造覆盖区域增长到10 a前的1.5倍，耕地和草地面积减少明显，同时森林面积增加，这说明城市生态及周边变化明显。通过变化检测结果更为精细地展示了覆盖类别变化和景观格局变化。

本文提出的方法利用现有的软件平台集成处理不同时相数据，获得了变化检测的结果。在变化检测处理的过程中，提出一种简化的计算分析策略，能够快速地分析出土地覆盖类别变化情况，包括量化转移和空间位置的变化。该分析方法对土地利用/土地覆盖变化分析具有较强的适应性。

志谢： 此次实验的数据获取得到了国家基础地理信息中心的支持，在此表示衷心的感谢！