EN中文

Journal of Geo-information Science >

2022 , Vol. 24 >Issue 2: 310 - 325

DOI: https://doi.org/10.12082/dqxxkx.2022.210034

Multi-temporal Cultivated Land Cover Extraction and Change Analysis: A Spatiotemporal Context Method Combining Remote Sensing and Spatial Statistics

  • XI Wenqiang ,
  • DU Shihong , * ,
  • DU Shouji
Expand
  • Institute of Remote Sensing and GIS, Peking University, Beijing 100871, China
*DU Shihong, E-mail:

Received date: 2021-01-20

  Request revised date: 2021-03-10

  Online published: 2022-04-25

Supported by

National Natural Science Foundation of China(41871372)

Copyright

Copyright reserved © 2022
Fold

Abstract

The extraction and change analysis of cultivated land covers based on multi-temporal images are important means to effectively manage and protect cultivated land resources. However, as far as the classification and extraction of multi-temporal cultivated land covers are concerned, the existing methods have limitations in the comprehensive expression of spatiotemporal features and the accurate modeling of spatiotemporal context relations among geographic objects, leading to the poor extraction accuracy of cultivated land covers. In addition, for the analysis of cultivated land change, the existing methods usually only focus on the statistical areal change of cultivated land cover based on administrative units, while consideration is seldom taken into the spatial correlation distribution characteristics of changes in cultivated land covers. Accordingly, first of all, this paper proposes a multi-temporal spatiotemporal context classification method, which comprehensively expresses and utilizes the multi-temporal spectral, texture, and spatial features of geographic objects, and models the contextual relations of features and semantics among geographic objects in both spatial and temporal dimensions of multi-temporal images, so as to improve the classification accuracy of cultivated land covers. Then, based on the extracted results of cultivated land covers, spatial statistical method of Geographic Information System (GIS) is used to analyze the spatial correlation characteristics of cultivated land changes in regular grids and administrative division units. Finally, Shunyi District of Beijing is taken as the study area while multi-temporal Sentinel-2 images in 2015-2019 are used as the data sources to conduct verification of the proposed method. The results show that, compared with the two existing common multi-temporal classification methods, the proposed method achieves the highest accuracies in the classification of multi-temporal cultivated lands. The average user's accuracy and producer's accuracy reach 91.21% and 90.53%, respectively, while the overall accuracy of all categories is 90.79%, indicating that the proposed method can accurately extract the multi-temporal cultivated land cover information. Furthermore, according to the analysis of the spatial distribution characteristics of cultivated land changes, this study found a phenomenon of regional aggregation of the cultivated land change in Shunyi District from 2015 to 2019, which mainly presents the characteristics of concentrated reduction. The aggregation reduction phenomena of cultivated land covers in Zhaoquanying Town, Gaoliying Town, Mulin Town, and Yang Town are especially obvious, indicating that the problems of cultivated land encroachment and reduction are quite serious in these areas.

Key words： cultivated land; multi-temporal data; spatiotemporal features; spatiotemporal context; conditional random field; image classification; spatial statistics; change analysis

Cite this article

XI Wenqiang , DU Shihong , DU Shouji . Multi-temporal Cultivated Land Cover Extraction and Change Analysis: A Spatiotemporal Context Method Combining Remote Sensing and Spatial Statistics[J]. Journal of Geo-information Science, 2022 , 24(2) : 310 -325 . DOI: 10.12082/dqxxkx.2022.210034

1 引言

耕地作为农业发展和粮食生产的保障,是关系人类生存和社会稳定的重要资源[1,2]。然而,城市化发展使耕地减少和退化问题日益严重[3,4]。因此,科学管理和保护耕地资源,有效获取耕地分布和变化信息尤为重要。因其大范围、重复一致的地表观测能力,遥感已成为监测耕地资源的重要手段[5]。 21世纪以来,遥感传感技术的发展（如MODIS, Landsat, Sentinel）,使传感器能在多个时相上获取高分辨率数据[6],从而能更精确地提取耕地信息以及分析变化,但同时也为多时相数据的有效处理和应用带来挑战。
在基于多时相遥感影像的耕地监测方面已有大量研究,可概括为2个方面：耕地覆盖提取与制图以及耕地覆盖变化分析。耕地覆盖提取主要基于多时相土地覆盖分类,现有研究主要存在2个问题。① 对影像时空特征的综合表达和利用上存在不足。多时相影像包含丰富的时空信息,如地物光谱表征、空间分布及不同时间上的分布信息。这些信息在不同维度上描述了地物的分布特点,从而有助于区分地物[7]。然而,现有方法往往仅利用了多时相光谱等特征而忽略地物的空间分布特征。例如,魏鹏飞等[8]、Slagter等[9]和阴海明[10]等仅使用多时相NDVI、EVI等光谱特征进行地物分类;Peña等[11]比较了多时相光谱、NDVI和NDWI等不同的组合对于多时相分类的影响,但忽略了空间几何和形状特征的表达,因此无法有效消除光谱的空间异质性影响,导致制图精度不佳。② 在多时相影像中,空间和时间上相邻的地物在特征和语义上存在着相关性,这种相关性可用时空上下文信息来度量,对耕地的分类有重要的意义[12]。现有研究一部分是采用传统面向对象方法分别对各时期影像单独分割分类处理,忽略了不同时相影像中上下文的表 达[13,14]。另有研究采用多时相的马尔科夫（MRF）和条件随机场（CRF）来建模上下文并在一定程度上提升了制图精度[15,16],但这些模型为了实现时空平滑效果,仅仅选择基于特征差异来计算上下文[17],而忽略了对地物的语义上下文特别是时间语义的变化关系建模。
对于耕地变化分析,主要是基于不同时相耕地提取结果,先统计耕地覆盖面积,后分析变化位置和因素。现有研究主要存在以下不足：① 以像素为分析单元,其结果易受像素随机噪声的影响[18]; ② 以市县或乡镇等行政区划为单元,往往尺度较大且无法分析行政区内的局部变化信息[19];③ 对于耕地变化的空间相关性,特别是局部空间自相关特点分析较少。因此,需要选取一个合适的时空一致性分析单元进行有效的耕地变化特点分析。
为解决上述问题,本文提出一种时空条件随机场模型来综合建模多时相影像的时空特征及语义上下文关系,提高耕地覆盖制图精度;同时基于耕地制图结果,构建规则格网作为变化分析单元,结合GIS空间分析对耕地面积变化规律和格局特点进行分析。最后采用2015—2019年北京市顺义区多时相Sentinel-2影像作为实验数据,对所提出的方法进行了验证分析。

