土地利用/土地覆盖变化(LUCC)一直被认为是全球环境变化与可持续发展研究的重要内容^[1–3], 在全球气候变化、食品安全、土地退化和生物多样性等关键问题研究中发挥重要作用^[4–6]. 相比以前依靠人力勘察为主的土地利用调查. 目前, 遥感技术是被公认能最快速、最有效地获取土地利用信息及其空间分布的手段^[7], 它具有更新周期快、可研究范围广、客观、技术发展成熟等优点. 近年来, 土地覆盖变化的研究成果多集中在土地利用模型模拟与预测^[8–10]、土地利用覆盖变化及驱动机制^[11–13]、土地利用景观格局变化^[14,15]等方面. 王书玉等^[16]基于随机森林的洪河湿地遥感影像分类研究, 结果表明在湿地信息提取中精度较高. 马玥等^[17]利用随机森林算法对农耕区土地利用分类研究, 较支持向量机和最大似然法总体分类精度分别高9.46%和5.27%. 马骊^[18]对随机森林算法优化改进提高了随机森林的分类性能.

景观格局分析是景观生态学发展的驱动力^[19]. 目前, 研究者主要从两个方面来进行景观格局分析: 模糊的景观格局分析和针对特定生态过程的景观格局分析^[19]. 模糊的景观格局分析是指通过计算各种景观指数, 静态的描述区域景观类型的组成和结构特征, 对于静态的空间格局分析已经做过大量研究^[20,21]. 针对特定生态过程的景观格局分析是在景观生态学家充分认识到格局-过程关系的重要性以后被提出来的, 是景观格局分析的更高层次, 它将在当前及未来的景观格局分析中占据主导地位^[22,23]. 土地利用变化是地表最显著的生态过程, 它受到人为和自然因素的共同影响, 因此, 研究LUCC和景观格局之间的关系对区域景观格局的构建具有重要意义.

2009年福建省根据国务院《关于支持福建省加快建设海峡西岸经济区的若干意见》精神, 决定设立平潭综合实验区. 2011年, “加快平潭综合实验区开放开发”写入国家“十二五”规划纲要和国务院批准的《海峡西岸经济区发展规划》, 平潭开放开发上升为国家战略. 一系列优惠政策促使海坛岛的建设与发展进入了飞速发展阶段, 本文借助遥感和GIS技术, 以海坛岛为研究对象, 定量分析1990~2017年土地利用变化对景观格局的影响. 以揭示近27年来该地区景观生态格局同生态过程和人为活动之间的关系, 从而为该地区土地资源的合理利用及可持续发展奠定基础.

平潭是福建省东部沿海的一个岛县, 位于福州市南部, 与台湾隔海相望. 全县陆地面积392平方公里, 海域面积6064平方公里, 平潭县主岛海坛岛, 为全省第一大岛, 中国第五大岛, 面积为267.13平方公里, 占平潭县总面积的72%, 平潭县东临台湾海峡, 是大陆与台湾本岛距离最近处, 平潭县地理坐标为东经119°32′~120°10′、北纬25°16′~25°44′. 年降水量1196.2毫米, 年均温19.5 ℃. 平潭县境内地势低平, 中部略高, 地形以海积平原为主, 南北有山丘, 海岸蜿蜒曲折. 有天然淡水湖三十六脚湖, 是福建省最大的天然淡水湖(图1).

