基于大数据的长时序国土空间地域功能结构变化检测方法转让专利
申请号 : CN202110620580.4
文献号 : CN113255808B
文献日 : 2021-12-21
发明人 : 王亚飞 , 申振宇
申请人 : 中国科学院地理科学与资源研究所
摘要 :
本发明涉及一种基于大数据的长时序国土空间地域功能结构变化检测方法，方法如下：首先，构建以开发和保护功能两个维度为指向的功能用地分类体系，基于Google Earth Engine云平台采用随机森林算法对时间序列影像进行初始分类。其次，利用时空一致性双向检测方法，对初始分类结果进行一致性逻辑纠错处理，并在此基础上识别时间序列变化的突变点。最后，以建设用地、农业用地、生态用地比例关系为基本参数，测算区域不同时期点‑轴结构和核心‑边缘等地域功能结构。这种方法提高了长时序功能用地分类与变化检测精度，尤其适用于多种区域政策作用下导致的各类地物间的多次转变状况，能够直接应用于国土空间规划战略格局刻画以及区域政策的辅助决策。
权利要求 :
1.一种基于大数据的长时序国土空间地域功能结构变化检测方法，包括以下步骤：步骤1、数据准备——借助遥感大数据平台Google Earth Engine收集指定时间范围和研究范围的Landsat影像数据，并进行影像拼接和裁剪，以获得时间序列影像，其中每一年对应一幅影像；
步骤2、样本点选取——通过实地采样或借助高分辨率影像选取指定时间范围内最后一年各类地物的样本点，并以此年为基准年，切换上一年参考影像，逐点比对，对地物类别发生改变的样本点进行修改，从而得到上一年的样本点，以此类推，直到完成全部年份的样本点选取；
步骤3、影像初始分类——将每幅影像的样本点按比例随机划分为训练集和验证集，用历年训练集对分类器进行训练，训练完成后，使用分类器对对应年份的影像进行分类，获得初始分类结果，分类结果具有相同的分辨率;步骤4、分类预处理——使用众数滤波工具对每幅影像的分类结果进行迭代处理，去除碎小图斑；
步骤5、分类结果合成——将所有年份的影像分类结果叠加，构建分类结果的时间序列栅格数据集，时间序列即为取该栅格数据集中某一行某一列像素的多年集合；
步骤6、双向时空一致性检测——在每个时间序列中分别构造两个种子窗口Ws和两个滑动检测窗口Wd，两个种子窗口Ws的初始位置位于时间序列的首末端，两个滑动检测窗口Wd的初始位置位于紧邻于对应的种子窗口Ws，按照以下步骤进行处理：a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；
b）、将滑动窗口置于种子窗口紧邻的内侧；
c）、统计滑动窗口内的主导地类Kd，窗口内某地类出现频率高于60%则称其为该窗口的主导地类Kd；若窗口内所有地类出现频率均低于60%，则该窗口没有主导地类；
d）、判断种子窗口Ws内地物类型Ks与紧邻的滑动检测窗口Wd内的主导地类Kd是否一致，若一致，则将滑动检测窗口Wd内首尾主导地类栅格之间所有地类均置为Kd，将种子窗口位置移动到滑动窗口内最后出现主导地类Kd的栅格位置，转至步骤b)直到两个滑动检测窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口相遇；
e）、当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口的像素地类一致的像素向左侧种子窗口方向移动，大窗口内与右侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年；
步骤7、获得修正后的分类结果——将时间序列栅格数据集按顺序拆分，同一年份的栅格单元构成该年份的分类结果，最终获得所有年份修正后的分类结果。
2.根据权利要求1所述的基于大数据的长时序国土空间地域功能结构变化检测方法，其特征在于：对分类器训练包含以下几个步骤：
1）、训练样本点分布优化，使用随机函数对训练样本点进行若干次随机位置分布，取精度最高的结果，作为最终训练样本；
2）、引入数种纹理特征，试验数种尺寸窗口，取精度最高的结果，作为最终纹理特征；
3）、引入数种气候要素特征，取精度最高的结果，作为最终气候特征；
4）、进行使用步骤1）的最终训练样本、步骤2）的最终纹理特征和步骤3）的最终气候特征，对分类器进行训练。
3.根据权利要求1所述的基于大数据的长时序国土空间地域功能结构变化检测方法，其特征在于：地物分类如下：以开发和保护功能两个维度为指向，将建设用地、农业用地作为开发指向的功能类型，将生态用地作为保护功能指向的功能类型；其中，农业用地是指耕地，生态用地包括林地、草地、灌木、水体和未利用地。
4.根据权利要求1所述的基于大数据的长时序国土空间地域功能结构变化检测方法，其特征在于：分类器选用随机森林分类器。
5.根据权利要求1所述的基于大数据的长时序国土空间地域功能结构变化检测方法，其特征在于：步骤7完成后对分类结果进行精度评价，具体的，使用验证集对分类结果进行精度评价，利用混淆矩阵计算生产者精度和用户精度对各类地物进行评价，并算出总体精度和kappa系数。
6.根据权利要求1所述的基于大数据的长时序国土空间地域功能结构变化检测方法，其特征在于：步骤7完成后，进行功能用地分类与结构变化检测，具体如下：
1）、长时序功能用地分类——以建设用地、农业用地、生态用地为一级分类进行功能类型合并，计算不同时期功能用地的面积及比重，进行长时序功能用地制图；
2）、长时序地域功能结构测算——以建设用地、农业用地、生态用地比例关系为基本参数，选取交通干线为轴线，以轴线和核心城市为中心设定不同缓冲区，统计不同缓冲区内建设用地、农业用地、生态用地比例关系，测算区域不同时期点‑轴结构和核心‑边缘地域功能结构。
说明书 :
基于大数据的长时序国土空间地域功能结构变化检测方法
技术领域
[0001] 本发明涉及一种遥感大数据平台Google Earth Engine（GEE）支持下的长时序功能用地检测方法，特别是涉及一种基于像素级时空一致性双向检测的长时序城镇建设、农
业、生态用地功能结构变化检测方法。
业、生态用地功能结构变化检测方法。
背景技术
[0002] 土地覆盖/利用功能是人类活动与地表自然生态系统长期相互作用的结果，其动态变化是社会‑生态复合系统可持续与否的重要表征，常作为国家和地方政府城镇化发展、
生态环境保护与监测、国土空间规划、可持续发展目标测度等基础数据和重要指标。利用多
时相遥感影像数据进行地表土地覆盖/利用功能的分类制图及变化检测一直是长期应用的
热点问题。随着遥感技术的发展以及各种类型遥感数据的历史积累，时间序列影像越来越
丰富，使得借助于地物季节性、长期稳定性等遥感影像时间序列特征的长时序功能用地变
化检测成为可能。
生态环境保护与监测、国土空间规划、可持续发展目标测度等基础数据和重要指标。利用多
时相遥感影像数据进行地表土地覆盖/利用功能的分类制图及变化检测一直是长期应用的
热点问题。随着遥感技术的发展以及各种类型遥感数据的历史积累，时间序列影像越来越
丰富，使得借助于地物季节性、长期稳定性等遥感影像时间序列特征的长时序功能用地变
化检测成为可能。
[0003] 近年来，新开发的地理空间数据分析云平台Google Earth Engine（GEE），改变了传统的遥感处理方法，其庞大的遥感影像数据集和高性能的计算能力，为长时期、大规模的
遥感变化监测分析提供了一种新途径。借助于GEE平台进行长时序地表覆盖与土地利用功
能分析已取得了一些突破性的进展，主要包括两类：一类是基于时间序列特征指数的单一
地物或同类地物研究，通过重构地物生产过程或利用季节性、物候等时间序列特征进行识
别与变化检测。如2019年Jin等在《International Journal of Remote Sensing》撰文
“Mapping the annual dynamics of cultivated land in typicalarea of the Middle‑
lower Yangtze plain using long time‑series of Landsat images based on Google
Earth Engine”，通过增强时空自适应反射率融合模型(ESTRAFM)分别构建年度植被差异指
数(NDVI)轨迹和合成NDVI，利用时序断点进行植被和耕地的制图与变化检测；2020年
Francisco等在《IEEE Transactions on Geoscience and Remote Sensing》撰文“EVI
Time‑Series Breakpoint Detection Using Convolutional Networks for Online
Deforestation Monitoring in Chaco Forest”，利用时序增强植被指数（EVI）和卷积神经
网络模型评估森林砍伐时间，进行森林制图与砍伐的变化检测。另一类是基于直接分类或
数据平滑的全覆盖的功能谱系分类与变化检测。如2020年Nishanta在《Remote Sensing》撰
文“A Comparison of Three Temporal Smoothing Algorithms to Improve Land Cover
Classification:A Case Study from NEPAL”，对比了三种平滑方案傅里叶变换平滑、惠特
克平滑和线性拟合，对尼布尔Nepal第一省2000‑2018的土地覆被数据进行了平滑，将整体
土地覆盖地图精度从79.18%提高到83.44%。在2020年Tzu‑Hsin Karen Chen在《Remote
Sensing》撰文“Mapping horizontal and vertical urban densification in Denmark
with Landsat time‑series from 1985 to 2018: A semantic segmentation solution”
提出了基于像素级语义分割框架的卷积神经网络（CNNs），对1985‑2018年间丹麦城市不透
水面进行水平和垂直方向两个维度进行遥感识别，与完全卷积网络(FCN)和随机森林(RF)
模型对比，在水平和垂直向，F1分数可以增加4个百分点和10个百分点。
遥感变化监测分析提供了一种新途径。借助于GEE平台进行长时序地表覆盖与土地利用功
能分析已取得了一些突破性的进展，主要包括两类：一类是基于时间序列特征指数的单一
地物或同类地物研究，通过重构地物生产过程或利用季节性、物候等时间序列特征进行识
别与变化检测。如2019年Jin等在《International Journal of Remote Sensing》撰文
“Mapping the annual dynamics of cultivated land in typicalarea of the Middle‑
lower Yangtze plain using long time‑series of Landsat images based on Google
Earth Engine”，通过增强时空自适应反射率融合模型(ESTRAFM)分别构建年度植被差异指
数(NDVI)轨迹和合成NDVI，利用时序断点进行植被和耕地的制图与变化检测；2020年
Francisco等在《IEEE Transactions on Geoscience and Remote Sensing》撰文“EVI
Time‑Series Breakpoint Detection Using Convolutional Networks for Online
Deforestation Monitoring in Chaco Forest”，利用时序增强植被指数（EVI）和卷积神经
网络模型评估森林砍伐时间，进行森林制图与砍伐的变化检测。另一类是基于直接分类或
数据平滑的全覆盖的功能谱系分类与变化检测。如2020年Nishanta在《Remote Sensing》撰
文“A Comparison of Three Temporal Smoothing Algorithms to Improve Land Cover
Classification:A Case Study from NEPAL”，对比了三种平滑方案傅里叶变换平滑、惠特
克平滑和线性拟合，对尼布尔Nepal第一省2000‑2018的土地覆被数据进行了平滑，将整体
土地覆盖地图精度从79.18%提高到83.44%。在2020年Tzu‑Hsin Karen Chen在《Remote
Sensing》撰文“Mapping horizontal and vertical urban densification in Denmark
with Landsat time‑series from 1985 to 2018: A semantic segmentation solution”
提出了基于像素级语义分割框架的卷积神经网络（CNNs），对1985‑2018年间丹麦城市不透
水面进行水平和垂直方向两个维度进行遥感识别，与完全卷积网络(FCN)和随机森林(RF)
模型对比，在水平和垂直向，F1分数可以增加4个百分点和10个百分点。
[0004] 相对于单一地物或同类地物在利用季节性、物候等时间序列特征进行识别与变化检测取得突出进展的同时，全覆盖的功能谱系分类与变化检测在处理不同地物变化与冲突
方面，仍有较大的探索空间。突出表现为：对于单向稳定变化的地物类型，比如城镇扩张、森
林砍伐等，通常被认为是侵占耕地或森林生态系统的单向过程。这些研究采用的主要处理
方式是，构建指数时间序列重构地物生长过程，通过检测时间突变点汇总地类转化情况。而
现有的全覆盖分类与变化检测分析，通常是利用人工智能的相关算法进行直接分类，一部
分研究直接忽略了不同地物间符合现实逻辑的转变方式，因此，分类的结果往往导致同一
地物在时间序列上产生较大的逻辑错误，比如不同年度间建设用地规模存在较大的波动，
且出现大面积的不合理的逆城市化过程，即最新年份的城镇建设用地面积反而小于较远年
份的情况。另一部分研究简化了这种地物间转变过程，通过全地类要素的平滑过程，使其逻
辑在某种程度上符合现实的情况，采用一种类似于单一地物类型变化检测的方法）。而更多
的现实是，有些区域由于受各种政策的影响，存在耕地、林草等不同地物双向、多次转变的
现象，甚至建设用地被恢复为生态用地的情况。对于这一类区域，以上这些方法就严重制约
了其分类与变化检测的精度。
方面，仍有较大的探索空间。突出表现为：对于单向稳定变化的地物类型，比如城镇扩张、森
林砍伐等，通常被认为是侵占耕地或森林生态系统的单向过程。这些研究采用的主要处理
方式是，构建指数时间序列重构地物生长过程，通过检测时间突变点汇总地类转化情况。而
现有的全覆盖分类与变化检测分析，通常是利用人工智能的相关算法进行直接分类，一部
分研究直接忽略了不同地物间符合现实逻辑的转变方式，因此，分类的结果往往导致同一
地物在时间序列上产生较大的逻辑错误，比如不同年度间建设用地规模存在较大的波动，
且出现大面积的不合理的逆城市化过程，即最新年份的城镇建设用地面积反而小于较远年
份的情况。另一部分研究简化了这种地物间转变过程，通过全地类要素的平滑过程，使其逻
辑在某种程度上符合现实的情况，采用一种类似于单一地物类型变化检测的方法）。而更多
的现实是，有些区域由于受各种政策的影响，存在耕地、林草等不同地物双向、多次转变的
现象，甚至建设用地被恢复为生态用地的情况。对于这一类区域，以上这些方法就严重制约
了其分类与变化检测的精度。
[0005] 对于监测时间跨度长、监测年份频率高的遥感影像时间序列来说，如何在消除分类误差对地物间转化过程干扰的同时，又合理保留地物间转化的趋势方向，往往成为决定
功能用地分类制图与变化检测精度的关键。为了能够较为准确的检测不同土地覆盖类型存
在的多次转变，本文提出了一种时空一致性双向检测支持下的长时序国土空间地域功能用
地变化检测方法，对遥感影像时间序列初始分类结果进行一致性逻辑纠错处理，使时间序
列的分类结果符合土地利用变化的逻辑，并在此基础上识别国土空间功能变化的突变点。
该方法适用于长时序功能用地的分类制图与变化检测，尤其是多种区域政策作用下导致的
各类地物的多次转变状况。
功能用地分类制图与变化检测精度的关键。为了能够较为准确的检测不同土地覆盖类型存
在的多次转变，本文提出了一种时空一致性双向检测支持下的长时序国土空间地域功能用
地变化检测方法，对遥感影像时间序列初始分类结果进行一致性逻辑纠错处理，使时间序
列的分类结果符合土地利用变化的逻辑，并在此基础上识别国土空间功能变化的突变点。
该方法适用于长时序功能用地的分类制图与变化检测，尤其是多种区域政策作用下导致的
各类地物的多次转变状况。
发明内容
[0006] 本发明要解决技术问题是：针对长时间序列、高检测频率土地覆盖/利用功能数据出现的时间和空间维度的逻辑冲突问题，对全要素地物分类结果进行修正，以保证长时间
序列功能用地的时空逻辑一致性，提高分类结果精度。
序列功能用地的时空逻辑一致性，提高分类结果精度。
[0007] 为了解决上述技术问题，本发明提出基于大数据的长时序国土空间地域功能结构变化检测方法，包括以下步骤：
[0008] 步骤1、数据准备——借助大数据平台Google Earth Engine收集指定时间范围和研究范围的Landsat影像数据，并进行影像拼接和裁剪，以获得时间序列影像，其中每一年
对应一幅影像；
对应一幅影像；
[0009] 步骤2、样本点选取——通过实地采样或借助高分辨率影像选取指定时间范围内最后一年各类地物的样本点，并以此年为基准年，切换上一年参考影像，逐点比对，对地物
类别发生改变的样本点进行修改，从而得到上一年的样本点，以此类推，直到完成全部年份
的样本点选取；
类别发生改变的样本点进行修改，从而得到上一年的样本点，以此类推，直到完成全部年份
的样本点选取；
[0010] 步骤3、影像初始分类——将每幅影像的样本点按比例随机划分为训练集和验证集，用历年训练集对大数据分类器进行训练，训练完成后，使用大数据分类器对对应年份的
影像进行分类，获得初始分类结果，分类结果具有相同的分辨率;
影像进行分类，获得初始分类结果，分类结果具有相同的分辨率;
[0011] 步骤4、分类预处理——使用众数滤波工具对每幅影像的分类结果进行迭代处理，去除碎小图斑；
[0012] 步骤5、分类结果合成——将所有年份的影像分类结果叠加，构建分类结果的时间序列栅格数据集，时间序列即为取该栅格数据集中某一行某一列像素的多年集合；
[0013] 步骤6、双向时空一致性检测——在每个时间序列中分别构造两个种子窗口Ws和两个滑动检测窗口Wd，两个种子窗口Ws的初始位置位于时间序列的首末端，两个滑动检测窗
口Wd的初始位置位于紧邻于对应的种子窗口Ws，按照以下步骤进行处理：
口Wd的初始位置位于紧邻于对应的种子窗口Ws，按照以下步骤进行处理：
[0014] a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；
[0015] b）、将滑动窗口置于种子窗口紧邻的内侧；
[0016] c）、统计滑动窗口内的主导地类Kd，窗口内某地类出现频率高于60%则称其为该窗口的主导地类Kd；若窗口内所有地类出现频率均低于60%，则该窗口没有主导地类；
[0017] d）、判断种子窗口Ws内地物类型Ks与紧邻的滑动检测窗口Wd内的主导地类Kd是否一致，若一致，则将滑动检测窗口Wd内首尾主导地类栅格之间所有地类均置为Kd，将种子窗
口位置移动到滑动窗口内最后出现主导地类Kd的栅格位置，转至步骤b)直到两个滑动检测
窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口
相遇；
口位置移动到滑动窗口内最后出现主导地类Kd的栅格位置，转至步骤b)直到两个滑动检测
窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口
相遇；
[0018] e）、当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口的像素地类一致的像素向左侧种子窗口方向移动，大窗口内与右
侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将
其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年；
侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将
其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年；
[0019] 步骤7、获得修正后的分类结果——将时间序列栅格数据集按顺序拆分，同一年份的栅格单元构成该年份的分类结果，最终获得所有年份修正后的分类结果。
[0020] 本发明方法在消除分类误差对地物间转化过程干扰的同时，又合理保留地物间转化的趋势方向，提高了长时序功能用地分类与变化检测精度。该方法适用于长时序功能用
地的分类制图与变化检测，尤其是多种区域政策作用下导致的各类地物的多次转变状况，
能够直接应用于国土空间规划战略格局刻画以及区域政策的辅助决策。
地的分类制图与变化检测，尤其是多种区域政策作用下导致的各类地物的多次转变状况，
能够直接应用于国土空间规划战略格局刻画以及区域政策的辅助决策。
附图说明
[0021] 下面结合附图对本发明作进一步的说明。
[0022] 图1是本发明实施例研究区概况图。
[0023] 图2 本发明方法的流程图。
[0024] 图3 是双向检测算法流程图。
[0025] 图4（a）是处理后分类结果图（1991年）。
[0026] 图4（b）是处理后分类结果图（2000年）。
[0027] 图5是长时序功能分类制图。
[0028] 图6是双向检测精度与初始分类精度对比。
[0029] 图7是1987‑2020不同时期功能变化比重图。
具体实施方式
[0030] 本实验采用的实验数据为1987—2020年间27年的Landsat系列影像，地表反射率产品，30m分辨率。如图1所示，研究区选取青海省湟水流域，湟水流域地处青藏高原东北部，
地理位置36°02′37°28′ N，100°42′103°04′ E，流域面积约16200km²。流域呈羽状水系，
~ ~
地势西北高、东南低，干流自西北流向东南。流域内地貌类型涉及河谷平原、黄土丘陵、高山
牧场以及中、高山地，最低点海拔1612m，最高点海拔4860m，具有较为明显的垂直差异性。湟
水流域总面积仅占青海省的2.24%，但人口、工农业生产总值等均占整个青海省的50%以上，
是青海省人口密度最大、经济最为发达地区。与此同时，湟水流域年降水量300 500mm，位于
~
干旱半干旱气候地区，是中国400mm等降雨量生态过渡地带与典型的农牧交错带。一方面，
这个区域处在不断集聚人口经济的过程中，由此产生的城镇扩张与耕地保护之间的矛盾重
重。另一方面，农牧交错带本身生态极为脆弱，是中国政府主导的退耕还林、退耕还草的政
策区域，耕地生产与生态脆弱区的保护之间的矛盾也极为尖锐。
地理位置36°02′37°28′ N，100°42′103°04′ E，流域面积约16200km²。流域呈羽状水系，
~ ~
地势西北高、东南低，干流自西北流向东南。流域内地貌类型涉及河谷平原、黄土丘陵、高山
牧场以及中、高山地，最低点海拔1612m，最高点海拔4860m，具有较为明显的垂直差异性。湟
水流域总面积仅占青海省的2.24%，但人口、工农业生产总值等均占整个青海省的50%以上，
是青海省人口密度最大、经济最为发达地区。与此同时，湟水流域年降水量300 500mm，位于
~
干旱半干旱气候地区，是中国400mm等降雨量生态过渡地带与典型的农牧交错带。一方面，
这个区域处在不断集聚人口经济的过程中，由此产生的城镇扩张与耕地保护之间的矛盾重
重。另一方面，农牧交错带本身生态极为脆弱，是中国政府主导的退耕还林、退耕还草的政
策区域，耕地生产与生态脆弱区的保护之间的矛盾也极为尖锐。
[0031] 本发明实施例时空一致性双向检测支持下的基于大数据的长时序国土空间地域功能结构变化检测方法，包括以下步骤，如图2所示：
[0032] 步骤1、数据准备。
[0033] 基于Google Earth Engine云平台的Landsat影像数据集的收集和预处理。借助Google Earth Engine云平台收集从1987‑2020年间的云量少于15%的Landsat数据，并完成
影像拼接、裁剪等处理。
影像拼接、裁剪等处理。
[0034] 步骤2、样本点选取。
[0035] 确定地物类型和分类体系，通过实地采样或借助高分辨率影像选取2020年不同地物的样本点，并以此年为基准年，切换2019年参考影像，逐点比对，对地物类别发生改变的
样本点进行修改，也可以视情况进行删除或新增样本点，从而得到2019年的样本点，以此类
推，直到完成全部年份的样本点选取。
样本点进行修改，也可以视情况进行删除或新增样本点，从而得到2019年的样本点，以此类
推，直到完成全部年份的样本点选取。
[0036] 本步骤中，功能用地分类。以开发和保护功能两个维度为指向，将建设用地、农业用地作为开发指向的功能类型，将生态用地作为保护功能指向的功能类型。其中，农业用地
主要指耕地，生态用地包括林地、草地、灌木、水体、未利用地等。
主要指耕地，生态用地包括林地、草地、灌木、水体、未利用地等。
[0037] 步骤3、影像初始分类。
[0038] 将每幅影像的样本点按比例（80%的样本点作为训练样本点，20%作为验证样本点）随机划分为训练集和验证集，用历年训练集对分类器进行训练，训练完成后，使用分类器对
对应年份的影像进行分类，获得初始分类结果，分类结果具有相同的分辨率。
对应年份的影像进行分类，获得初始分类结果，分类结果具有相同的分辨率。
[0039] 本实施例中，使用随机森林分类器对遥感影像进行分类。
[0040] 训练方法如下：引入以下步骤，以获得较好的初始分类结果。
[0041] （1）训练样本点分布优化，使用GEE内置随机函数对训练样本点进行10次随机位置分布，取精度最高的结果，作为最终训练样本；
[0042] （2）引入数种纹理特征，试验数种尺寸窗口，取精度最高的结果，作为最终纹理特征；
[0043] （3）引入数种气候要素特征，取精度最高的结果，作为最终气候特征；
[0044] （4）进行使用步骤（1）的最终训练样本、步骤（2）的最终纹理特征和步骤（3）的最终气候特征，对分类器进行训练。
[0045] 使用训练完成的分类器对对应年份的影像进行分类，获得初始分类结果，从GEE云平台下载到本地，进行后续处理。
[0046] 步骤4、分类预处理。
[0047] 初步分类结果会存在“椒盐”现象，使用ArcGIS软件中众数滤波工具对其进行迭代处理，去除碎小图斑，直至结果不再变化，保持相对稳定。
[0048] 步骤5、构建分类结果时序数据集。
[0049] 使用ArcGIS软件中波段合成工具对预处理后历年分类结果进行合成，1987为第一波段，2020为最后一个波段，构造1987‑2020时间序列数据集。可以理解为将所有年份的影
像分类结果叠加，构建分类结果的时间序列栅格数据集，时间序列即为取该栅格数据集中
的某一行某一列像素的多年集合。
像分类结果叠加，构建分类结果的时间序列栅格数据集，时间序列即为取该栅格数据集中
的某一行某一列像素的多年集合。
[0050] 步骤6：双向时空一致性检测。
[0051] 在每个时间序列中分别构造两个种子窗口Ws和两个滑动检测窗口Wd，两个种子窗口Ws的初始位置位于时间序列的首末端，两个滑动检测窗口Wd的初始位置位于紧邻于对应
的种子窗口Ws，按照以下步骤进行处理，如图3所示：
的种子窗口Ws，按照以下步骤进行处理，如图3所示：
[0052] a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；
[0053] b）、将滑动窗口置于种子窗口紧邻的内侧；
[0054] c）、统计滑动窗口内的主导地类Kd，窗口内某地类出现频率高于60%则称其为该窗口的主导地类Kd；若窗口内所有地类出现频率均低于60%，则该窗口没有主导地类；
[0055] d）、判断种子窗口Ws内地物类型Ks与紧邻的滑动检测窗口Wd内的主导地类Kd是否一致，若一致，则将滑动检测窗口Wd内首尾主导地类栅格之间所有地类均置为Kd，将种子窗
口位置移动到滑动窗口内最后出现主导地类Kd的栅格位置，转至步骤b)直到两个滑动检测
窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口
相遇；
口位置移动到滑动窗口内最后出现主导地类Kd的栅格位置，转至步骤b)直到两个滑动检测
窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口
相遇；
[0056] e）、当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口地类一致的像素向左移动，右侧同理，若还有其他地类像素，则将
其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年。
其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年。
[0057] 测试序列：[2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2]
[0058] 输出序列：[2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
[0059] 其中，1代表耕地，2代表灌木，可见在测试序列中,某一地块的耕地和灌木存在频繁的地物类型转换，这与实际情况是不符的，而经过双向时空检测滤波处理后，重新梳理了
地物类型转换趋势，消除了某些年份孤立存在的地物点，确定了地物转换的年份。
地物类型转换趋势，消除了某些年份孤立存在的地物点，确定了地物转换的年份。
[0060] 步骤7：获得修正后的分类结果。
[0061] 将时间序列栅格数据集按顺序拆分，同一年份的栅格单元构成该年份的分类结果，最终获得所有年份修正后的分类结果，如图4 、图5所示。
[0062] 精度评价。
[0063] 使用验证样本对结果进行精度评价，利用混淆矩阵计算生产者精度和用户精度对各类地物进行评价，并算出总体精度和kappa系数。
[0064] 选取完整率和正确率两个指标，其中完整率为某地类分类所得像元数与地类实际像元总数之比，对应的是漏分；正确率为正确分类像元数与地类实际像元总数之比，对应的
是错分。
是错分。
[0065] 研究发现，初步分类结果的总体准确率为76.38%‑89.07%，总体平均准确率为82.64%；双向检测结果的总准确率为84.20%‑91.74%，总平均准确率为88.25%。很明显，双向
检测方法的总体变化检测精度明显优于初始分类检测,总体平均精度提高了5.61%，如图6
所示。不同年份的精度提高到不同程度，几年来整体分类准确率显著提高。与初始分类结果
相比，精度提高0.92% 10.22%，2006、2009、1996年精度提高最显著，增加10.22%、10.00%和
~
9.17%。
检测方法的总体变化检测精度明显优于初始分类检测,总体平均精度提高了5.61%，如图6
所示。不同年份的精度提高到不同程度，几年来整体分类准确率显著提高。与初始分类结果
相比，精度提高0.92% 10.22%，2006、2009、1996年精度提高最显著，增加10.22%、10.00%和
~
9.17%。
[0066] 长时序功能制图与变化检测。
[0067] 1）、长时序功能用地分类——以建设用地、农业用地、生态用地为一级分类进行功能类型合并，计算不同时期功能用地的面积及比重，进行长时序功能用地制图；
[0068] 2）、长时序地域功能结构测算——以建设用地、农业用地、生态用地比例关系为基本参数，选取主要交通干线为轴线，以主要轴线和核心城市为中心设定不同缓冲区，统计不
同缓冲区内建设用地、农业用地、生态用地比例关系，测算区域不同时期点‑轴结构和核心‑
边缘等地域功能结构。
同缓冲区内建设用地、农业用地、生态用地比例关系，测算区域不同时期点‑轴结构和核心‑
边缘等地域功能结构。
[0069] 本实施例中，按城镇建设用地、农业用地和生态用地三类功能用地进行类型归类。通过汇总分析可知，1987年、2000年、2010年、2020年三类功能用地比重分别为：1.23%、
27.76%、71.01%；2.75%、27.55%、69.70%；3.84%、24.71%、71.45%；6.38%、23.47%、70.15%。尤
2
以农业用地转变为建设用地最为显著，1987‑2020年转变为737.21km ，占农业用地总流失
面积的39.54%。其中，农业用地与生态用地在1987‑2020年间存在多次的转变关系，受退耕
还林等政策影响，1987‑2000、2000‑2010、2010‑2000三个时期，农业用地流向生态用地面积
2 2 2
分别为：731.91km 、922.06km 、554.29km ；而生态用地流向农业用地面积分别为
2 2 2
905.37km、557.86km 、632.69km 。，图7展示了三个时期各用地类型面积的流入和流出情
况。
27.76%、71.01%；2.75%、27.55%、69.70%；3.84%、24.71%、71.45%；6.38%、23.47%、70.15%。尤
2
以农业用地转变为建设用地最为显著，1987‑2020年转变为737.21km ，占农业用地总流失
面积的39.54%。其中，农业用地与生态用地在1987‑2020年间存在多次的转变关系，受退耕
还林等政策影响，1987‑2000、2000‑2010、2010‑2000三个时期，农业用地流向生态用地面积
2 2 2
分别为：731.91km 、922.06km 、554.29km ；而生态用地流向农业用地面积分别为
2 2 2
905.37km、557.86km 、632.69km 。，图7展示了三个时期各用地类型面积的流入和流出情
况。
[0070] 除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。