2 研究方法

研究方法包括3个部分（图1）：① 多时相特征提取;② 耕地覆盖信息提取;③ 耕地面积变化分析。
图1 多时相时空上下文耕地提取和变化分析方法框架

Fig. 1 The framework of multi-temporal cultivated land extraction and change analysis based on spatiotemporal context

2.1 多时相特征提取

2.1.1 多时相光谱和纹理特征
多时相光谱和纹理特征记录了不同时间点地物的物理反射/辐射属性。本文所定义的多时相光谱特征包括：全波段的光谱特征（即红、绿、蓝、近红外和短波红外波段）;3个指数特征,包括归一化植被指数（NDVI）、增强植被指数（EVI）和地表水体指数（LSWI）。多时相纹理特征由灰度共生矩阵（GLCM）计算的5个特征值组成,窗口大小为 5×5（表1）。
表1 多时相分类所采用的特征介绍

Tab. 1 Introduction of features used in multi-temporal classification

类型 维度 描述
光谱特征 1—8 R, G, B, NIR, SWIR, NDVI, EVI, LSWI
纹理特征 9—12 GLCM-Homogeneity
13—16 GLCM-Correlation
17—20 GLCM-Energy
21—24 GLCM-Entropy
25—28 GLCM-Contrast
空间特征 29—36 面积、边界长度、长宽比、矩形度、椭圆度、紧凑度、光滑度和形状指数
2.1.2 多时相空间特征
多时相空间特征反映了地物在不同时间点的尺寸、形状和几何分布特点。但空间特征无法直接基于像素单元来表达,已有面向对象方法（MRS）[20]虽易于表达空间特征,但缺乏多时相影像统一分割方法,单独分割会在各时相生成不一致对象,从而无法建立对象的时间上下文关联,因此本文利用多时相像素的空间化描述来表达空间特征。首先对多时相影像的每景进行分别分割,对于时间t中的某个像素 p j t,将包含像素 p j t的对象 O i t的空间特征 F S ( O i t )赋予 p j t（式（1）,图2）。本文表达了多边形面积、边界长度、长宽比等8个主要空间特征（表1）。
F S p j t = F S O i t , O i t = ∑ j = 1 m p j t
图2 多时相像素空间化描述

Fig. 2 Spatialization of multi-temporal pixels

2.1.3 多时相时空上下文特征
多时相影像光谱、纹理和空间特征在不同时间的分布即可构成时间特征,可表达为地物在时间上下文中的特征依赖。本文将在时间和空间上分别表达地物的特征及语义上下文依赖,具体如2.2节所示。

2.2 耕地提取：时空上下文分类

时空上下文指空间和时间维度上相邻地物间的相关性。为了有效表达多时相影像分类中地物的时空上下文关系从而指导地物分类,本文将传统单时相CRF扩展至多时相域构建时空上下文CRF模型。与单时相CRF只考虑空间上下文不同,多时相影像分析必须同时考虑空间和时间维度的上下文关系。在具体定义中,对于包含N个影像的多时相数据,令所有像素位置集合为 P s,整个多时相特征集 F = f t | t = 1,2 , … , N, f t = { f i t ∈ R d | i ∈ P s }, f i t为时间点 t影像中位置 i像素的特征向量,维度为 d。所有像素的类别标记集合为 L = l i t | t ∈ T , i ∈ P s, l i t为时间 t影像中位置 i像素的类别, l i t ∈ C且 C = { 1,2 , … , M },共包含 M个类别。为建模地物在空间和时间上的上下文依赖,扩展的时空上下文CRF模型可由式（2）所示。
P L | F = 1 Z exp ∑ A l i t , f t + ∑ ∑ SC l i t , l j t , f t + ∑ ∑ TC l i t , l i u , f t , f u
式中： Z为归一化常量也称划分函数; A为特征关联,反映像素本身特征对类别的影响; SC（Spatial Context Correlation）为空间上下文关联; S i为当前像素 i的空间邻域像素集合; TC（Temporal Context Correlation）为时间上下文关联作用, T i为当前像素 i的时间邻域像素集合。对于当前时间 t的像素,空间邻域由时间 t内其前后左右4个最邻近像素组成;时间邻域由最邻近前后时间的2个像素组成（图3）。
图3 多时相时空邻域