2009年以前海坛岛土地利用变化并不显著, 但2009年设立平潭综合实验区后, 大规模的经济投入致使土地利用变化显著. 本文选择时间序列长且连续的Landsat数据作为研究的主要遥感数据源, 时间跨度为1990~2017年, 共分4期, 影像均下载自地理空间数据云平台( http://www.gscloud.cn/). 为保证影像质量, 本文仅选取云量小于20%的影像. 其中1990年为Landsat4影像、2000年和2010年为Landsat5影像、2017年为Landsat8影像, 影像下载时已经做过几何校正和地形校正, 后期进行了定标、大气校正和裁剪等预处理, 辅助数据包括海坛岛DEM数据和边界矢量数据. 海坛岛位于湿润亚热带地区, 地物(特别是影响分类的植被)四季变化不明显, 可选择不同季节的影像进行土地利用信息的提取, 研究所用软件有ArcGIS 10.3、FRAGSTATS 4.3和ENVI 5.3等. 具体数据源信息(表1).

图 1

Fig. 1

图 1 研究区域

表 1(Table 1)

表 1 数据源信息 成像时间 传感器类型 条带号 中新经纬度 多光谱分辨率(m) 1990年7月5日 Landsat4 MSS 118/42 E26°014′, N120°386′ 30 2000年7月24日 Landsat5 TM 118/42 E26°000′, N120°429′ 30 2010年7月20日 Landsat5 TM 118/42 E26°001′, N120°445′ 30 2017年10月27日 Landsat8 OLI 118/42 E25°993′, N120°435′ 30

表 1 数据源信息

随机森林是2001年被提出的一种机器学习方法, 在原始训练样本中通过bootstrapping方法随机有放回抽取n个样本, 构成训练样本集, 每个训练样本集可构造一个决策树模型, 这样可生成n颗决策树模型. 在每颗决策树的分裂节点处随机选取m个特征变量作为预测变量, 根据节点不纯度最小原则来选取最优的特征变量作为该分裂节点的分裂特征, 直到每个分裂节点不纯度最小时决策树停止生长^[24]. 生成所有决策树后, 用投票的方法对所有决策树的分类结果进行综合, 得出最终结果^[25,26]. 森林中每一棵树都依赖于一个随机向量, 森林中的所有向量都是独立分布的. 随机森林的预测精度与单个树的强度和树间的相关性有关^[27].

应用随机森林模型分类器需要设定2个关键参数: ① 生成树个数(number-of-trees), 该值决定了随机森林总体规模. 取值一般在100~10 000之间, 值越大模型越收敛, 模型运行时间增加, 且当树的个数不断增加时, 模型会出现过饱和现象. ② 随机抽取最大特征数(max-features), 该值表示生成每一棵决策树时从特征空间中随机抽取特征的最大数. 该值越大模型中每棵决策树强度越大, 但决策树间相关性也增大. 因此, max-features需根据oob误差率来调优, 以达到一定精度^[28].

本文采用2017年发布并实施的土地利用现状分类国家标准, 选取4个一级类(林地、耕地、水域和建设用地), 2个二级类(沙地和裸土地). 通过ENVI 5.3插件随机森林分类器可实现分类.

随机森林分类器的优点为^[27]; ① 算法精度高; ② 可以处理大数据集, 无需对大量的输入特征变量进行删减; ③ 可以给出变量的重要性估计; ④ 在模型建立过程中可以产生一个对一般误差的无偏估计, 不会过度估计; ⑤ 可以有效处理缺失数据的情况; ⑥ 产生的森林模型易于保存和未来重复利用; ⑦ 可以扩展到无类别数据中, 进行非监督分类.

在排除高相关性指标的基础上, 研究在类型级别上选用了4个指标, 分别为最大斑块所占景观面积的比例、边界密度、聚合度和周长-面积分维数. 在景观级别上选用了4个指标蔓延度、香农均匀度、香农多样性和景观形状. 运用FRAGSTATS 4.3软件, 采用eight-cell-neighborhood规则计算景观指标, 指标的描述如表2所示.

基于随机森林理论, 选取334个训练样本(1990年85个、2000年79个、2010年8个、2017年87个)将地物分为林地、耕地、水域、建设用地、沙地和裸土地6类, 选择分类树数目500对研究区四期的Landsat影像进行分类, 分类结果如图2所示.

表 2(Table 2)

表 2 景观格局指标 景观类型 指标 指标缩写 土地利用格局意义 斑块类型 最大斑块所占景观面积的比例 LPI 反映景观格局的规模或景观中的优势斑块类型等, 间接反映出土地利用变化的方向和强弱 边缘密度 ED 反映各种土地利用类型的基本特征 聚合度 AI 反映土地利用格局中同类型像元的聚集程度 周长-面积分维数 PAFRAC 反映不同景观类型的性状的复杂性 景观格局 蔓延度 CONTAG 反映区域土地利用格局中不同土地利用类型的聚集程度 均匀度 SHEI 反映各种土地利用类型分配的均匀程度, 其值越大, 表明各土地利用类型分配越均匀 多样性 SHDI 反映区域中土地利用/景观类型的多少和各类型所占比例的变化 景观形状 LSI 反映区域中土地利用类型/景观类型整体几何复杂程度

表 2 景观格局指标

图 2

Fig. 2

图 2 影像分类结果

表3统计各地类面积占比, 1990年和2000年, 海坛岛土地利用类型面积从大到小排序为耕地>林地>建设用地>水域>裸土地>沙地. 2010年, 海坛岛土地利用类型面积从大到小排序为耕地>林地>建设用地>水域>沙地>裸土地. 2017年, 海坛岛土地利用类型面积从大到小排序为耕地>林地>沙地>水域>裸土地. 以2017年为例, 耕地占40.00%, 林地占 24.96%, 建设用地占28.43%, 三者之和占92.00%以上, 其余土地利用类型面积占比不到8.00%, 因此研究区耕地、建设用地和林地为主要景观. 1990~2017年建设用地面积持续增加14.12%; 水域和裸土地面积分别减少了10.85%、2.68%; 林地面积先减少后增加, 最大变化量为7.29%; 沙地面积先增加后减少, 最大变化量为4.56%.

表 3(Table 3)

表 3 1990~2017年土地利用类型面积占比统计表(%) 年份 林地 水域 建设用地 耕地 沙地 裸土地 1990 27.23 13.16 14.11 40.99 1.38 3.13 2000 23.63 12.37 17.38 42.44 1.07 3.11 2010 19.94 9.58 18.96 42.92 5.63 2.96 2017 24.96 2.31 28.43 40.00 3.85 0.45

表 3 1990~2017年土地利用类型面积占比统计表(%)

分类后利用计算机按土地利用类型随机生成若干个检验点(1990年256个、2000年263个、2010年274个、2017年264个). 通过目视解译的方法对检验点进行分析判读, 采用误差矩阵和Kappa系数对随机森林分类结果进行精度评价计算出的误差矩阵(表4). 从4期分类误差矩阵看, 制图精度中水域最高, 沙地最低; 用户精度中水域最高, 裸土地最低, 整体来看, 4期总分类精度在87.07%以上, 总Kappa系数在0.84以上.

表 4(Table 4)

表 4 1990~2017随机森林分类器的误差矩阵计算结果 年份 类别 林地 水域 建设用地 耕地 沙地 裸土地 行和 制图精度(%) 用户精度(%) Kappa系数 1990 林地 36 1 1 2 0 0 40 90.00 90.00 0.8815 水域 0 40 0 1 0 0 41 95.24 97.56 0.9708 建设用地 1 0 50 1 1 2 55 90.91 90.91 0.8842 耕地 2 0 1 45 2 1 51 86.54 88.24 0.8524 沙地 0 1 1 1 24 3 30 77.42 80.00 0.7724 裸土地 1 0 2 2 4 30 39 83.33 76.92 0.7315 列和 40 42 55 52 31 36 256 总分类精度=87.89 总Kappa=0.8536 2000 林地 38 1 2 1 0 0 42 88.37 90.48 0.8861 水域 1 45 0 0 0 0 46 97.83 97.83 0.9737 建设用地 1 0 47 2 1 3 54 90.38 87.04 0.8384 耕地 3 0 1 40 3 1 48 83.33 83.33 0.7961 沙地 0 0 1 2 27 3 33 77.14 81.82 0.7903 裸土地 0 0 1 3 4 32 40 82.05 80.00 0.7652 列和 43 46 52 48 35 39 263 总分类精度=87.07 总Kappa=0.8443 2010 林地 40 1 0 2 0 0 43 88.89 93.02 0.9165 水域 0 48 0 1 0 0 49 97.96 97.96 0.9751 建设用地 1 0 45 1 1 2 50 90.00 90.00 0.8777 耕地 4 0 2 43 5 1 55 81.13 78.18 0.7295 沙地 0 0 1 2 28 2 33 75.68 84.85 0.8248 裸土地 0 0 2 4 3 35 44 87.50 79.55 0.7605 列和 45 49 50 53 37 40 274 总分类精度=87.23 总Kappa=0.8462 2017 林地 53 0 1 2 0 1 57 91.38 92.98 0.9101 水域 1 41 0 0 1 0 43 95.35 95.35 0.9444 建设用地 1 0 47 2 2 2 54 90.38 87.04 0.8386 耕地 2 0 1 40 2 2 47 86.96 85.11 0.8196 沙地 1 1 2 1 20 1 26 74.07 76.92 0.7429 裸土地 0 1 1 1 2 32 37 84.21 86.49 0.8421 列和 58 43 52 46 27 38 264 总分类精度=88.26 总Kappa=0.8576

表 4 1990~2017随机森林分类器的误差矩阵计算结果

每次计算以前一时期的土地利用分类结果为标准, 求出4个时期土地利用类型间的相互转化量, 得到土地利用转移矩阵(表5). 1990~2000年间, 建设用地面积增加9.36 km², 主要源于耕地和林地的减少; 林地面积减少10.29 km², 流失的林地主要转化为耕地和建设用地; 水域面积减少2.25 km², 流失的水域面积主要转化为建设用地; 耕地面积增加4.15 km², 主要源于林地的减少; 沙地面积减少0.90 km², 流失的沙地主要转化为耕地和裸土地; 裸土地面积减少0.07 km², 流失的裸土地主要转化为建设用地和耕地. 1990~2000年间土地利用变化不明显, 主要是当时经济发展速度有限, 大部分土地都保持原有面貌. 但随着人口的缓慢增加, 建筑用地和耕地也会略有增加. 2000~2010年间, 建设用地面积增加4.53 km², 主要源于耕地和林地的减少; 林地面积减少10.55 km², 流失的林地主要转化为耕地和建设用地; 水域面积减少7.98 km², 流失的水域面积主要转化为建设用地和沙地; 耕地面积增加1.38 km², 主要源于林地的减少; 沙地面积增加13.04 km², 增加的沙地主要源于水域和林地减少; 裸土地面积减少0.41 km², 流失的裸土主要转化为建设用地和沙地. 2000~2010年间土地利用变化较为缓慢, 其主要原因与前一时期类似.

表 5(Table 5)

表 5 1990~2017土地利用转移矩阵(km2) 2000年 1990年 林地 水域 建设用地 耕地 沙地 裸土地 林地 50.33 0.82 7.94 6.64 0.45 1.43 水域 2.30 31.98 0.72 0.28 0.02 0.09 建设用地 8.07 3.78 22.99 11.85 0.60 2.43 耕地 16.40 0.51 6.07 95.58 1.11 1.75 沙地 0.11 0.17 0.36 0.60 1.02 0.80 裸土地 0.69 0.39 2.30 2.32 0.75 2.45 2010年 2000年 林地 水域 建设用地 耕地 沙地 裸土地 林地 44.44 0.66 3.36 8.45 0.04 0.10 水域 0.98 23.40 1.81 0.94 0.12 0.16 建设用地 4.25 3.33 28.78 13.59 0.61 3.71 耕地 15.73 0.66 12.31 92.49 0.43 1.18 沙地 1.20 7.30 2.12 1.28 1.69 2.51 裸土地 1.02 0.05 1.35 4.68 0.16 1.23 2017年 2010年 林地 水域 建设用地 耕地 沙地 裸土地 林地 44.82 3.13 3.00 18.15 1.68 0.63 水域 0.08 5.63 0.34 0.10 0.31 0.14 建设用地 4.32 13.94 36.16 17.15 7.57 2.20 耕地 7.21 2.10 12.99 84.66 2.47 5.00 沙地 0.50 2.49 1.66 2.01 4.05 0.29 裸土地 0.11 0.11 0.10 0.73 0.01 0.22 2017年 1990年 林地 水域 建设用地 耕地 沙地 裸土地 林地 44.60 5.25 6.31 13.25 0.56 1.44 水域 0.74 4.99 0.43 0.33 0.02 0.09 建设用地 12.25 18.14 20.30 25.81 1.28 3.56 耕地 19.13 3.85 12.23 74.88 1.33 3.02 沙地 0.94 5.32 0.97 2.28 0.76 0.75 裸土地 0.24 0.09 0.12 0.71 0.01 0.10

表 5 1990~2017土地利用转移矩阵(km2)

2010~2017年间, 建设用地面积增加27.09 km², 主要源于水域、沙地和耕地的减少; 林地面积增加14.35 km², 增加的林地主要源于水域和耕地的转化; 水域面积减少20.80 km², 流失的水域面积主要转化为建设用地和林地; 耕地面积减少8.36 km², 流失的耕地主要转化为建设用地和林地; 沙地面积减少5.08 km², 减少的沙地主要转化为建设用地和林地; 裸土地面积减少7.20 km², 流失的裸土主要转化为建设用地.

1990~2017年海坛岛土地利用类型与景观类型级别指标关系如图3所示. 从图3可见, 1990~2000年, 耕地的最大斑块占景观面积比例(LPI)和聚合度(AI)呈上升趋势. 因为耕地在这一时期是最主要的景观类型, 分布规模最大且面积增加, 耕地分布更加聚集. 聚集的耕地使边界周长减小, 在面积相同的情况下, 导致耕地的边缘密度(ED)和周长-面积分形维数(PAFRAC)下降. 林地的LPI、ED和PAFRAC呈下降趋势, 因为林地面积减少, 林地的规模在缩减, 且减少的林地位于林地的集中分布区, 研究区的东北部(图2). 建设用地LPI、ED、PAFRAC和AI均呈上升趋势, 主要是建筑设用地面积增加, 且在上一期的建设用基础上向外扩散造成的.

2000~2010年, 耕地的最大斑块占景观面积比例(LPI)和聚合度(AI)呈下降趋势. 因为耕地在这一时期面积增加缓慢, 且耕地最大斑块面积减小, 由图2可知, 耕地分布聚集度下降. 聚集度下降的耕地使边界周长增加, 则在面积相同的情况下, 导致耕地的边缘密度(ED)和周长-面积分形维数(PAFRAC)上升. 林地的LPI和PAFRAC呈上升趋势, 因为林地最大斑块面积增加. 随着林地面积的减少, ED和AI均下降.

2010~2017年, 建设用地的最大斑块占景观面积比例(LPI)、聚合度(AI)和边缘密度(ED)呈上升趋势. 因为建设用地在这一时期面积显著增加, 分布规模增大. 边缘长度增加明显大于面积增大. 但建设用地整体形状复杂性下降, 即周长-面积分形维数(PAFRAC)下降. 水域的LPI、ED、AI和PAFRAC均呈下降趋势, 因为水域面积急剧减少, 水域形状由复杂变简单且分布零散, 聚集度下降. 林地的面积有所增加, 人工植树造林导致分布集中形状复杂. 因此林地的LPI、ED、PAFRAC和AI均呈上升趋势.

图 3

Fig. 3

图 3 1990~2017类型级别指标与土地利用关系图

表6统计景观级别指标变化, 1990~2000年, 蔓延度(CONTAG)上升0.713, 说明这10年间不同土地利用类型之间的聚集程度增加. 景观形状(LSI)下降3.315, 说明景观的几何复杂性在减小, 同时表明相同土地利用类型之间分布更加离散. 1990~2010年, 多样性(SHDI)和均匀度(SHEI)变化不明显, 说明3期的土地利用类型分配的均匀程度相似, 各土地利用类型所占比例变化不明显, 1990~2010年间的土地利用变化也可以看出. 2000~2010年, CONTAG下降1.936, 说明这10年间不同土地利用类型之间的聚集度程度降低. LSI上升7.812, 说明景观的复杂性在增加, 同时表明相同土地利用类型之间分布更加集聚. 2010~2017年, CONTAG上升2.373, 说明近10年间不同土地利用类型之间的聚集程度增加较1990~2010年显著, LSI下降0.783, 说明景观的几何复杂性在减小, 同时表明相同土地利用类型之间分布更加离散. SHDI和SHEI分别增加了0.081和0.042, 说明不同土地利用类型分配的均匀程度增加, 主要土地利用类型占比增加, 2010~2017年土地利用急剧变化, 建设用地增加, 水域减少是导致SHDI和SHEI变化的主要原因.

表 6(Table 6)

表 6 1990~2017年景观级别指标统计表 年份 CONTAG SHDI SHEI LSI 1990 57.6567 1.3191 0.6779 38.8408 2000 58.3701 1.3160 0.6763 35.5256 2010 56.4344 1.3457 0.6915 43.3376 2017 58.8907 1.2644 0.6498 42.5551

表 6 1990~2017年景观级别指标统计表

运用随机森林法进行遥感影像分类, 分类精度很大程度上受制于训练样本的选择. 本文通过目视解译选取较纯净的地类像元作为训练样本, 经过多次实验获得较高的土地利用分类精度(4期遥感影像分类的总体精度均在87％以上, Kappa系数均在0.84以上).

2010年以前, 研究区土地利用覆盖整体上变化不明显, 但在2010~2017年间土地利用变化显著, 主要是2009年设立平潭综合试验区后, 大量的资金注入以及政策支持造成, 海坛岛开发进入快车道, 为了获取足够的发展空间, 岛内进行了大规模的围垦造陆, 使得近海区域的滩涂、湿地和水域等景观被大规模开发利用, 其中大部分被规划为建设用地使用. 研究显示, 2010~2017年耕地面积变化很小, 说明国家对耕地的保护力度很大. 林地面积在2010~2017年略微增加, 说明在开发的同时考虑到生态环境的重要性.

1990~2017年, 持续增加的建设用地致使斑块类型级别的最大斑块占景观面积比例(LPI)、聚合度(AI)和边缘密度(ED)景观指数呈上升趋势, 其中LPI受建设用地增加的影响最为显著. 景观多样性(SHDI)和景观形状(LSI)呈下降趋势, 表明研究区整体景观破碎化程度和景观异质性降低. 因此, 海坛岛的发展应该注重景观格局的演变效应, 提出适合海岛型城市发展的模式, 从而实现城市发展与景观格局优化的共赢.