Fig. 3 Multi-temporal spatiotemporal neighbourhoods

2.2.1 多时相特征关联
特征关联度量了语义标记和特征间的关系,因此 A ( l i t , f t )可表示为像素 i在特征 f t下取类别 l i t为 c的概率。通常而言,所有以概率输出为结果的判别式分类模型均可来建模 A ( l i t , f t )。本文选取概率输出的BP神经网络（BPNN） P BPNN来表达 A ( l i t , f t )[21],定义如式（3）所示。
A ( l i t , f t ) = - ln P l i t = c f t = - ln P BPNN l i t = c f t
2.2.2 多时相空间上下文关联
空间上下文指空间相邻像素在特征和语义上的相关性。目前常用方法是将空间上邻近地物分配相同类别,故也称空间平滑项。然而,空间相邻地物也存在类别变化,因此需要考虑语义关系。为综合考虑空间上特征和语义的关联,采用一种扩展的差异敏感性Potts模型[22]来建模空间上下文 （式（4））,在支持空间平滑效应同时,也兼顾了空间邻域的语义类别差异。
SC l i t , l j t , f t = δ s ∙ 1 - exp - ε s D ij t ( f t ) 2 R l i t = l j t δ s ∙ exp - ε s D ij t ( f t ) 2 R l i t ≠ l j t
式中： D ij t ( f t )为 D ij t ( f t )的欧式范式,表达空间相邻像素的特征差异,因此 D ij t ( f t ) = f i t - f j t;权重参数 δ s度量空间上下文在分类中的权重;交互系数 ε s调节特征依赖在 SC ( l i t , l j t , f t )中的影响; R表示 f t的维度数,用以使每个维度特征对于 SC ( l i t , l j t , f t )的贡献一致。
2.2.3 多时相时间上下文关联
时间上下文指时间相邻像素在特征和语义上的相关性。地物在不同时间的类别分布和变化往往呈现一定规律。例如耕地由于弃耕或者侵占更有可能转换为裸地或者建筑物而不是水体。这种时间语义关联能为耕地分类提供重要指导,因此需要在时间上下文中得以考虑。本文通过综合时间类别转换矩阵和特征差异构建了 TC ( l i t , l i u , f t , f u )（式（5））。
TC ( l i t , l i u , f t , f u ) = δ t ∙ ( 1 - TM ( l i t , l i u ) ) ∙ 1 - exp - ε t D l i t , l i u , f t , f u
式中： δ t和 ε t分别为权重和交互系数; TM ( l i t , l i u )为转换矩阵,每个元素都对应一个条件概率 P l i t = c i | l i u = c u,表示当前像素取类别 c i条件下,其时间邻域像素取类别 c u的概率; D度量了时间邻域间的特征差异,定义如式（6）所示。
D l i t , l i u , f t , f u = f i t - f i u - f t ¯ ( c i ) - f u ¯ ( c u ) R
式中： f t ¯ ( c i )和 f u ¯ ( c u )分别表示样本中在时间 t和 u上类别为 c i和 c u的所有像素的特征向量中心,因此 D在特征上惩罚了从类别 c i变换到 c u的误差,误差值越大则 D值越大。
2.2.4 模型参数的训练和推理
对于时空上下文CRF模型的训练,基于样本首先训练 A ( l i t , f t )中的BPNN,在权衡效率和精度下,设置了BPNN参数。包括：网络层数为3,隐层神经元个数为15,目标误差0.01。同时,基于交叉验证和迭代实验,设置 SC ( l i t , l j t , f t )中的 δ s = 0.75, ε s = 2.00, TC ( l i t , l i u , f t , f u )中 δ t = 1.50, ε t = 0.75,因为这些参数对应的迭代实验精度最高。考虑到 TM中参数的确定需大量有代表性的时间变化样本数据,本文基于顺义区地物变化的专家知识来设置 TM（表2）。 TM中较大值所对应的2类地物更可能发生相互变化,例如耕地更有可能变为裸土或建筑物而不是其他类别,整体而言大部分地物未发生类别变化,因此表2对角线的值最大。
表2 类别转换矩阵 TM中参数设置

Tab. 2 Parameters setting in class transition matrix TM

t t+1
耕地 林地 草地 裸土 建筑物 道路 水体
耕地 0.90 0.04 0.04 0.10 0.15 0.01 0.01
林地 0.02 0.90 0.02 0.01 0.01 0.01 0.01
草地 0.02 0.02 0.90 0.02 0.01 0.01 0.01
裸土 0.02 0.02 0.02 0.90 0.10 0.02 0.01
建筑物 0.01 0.01 0.01 0.08 0.90 0.05 0.05
道路 0.01 0.01 0.01 0.02 0.01 0.90 0.01
水体 0.01 0.01 0.01 0.01 0.01 0.01 0.90
推理的目的是在CRF关联函数的能量全局最小化下寻找最优的类别配置,本文采用基于α-膨胀算法的图割模型(Graph-cut model)进行推理[23]。

2.3 耕地时空变化分析

2.3.1 耕地变化分析单元和指标
由于已有行政区划和像素单元在尺度上存在着过大和过小的局限性,本文根据顺义区特点,构建了大小为1.55 km×1.55 km共468个规则格网,作为耕地变化分析的一致性单元。基于耕地覆盖结果,统计每年每个单元中的耕地面积,计算2015—2019年逐年各单元耕地面积变化率（式（7））,并将其作为变化分析指标。
R Crop T 1,2 = A T 2 - A T 1 A T 1
式中： R Crop T 1,2为相邻2年 T 1和 T 2的耕地面积变化率; A T 1和 A T 2为各年的单元耕地面积。
2.3.2 耕地变化分析方法
为挖掘耕地变化的空间分布,采用空间统计中的全局和局部空间自相关指标（Moran's I）来进行变化分析[24]。其中,全局自相关（ I Global,式（8））描述研究区整体耕地变化分布是否聚集,而局部自相关（ I Local i,式（9））能度量局部空间的聚集性。
I Global = N ∑ i ∑ j W ij ∙ ∑ i ∑ j W ij ( X i - X ̅ ) ( X j - X ̅ ) ∑ j ( X j - X ̅ ) 2
I Local i = X i - X ̅ ∑ j ≠ i N X j 2 ( N - 1 ) - X ̅ 2 ∑ j = 1 N W ij ( X j - X ̅ )
式中： N为格网总数; X i、 X j分别为格网 i和 j处的耕地面积变化率; X ̅为所有格网面积变化率的平均值; W ij为空间邻近度量指数。当格网 i和 j相邻时, W ij = 1;否则, W ij = 0。 I Global和 I Local i取值范围为(-1, 1),取正值且越大说明相邻格网正相关性越强,取负值且越小说明负相关越强,接近0时表明变化为随机分布特点。
此外,利用Z-Score模型进行Moran's I指数的标准化[25],以检验空间相关性在统计上的显著性。

3 实验区概况与数据来源

3.1 实验区概况

本文选取北京市东北部的顺义区为实验区 （图4）,范围为116°28′E—116°58′E,40°00′N—40°18′N,总面积1021.02 km2,涵盖耕地、建筑、森林、草地、水体等典型地物类型。研究区位于华北平原北端,农业生产条件优越,同时作为北京东北部发展带重要节点和发展新城之一,近年来城市建设和经济发展引起的耕地变化较为明显。
图4 北京市顺义区位置

注：子区域A、B将用于实验结果对比。

Fig. 4 Location of Shunyi District, Beijing

3.2 实验数据

所使用的多时相数据为北京市顺义区2015—2019年5景连续的Sentinel-2影像（数据来源于欧空局哥白尼数据中心 https://scihub.copernicus.eu/dhus/#/home[26]）,时间间隔1年,空间分辨率10 m（表3）。每景包含 4430 像素 × 3042像素,都保证云量尽可能少。影像的成像时间均选在8月,在这个时间段,研究区的主要农作物（玉米、高粱、小麦等）均处于生长最旺盛的抽穗和开花等时期,避免了作物收割造成的耕地误分影响,因此是识别耕地的最佳时间。所有影像都经过了几何校正和大气校正处理,有效消除了几何和辐射误差。
表3 2015—2019年多时相Sentinel-2数据信息

Tab. 3 Information of multi-temporal Sentinel-2 data from 2015 to 2019

参数 数值
空间分辨率 10 m
光谱波段 B2 “Blue” (490 nm)
B3 “Green” (560 nm)
B4 “Red” (665 nm)
B8 “Near-Infrared” (842 nm)
B6 “Short Wave Infrared 1” (1610 nm)
成像时间 2015-08-13, 2016-08-21, 2017-08-06, 2018-08-16, 2019-08-18

4 实验结果及分析

4.1 耕地覆盖提取结果

为验证方法的有效性,对比了本文提出的多时相时空上下文分类法(ST-PC)和已有的2种主流多时相制图方法：① 多时相面向对象分类法(NOC)[13],对多时相影像进行单独的对象分割和特征提取,并采用浅层全连接BP神经网络分类器进行分类; ② 未考虑空间特征的多时相像素上下文分类法(PCC)[15,16], 采用多时相CRF建模像素上下文。对于5个时期的每景影像,基于Google Earth、调查资料和专家知识的目视解译,选取7类监督样本。所选样本随机分布,并保持类别比例和实际相对一致（2017年的样本分布如图5所示）,总体划分的训练样本数量约为测试样本的2倍（表4）。最后,采用整体精度、每类的用户、制图精度和Kappa系数定量评价分类结果[25]。
图5 2017年实验区影像的样本空间分布

Fig. 5 Spatial distribution of samples of experimental area in 2017

表4 2015—2019年5个时相影像各类样本的总数

Tab. 4 The total number of samples of each class for 5 images from 2015 to 2019

样本/像素 类别
耕地 林地 草地 裸土 建成区 道路 水体
训练 1 382 063 478 291 165 725 115 176 1 359 133 448 467 171 082
测试 690 474 238 686 80 637 55 349 678 865 224 026 83 813
采用3种方法对5个时相影像做了分类。为便于分析,选取2个代表性子区域（图4中A和B）,对A区域2018年和B区域2019年的分类结果进行对比分析（图6）。对2个结果进行目视对比分析,面向对象的分类结果中存在普遍的植被类型误分现象,如耕地、林地和草地间的混淆。这是由于多时相面向对象方法单独分割每个时相,无法表达不同时相对象间的关联,因此难以借助地物的时间关联特征来有效分类光谱相似的植被类型。多时相像素上下文方法能有效表达地物的时间关联,因此能有效分类耕地、林地和草地,但由于未考虑地物空间分布,因此结果中存在着随机椒盐噪声,如耕地和建筑的椒盐破碎现象。整体而言,本文方法获得的多时相制图结果从目视上看效果最好。本文方法通过建模时空上下文信息消除了不同地物间的混淆分类现象,同时考虑空间特征有效地克服了椒盐噪声。
图6 3种方法在A区域2018年和B区域2019年的分类结果

Fig. 6 Classification results of three methods in 2018 of sub-region A and 2019 of sub-region B

表5为3种方法对5个时相影像的用户和制图精度,表6列出了3种方法对5个时相影像的整体分类精度和Kappa系数。面向对象方法对耕地、林地和草地上取得了最低的用户和制图精度,而在像素上下文方法中精度得以提升,证实了时间关联的重要性;但由于未考虑空间特征,像素上下文分类在空间异质性较强的建成区和道路上的精度相对最低。本文方法在耕地上取得了最高的用户（91.21%）和制图精度（90.53%）,同时除了水体的用户精度外,其他所有类型均取得了最高精度。此外,本文方法在每个时相上都取得了最高的整体精度和Kappa系数（表6）,5个时相平均整体精度为90.79%,比像素上下文方法高3.03%,比对象方法高6.13%。这表明本文方法比已有常用的2种方法精度更高,能更准确地获取多时相耕地覆盖结果。
表5 2015—2019年5个时相上各土地利用类别的平均用户精度和制图精度

Tab. 5 Average user's and producer's accuracies of each class in 5 times from 2015 to 2019 (%)

类别 分类方法
NOC PPC ST-PC
ST-CRF
用户
精度
耕地 81.27 88.14 91.21
林地 80.01 87.96 90.03
草地 80.36 87.32 88.92
裸土 87.19 89.25 90.85
建成区 88.82 84.16 91.84
道路 89.73 83.58 90.41
水体 96.79 95.06 96.62
制图
精度
耕地 80.83 87.75 90.53
林地 80.94 87.24 89.72
草地 81.59 88.38 90.98
裸土 86.93 88.62 91.06
建成区 88.95 84.76 91.17
道路 88.57 84.19 90.64
水体 96.02 95.86 96.93
表6 2015—2019年5个时相上3种分类的整体精度和Kappa系数

Tab. 6 The overall accuracies and Kappa coefficients of five times for three classifications from 2015 to 2019

年份 分类方法
NOC PPC ST-PC
OA K OA K OA K
2015 85.11 0.8462 87.37 0.8672 90.52 0.8993
2016 84.23 0.8359 88.28 0.8713 91.05 0.9113
2017 84.53 0.8427 87.54 0.8705 90.43 0.8947
2018 85.26 0.8471 88.32 0.8746 90.81 0.9016
2019 84.15 0.8303 87.27 0.8657 91.12 0.9092
平均值 84.66 0.8404 87.76 0.8699 90.79 0.9032
采用本文方法对5个时相影像进行耕地分类,得到最佳耕地覆盖结果。选取2018和2019年的结果进行展示分析（图7）。顺义区主要覆盖类型为耕地和建成区,耕地在全区呈现较为均匀的分布状态,林地主要分布在东北部。对耕地分析可得,2018—2019年间,顺义区耕地覆盖整体有所减少,减少较为明显的区域如图中红色框所示,减少的耕地在2019年主要转换为建成区,少部分转换为裸土。这些耕地覆盖结果将被用于后面的定量化耕地变化分析。
图7 2018年和2019年实验区耕地覆盖结果

Fig. 7 Results of cultivated land covers in 2018 and 2019

4.2 耕地变化空间格局分析

4.2.1 基于行政区划的变化分析
将2015—2019年顺义区耕地分类结果与顺义区的行政区划矢量图叠加,进行各乡镇和街道的耕地面积变化分析。表7和图8为顺义区各乡镇耕地面积统计结果。5年间各乡镇的耕地面积随时间都呈逐步下降趋势,其中杨镇地区、赵全营镇、李桥镇逐年下降幅度较大,有些地区耕地变化较小或者偶有面积增加的情况（如木林镇2016年到2017年面积有所增加）。按整个区域看（图9）,顺义区的整体耕地面积随时间不断下降,5年时间总面积减少了8270.20 hm2。结合分类结果分析可知,减少的耕地主要转换成建成区,少部分由于弃耕转成裸土,说明耕地下降主要由建设用地侵占引起。下降速率逐步放缓,说明顺义区的城市扩张减速,耕地侵占减弱。
表7 2015—2019年顺义区各乡镇耕地面积

Tab. 7 Cultivated land area of towns in Shunyi District from 2015 to 2019 （hm2）

乡镇 耕地面积
2015年 2016年 2017年 2018年 2019年
龙湾屯镇 955.91 875.59 675.13 695.92 642.11
北石槽镇 453.43 421.88 364.89 349.42 348.31
牛栏山地区 424.53 403.28 327.05 327.88 309.36
后沙峪地区 157.90 151.11 128.06 122.63 114.01
木林镇 2066.39 1624.53 1902.75 1778.64 1696.36
南法信地区 71.40 69.44 64.48 64.30 63.60
杨镇地区 5023.89 4761.19 4188.13 3803.21 3699.04
高丽营镇 2167.58 1740.40 1573.42 1554.92 1422.35
仁和地区 99.62 90.32 83.44 82.18 79.06
北小营镇 1007.41 922.73 636.51 635.73 639.36
南彩镇 629.17 628.17 518.16 517.08 518.78
李遂镇 2046.14 1975.79 1731.35 1722.20 1607.75
李桥镇 3170.52 2942.48 2465.86 1887.23 1878.01
北务镇 2409.29 2239.89 2136.82 2087.35 2035.42
天竺地区 53.52 47.11 42.96 42.93 42.74
大孙各庄镇 2747.66 2796.68 2421.17 2318.53 2103.64
赵全营镇 3171.98 2532.47 2088.10 1672.52 1501.54
张镇 1287.64 1272.60 1150.34 1061.18 987.65
马坡地区 124.73 124.25 121.81 121.35 110.45
街道区域 10.06 9.71 9.62 9.60 9.07
总面积 28 078.79 25 629.62 22 630.05 208 54.81 19 808.59
图8 2015—2019年顺义区各乡镇耕地面积

Fig. 8 Cultivated land area of towns in Shunyi District from 2015 to 2019

图9 2015—2019年顺义区耕地面积变化

Fig. 9 Change of cultivated land area in Shunyi District from 2015 to 2019

4.2.2 基于规则格网的变化分析
依据3.3.1节与3.3.2节中的方法,利用ArcGIS 10.3中的空间统计工具计算全局和局部的Moran's I值。表8为顺义区各格网耕地面积变化率的全局Moran's I结果。2015—2019年逐年耕地面积变化率全局Moran's I值均大于0,且Z-Score值均大于1.96,说明各格网间耕地变化存在显著的空间正相关性,面积变化整体具有聚集现象。2015—2016年的Moran's I值最高,说明这一年变化的聚集程度最大;随后2016—2017年和2017—2018年变化的聚集度减小,2018—2019年聚集度又呈现增加趋势。
图10为顺义区耕地面积变化率局部空间自相关的Z-score标准化结果。其中,Z-score值为 ( - 1.65,1.65 )时,表明格网单元的耕地面积变化没有相关性;值为 ( - 1.96 , - 1.65 )或 ( 1.65,1.96 )时,在0.10水平上呈现显著空间相关性;值为 ( - 2.58 , - 1.96 )或 ( 1.96,2.58 )时,在0.05水平显著;值小于 - 2.58或大于 2.58时,在0.01水平显著。同时,Z-score为正且越大时正相关性越强,说明单元间的耕地变化情况一致,如红色区域所示;Z-score为负且越小时负相关性越强,如蓝色区域所示。
图10 2015—2019年顺义区耕地面积变化局部Z-score分布

Fig. 10 Local Z-score distribution of cultivated land area change in Shunyi District from 2015 to 2019

2015—2016年顺义区的耕地面积变化存在 4个正相关聚集区,分布在赵全营镇、高丽营镇、木林镇和杨镇地区;2016—2017年耕地变化存在6个正相关聚集区,分别为赵全营镇、后沙峪地区、北小营镇、李桥镇、木林镇以及杨镇地区;2017—2018年的2个主要正相关聚集区分布在赵全营镇和李桥镇,同时在杨镇地区存在一个负相关聚集区;2018—2019年耕地变化的6个正相关聚集区主要分布在赵全营镇、高丽营镇、龙湾屯镇、张镇、大孙各庄镇和李遂镇。
同时获取的还有顺义区耕地面积变化率的LISA聚集结果（图11）,其中HH表示格网耕地面积增加的位置其邻域也增加,LL表示面积减少的位置其邻域也减少,HL表示面积增加的位置其邻域减少,LH表示面积减少的位置其邻域增加。由结果可知,5年间,顺义区耕地面积整体呈现集中减少的趋势,各年份的LISA图中大部分呈现为LL聚集,如蓝色区域所示。集中减少较为明显的区域为赵全营镇、高丽营镇、木林镇和杨镇。以赵全营镇为例,其2015—2016年耕地覆盖及变化结果如图12所示,分析可知此期间赵全营镇减少的耕地主要转换为建成区（如图12中红色区域所示）。由于顺义区工业化的发展,大量园区外的村镇企业不断涌现侵占农田,以赵全营镇、木林镇和杨桥镇工业发展较为迅速,因此耕地面积集中减少。同时部分地区在某些年份也会出现HH聚集现象,如2015—2016年的大孙各庄镇、2016—2017年和2017—2018年的木林镇等。说明这些地区对于耕地保护管理工作做的较好,使耕地的面积有小幅增大,但整体远小于耕地减少的速度。
图11 2015—2019年顺义区耕地面积变化LISA聚集分布

Fig. 11 LISA distribution of cultivated land area change in Shunyi District from 2015 to 2019

图12 2015—2016年赵全营镇耕地覆盖及变化结果

Fig. 12 Results of cultivated land cover and change in Zhaoquanying town from 2015 to 2016

图13为顺义区各乡镇逐年耕地面积变化率结果。从变化率柱状图的起伏程度分析,5年间耕地面积变化幅度较大的是木林镇、高丽营镇、北小营镇、李桥镇和赵全营镇,与局部空间自相关结果相符。在大孙各庄镇（2015—2016年）和木林镇（2016—2017年）中存在变化率大于零的情况,说明这两地的耕地面积偶有增加,耕地的保护取得一定成效,这也与局部空间自相关结果相符。
图13 2015—2019年顺义区各乡镇逐年耕地面积变化率

Fig. 13 Annual change rate of cultivated land area of Shunyi District from 2015 to 2019

5 结论与讨论

基于多时相遥感影像的耕地分类和变化研究可准确、有效地获取不同时期耕地的覆盖和变化信息。现有方法在多时相耕地的分类和变化分析上都存在一定的局限性,导致分析精度不佳。因此,本文提出了一种多时相时空上下文的耕地覆盖分类方法,综合表达和利用了多时相遥感影像中地物的光谱、纹理、空间特征以及时空上下文关联信息,获得了精确的多时相耕地提取结果。同时基于耕地覆盖结果,在规则格网和行政区划单元上,结合GIS空间统计对耕地变化的空间分布和格局特点进行了分析。研究以北京市顺义区2015—2019年的耕地覆盖和变化为切入点,以多时相的Sentinel-2影像为数据源进行了实验验证,得到如下结论：
（1）在多时相影像分类中,利用空间特征能有效地消除空间异质性引起的椒盐随机噪声。如本文方法与已有未考虑空间特征的多时相像素上下文分类方法相比,有效地消除了椒盐噪声,耕地的平均制图精度也从87.75%提升至90.53%。
（2）地物的时空语义和特征上下文信息能有效地提升多时相制图精度,特别是对于光谱相似的植被覆盖类型。与未考虑上下文的面向对象方法相比,像素上下文和本文方法显著提升了耕地、林地和草地的分类精度,幅度为6.30%~9.70%。
（3）本文提出的方法优于已有的2种常见多时相影像分类方法（多时相面向对象分类法和未考虑空间特征的像素上下文分类法）。与已有2种方法对比,本文方法在耕地上获得了最高的用户精度（91.21%）和制图精度（90.53%）,整体精度也最高（90.79%）。
（4）基于规则格网的空间统计方法能有效地挖掘出顺义区耕地变化的空间规律和格局特点。结果表明,2015—2019年顺义区耕地面积整体呈现逐年下降的趋势,变化存在区域聚集现象,且主要为集中减少特点。其中集中减少较为明显的区域为赵全营镇、高丽营镇、木林镇和杨镇地区。虽然部分地区的耕地面积偶有小幅增加现象,但整体远小于耕地减少的速度。
因此,建议顺义区政府在城市化发展的同时,应实施更为严格的耕地保护政策,特别是对于赵全营镇、高丽营镇、木林镇和杨镇等几个较为明显的耕地减少地区,要加大国家耕地占补平衡和土地开发整理的力度,实现耕地资源的科学保护和管理。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
Smil V. China's energy and resource uses: continuity and change[J]. The China Quarterly, 1998, 156:935-951.

DOI

[2]
毕于运. 中国耕地[M]. 北京: 中国农业科技出版社, 1995:1-9.

[ Bi Y Y. China cultivated land[M]. Beijing: China Agricultural Science and Technology Publishing, 1995:1-9. ]

[3]
赵其国, 周生路, 吴绍华, 等. 中国耕地资源变化及其可持续利用与保护对策[J]. 土壤学报, 2006, 43(4):662-672.

[ Zhao Q G, Zhou S L, Wu S H, et al. Cultivated land resources and strategies for its sustainable utilization and protection in China[J]. Acta Pedologica Sinica, 2006, 43(4):662-672. ]

[4]
张国平, 刘纪远, 张增祥. 近10年来中国耕地资源的时空变化分析[J]. 地理学报, 2003, 58(3):323-332.

[ Zhang G P, Liu J Y, Zhang Z X. Spatial-temporal changes of cropland in China for the past 10 years based on remote sensing[J]. Acta Geographica Sinica, 2003, 58(3):323-332. ]

[5]
Lillesand T, Kiefer R W, Chipman J. Remote sensing and image interpretation[M]. John Wiley & Sons, 2015.

[6]
Giuliani G, Dao H, De Bono A, et al. Live Monitoring of Earth Surface (LiMES): A framework for monitoring environmental changes from Earth Observations[J]. Remote Sensing of Environment, 2017, 202:222-233.

DOI

[7]
Xi W, Du S. What information is important? A spatiotemporal inference for classification of satellite image time series[C]. IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium. IEEE, 2019:2415-2418.

[8]
魏鹏飞, 徐新刚, 杨贵军, 等. 基于多时相影像植被指数变化特征的作物遥感分类[J]. 中国农业科技导报, 2019, 21(2):54-61.

[ Wei P F, Xu X G, Yang G J, et al. Crop remote sensing classification based on vegetation index change characteristics of multi temporal images[J]. Chinese Journal of Agricultural Science and Technology, 2019, 21(2):54-61. ]

[9]
Slagter B, Tsendbazar N E, Vollrath A, et al. Mapping wetland characteristics using temporally dense Sentinel-1 and Sentinel-2 data: A case study in the St. Lucia wetlands, South Africa[J]. International Journal of Applied Earth Observation and Geoinformation, 2020, 86:102009.

DOI

[10]
阴海明, 王立辉, 董明霞, 等. 基于多时相Sentinel-2遥感影像的江汉平原夏收作物提取方法[J]. 福建农林大学学报(自然科学版), 2021, 50(1):16-22.

[ Ying H M, Wang L H, Dong M X, et al. Extraction method of summer harvest crops in Jianghan plain based on multi-temporal Sentinel-2 remote sensing images[J]. Journal of Fujian Agriculture and Forestry University (Natural Science Edition), 2021, 50(01):16-22. ]

[11]
Peña M A, Brenning A. Assessing fruit-tree crop classification from Landsat-8 time series for the Maipo Valley, Chile[J]. Remote Sensing of Environment, 2015, 171:234-244.

DOI

[12]
Xi W, Du S, Wang Y C, et al. A spatiotemporal cube model for analyzing satellite image time series: Application to land-cover mapping and change detection[J]. Remote Sensing of Environment, 2019, 231:111212.

DOI

[13]
Blaschke T. Object based image analysis for remote sensing[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2010, 65(1):2-16.

DOI

[14]
Petitjean F, Kurtz C, Passat N, et al. Spatio-temporal reasoning for the classification of satellite image time series[J]. Pattern Recognition Letters, 2012, 33(13):1805-1815.

DOI

[15]
Hoberg T, Rottensteiner F, Feitosa R Q, et al. Conditional random fields for multitemporal and multiscale classification of optical satellite imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 53(2):659-673.

DOI

[16]
Wehmann A, Liu D. A spatial-temporal contextual Markovian kernel method for multi-temporal land cover mapping[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 107:77-89.

DOI

[17]
Liu D S, Cai S S. A spatial-temporal modeling approach to reconstructing land-cover change trajectories from multi-temporal satellite imagery[J]. Annals of the Association of American Geographers, 2012, 102(6):1329-1347.

DOI

[18]
Giri C, Muhlhausen J. Mangrove forest distributions and dynamics in Madagascar (1975—2005)[J]. Sensors (Basel, Switzerland), 2008, 8(4):2104-2117.

DOI

[19]
柏林川, 武兰芳, 宋小青. 1995—2010年山东省粮食单产变化空间分异及均衡增产潜力[J]. 地理科学进展, 2013, 32(8):1257-1265.

[ Bo L C, Wu L F, Song X Q. Spatial difference of grain yield changes during 1995-2010 and balanced potential output to increase in Shandong Province[J]. Progress in Geography, 2013, 32(8):1257-1265. ]

[20]
Baatz M. Multi resolution segmentation: An optimum approach for high quality multi scale image segmentation[C]. Beutrage zum AGIT-Symposium. Salzburg, Heidelberg, 2000:12-23.

[21]
Ding S F, Su C Y, Yu J Z. An optimizing BP neural network algorithm based on genetic algorithm[J]. Artificial Intelligence Review, 2011, 36(2):153-162.

DOI

[22]
Shotton J, Winn J, Rother C, et al. Texton Boost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[J]. International Journal of Computer Vision, 2009, 81(1):2-23.

DOI

[23]
Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graph cuts[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2001, 23(11):1222-1239.

DOI

[24]
Fotheringham A S, Brunsdon C, Charlton M. Quantitative geography: Perspectives on spatial data analysis[M]. London: SAGE Publication, 2000.

[25]
Foody G M. Classification accuracy comparison: Hypojournal tests and the use of confidence intervals in evaluations of difference, equivalence and non-inferiority[J]. Remote Sensing of Environment, 2009, 113(8):1658-1663.

DOI

[26]
Suhet. Sentinel-2 User Handbook[EB/OL]. https://sentinel.esa.int/documents/247904/685211/Sentinel-2_User_H and book, 2013-9-1.

Options
Outlines

/

〈 〉