CN115909044B - 一种国土空间结构时空演变模式挖掘方法 - Google Patents

Abstract

本发明公开一种国土空间结构时空演变模式挖掘方法，方法如下：首先，借助于Google Earth Engine多源地理大数据平台，利用随机森林与时空立方体一致性方法，实现国土空间地域功能长时序精确分类；其次，借助于时空立方体网格化分割与长时序功能参数计算，实现网格化单元的空间聚类，并利用时空变化点检测识别有效变化点。最后，依据关键变化点所处位置，结合聚类结果进行重新分类，实现时空演变不同模式识别与分类，解析不同模式的应用潜力。该方法提高了国土空间地域功能分类与空间结构时空演变模式识别的精度，尤其适用于大范围、长时序、人类活动引发的区域国土空间地域功能与空间结构演变识别与分类，能够直接应用于国土空间规划与空间治理的辅助决策。

Description

一种国土空间结构时空演变模式挖掘方法

技术领域

本发明涉及空间结构演变技术领域，具体为一种大数据支持下的国土空间结构时空演变模式挖掘方法。

背景技术

国土空间结构是指不同开发与保护功能指向的国土空间面积、比例结构的空间组织特征，与自然资源空间分布、人类生产生活活动空间组织密切相关。以往受制于数据获取方式、数据时长、空间范围的约束，难以对长时序、大范围的国土空间结构时空演变模式进行识别与挖掘。随着多源遥感数据的长期积累以及处理技术的发展，利用大范围、高频率、长时序的多源遥感大数据，挖掘国土空间地域功能与结构演变模式，揭示人类生产生活活动空间组织与地表作用过程，成为可能。

现有技术的不足：

以往传统国土空间或土地利用功能时空演变模式研究，主要是借助于遥感影像分类数据或土地调查数据，从数量、比例关系、空间集聚等方面进行识别与挖掘。由于数据量大、获取难度大，以往研究多采用5-10年的时间间隔进行分析，通过计算景观结构等指数，揭示不同功能用地的时空演变特征。尽管这样能简化数据获取和处理过程，但也会遗漏掉可能存在的突变节点，忽略一些有价值的国土空间结构演变模式细节。然而，当进行大尺度长时序、高频率的遥感长时序影像分析时，受制于处理平台的计算能力和储存设备，且需要付出巨大的时间成本，传统遥感手段显然不是理想平台。

发明内容

本发明的目的在于提供一种国土空间结构时空演变模式挖掘方法，以解决上述背景技术中提出的问题。

一种国土空间结构时空演变模式挖掘方法，本挖掘方法包括以下步骤：S1、国土空间地域功能的分类分级；S2、国土空间地域功能长时序分类；S3、基于哈希表的像素级时空立方体构建；S4、基于像素级时空立方体的时间一致性精校正；S5、基于像素级时空立方体的空间一致性精校正；S6、时空立方体的空间网格化分割与地域功能参数计算；S7、网格化单元长时序地域功能参数的自动聚类；S8、长时序功能参数的变化点识别；S9、时空演变模式识别与分类；S10、不同时空演变模式的特征分析与应用，其中：

S1、国土空间地域功能的分类分级：以地域功能生成机理为基点，从人类生产、生活、生态三类基本功能为划分依据，考虑土地利用方式与调查特征，将国土空间划分为建设用地、农业用地、生态用地三个一级分类；

S2、国土空间地域功能长时序分类：借助GEE平台，采用随机森林算法对国土空间地域功能进行长时序分类；

S3、基于哈希表的像素级时空立方体构建：以像素为基本单元，根据像素级时间序列排序类型，对研究区域所有像素进行类型标记，哈希表每个唯一的key值对应存储不同序列类型；

S4、基于像素级时空立方体的时间一致性精校正：构建时空立方体的滑动窗口W_m和种子窗口W_s，滑动窗口W_m的初始位置位于紧邻于对应的种子窗口W_s的外侧，对研究区国土空间地域功能长时序分类数据进行时间一致性校正；

S5、基于像素级时空立方体的空间一致性精校正：将每一个像素级时间序列作为一个处理单元，利用k*k滑动窗口，采用众数滤波方法，对窗口中心像元时间序列分类数据进行整体性纠正，直至对所有像素级时空立方体完成精校正；

S6、时空立方体的空间网格化分割与地域功能参数计算：参考长时序地域功能分类数据的空间分辨率r与研究区范围大小S，采用网格化单元进行分割，网格化单元k，应满足k介于100r～1000r之间，k²介于S/10000～S/1000，k取值过程采用取值区间的二分法进行试验；

以不同功能类型比例关系为地域功能参数，逐年计算第t年第m个网格化单元内第q类国土空间的地域功能参数P_t ^m,q：

其中，

表示对应地域功能类型所占面积，s为网格化单元面积；

S7、网格化单元长时序地域功能参数的自动聚类：根据功能参数时间序列的特征，采用欧氏距离作为相似度，利用K-medoids聚类方法，对所有网格化单元地域功能参数长时序数据进行自动聚类；

S8、长时序地域功能参数的变化点识别：针对整个研究区域和每一个网格化单元，识别变化点；

S9、时空演变模式识别与分类：对所有网格化单元的时空演变模式进行识别与分类；

S10、不同时空演变模式的特征分析与应用：结合国土空间类型的具体分析，解析不同网格化单元类型所代表的具体现实含义，并分析其在国土空间开发保护格局、规划与决策中所起的作用。

作为本发明的进一步改进，本方法中步骤S2包括有以下步骤：

a1、数据搜集与预处理：利用地理云计算服务平台收集指定时间范围和研究范围的所有云量少于15％的Landsat数据，以及VIIRS夜间灯光数据、SRTM数字高程数据，进行影像拼接与裁剪的处理；

a2、样本点选取与分布优化：借助于Google Earth高分辨率历史影像，选取指定时间范围内每一年的不同功能地物的样本点，将其中80％的样本点作为训练样本点，20％作为验证样本点，使用地理云计算服务平台内置随机函数对训练样本点进行20次随机位置分布，取精度最高的结果，作为最终训练样本；

a3、特征向量计算与遴选：根据研究区域的特征，计算归一化植被指数(NDVI)、归一化建筑指数(NDBI)、归一化差水指数(NDWI)特征变量，引入GEE平台可计算的所有纹理特征，与VIIRS夜间灯光数据、数字高程数据(SRTM)构成不同特征变量的组合方式，对比不同组合方式的精度结果，筛选出最优的特征向量；

a4、基于随机森林分类器的长时序分类：利用最终训练样本、最优特征变量组合对随机森林分类器进行训练，完成训练后，逐年对相应年份的影像进行分类；

a5、初始精度评价：利用验证样本，依次对每一年国土空间地域功能分类结果进行初始精度评价，计算总体精度和kappa系数，将多年平均总体精度和kappa系数作为初始精度评价指标。

作为本发明的进一步改进，本方法中步骤S4包括有以下步骤：

b1、针对每一个key值存储的时间序列，分别构造两个种子窗口W_s和两个滑动检测窗口W_m，两个种子窗口W_s的初始位置位于时间序列的首末端，设置滑动时间窗口步长n，n为大于3的自然数；

b2、将滑动窗口W_m置于种子窗口W_s外侧；

b3、计算滑动窗口W_m内的主导地类f_m；

b4、判断种子窗口W_s内地物类型f_s与紧邻的滑动检测窗口W_m内的主导功能地类f_m是否一致，若一致，则将滑动窗口W_m内首尾主导地类栅格之间所有地类均置为f_s，并将种子窗口位置移动到滑动窗口内最后出现主导地类f_s的栅格位置，转至步骤b2直至时序首或尾年份，处理完哈希表存储的所有时空立方体；若不一致，则将种子窗口向内移动一位，转至步骤b2直至时序首或尾年份，处理完哈希表存储的所有时空立方体。

作为本发明的进一步改进，本方法中步骤S7相似度采用如下公式计算：

其中，dist(i,j)表示第i、j网格化单元的相似度，n表示时间序列长度，t表示第t个时间节点，

分别表示第t个时间节点第i、j个格网的地域功能参数值。

作为本发明的进一步改进，方法中步骤S8识别变化点包括以下步骤：

c1、时间序列地域功能参数的平稳性分析，按照以下公式，将第t个时间节点(t∈(T₀,T_n))的地域功能参数

与其前一时刻的功能参数/>

相减，然后对差值取对数，去除时间序列的趋势性，转变为平稳性时间序列；

取对数时加1是为了避免当面积差不为正数时取对数产生缺失值；

c2、变化节点数量k确定，根据预设的时间演变类型数目，对应设置k个变化节点，分别计算k+1时间序列不同分段的标准差，如果不同分段标准差在阈值以内，k设置为0；

c3、当时间序列中有k个变化点(k∈[1,n-1])，它们的位置集合为τ_0:k＝(τ₀,…,τ_k+1)，其中τ₀＝0，τ_k+1＝n，k个变化点会将序列分为k+1段，其中第i段为

对包含k个变化点的序列/>

进行分割，其最小成本记为/>

则/>

可表示为：

其中

是分割点的成本函数，采用负对数似然的两倍表示。

在子序列

中(t∈(1,n))，最优分割/>

按如下方法被递归地求解：

其中τ_k∈{k,…,t-1}，从最优分割中提取变化点的位置通过对序列的从后向前递归实现，首先找到

中最后一个变化点的位置，如果该位置在时刻τ上，那么就在/>

中寻找下一个变化点，由此重复，直到提取出指定数量的变化点，如果相邻的t、t+1时刻被识别为变化点，则留取分割成本较低的变化点作为该序列中唯一的变化点。

作为本发明的进一步改进，本方法中步骤S9包括有以下步骤：

d1、确定时间分段，以整个研究区域的变化点为划分标准，将时间序列划分为不同阶段；

d2、网格化单元关键节点确定与归类：依次计算变化点前后分段的功能参数均值，取变化点前后均值绝对值最大的节点，作为关键节点。根据变化关键节点所处的时间位置，将网格化单元归类到整体研究区域对应的时间分段；

d3、将自动聚类结果与关键变化节点分段结果进行集成，结合研究区域国土空间结构的总体特征。

作为本发明的进一步改进，本方法步骤S1考虑遥感影像分类特点进行二级分类，将生态用地进一步划分为林地、草地、水体、裸地，建设用地为不透水面，农业用地为耕地。

作为本发明的进一步改进，本方法中地理云计算服务平台为GEE平台。

作为本发明的进一步改进，本方法步骤b3中滑动窗口W_m内某地类出现频率高于r则称该地类为该滑动窗口W_m的主导功能地类f_m，若滑动窗口内所有地类出现频率均低于r，则该滑动窗口没有主导功能地类，不做接下来的处理，r取值在0.7-1.0。

作为本发明的进一步改进，本方法步骤d3中将所有网格化单元划分为稳定高值类、前期启动高值类、后期启动高值类、前期启动中值类、后期启动中值类、前期启动低值类、后期启动低值类、稳定低值类等多种不同类型。

与现有技术相比，本发明的有益效果是：

本发明采用开源地理信息大数据平台，利用多源大数据实现对国土空间地域功能的长时序分类，通过像素级时空立方体构建与网格化分割，在此基础上对网格化单元进行聚类与变化点识别，实现长时序国土空间结构演变模式的识别、分类与挖掘，适用于大范围、长时序、高频率的国土空间地域功能与空间结构识别与分类，尤其是适用于快速扩张的城市群或城市化连绵区的国土空间结构演变模式挖掘。

附图说明

图1为本发明一种国土空间结构时空演变模式挖掘方法实施例研究区概况图；

图2为本发明一种国土空间结构时空演变模式挖掘方法的总体技术流程图；

图3为1990-2020年国土空间地域功能分类结果图；

图4为时空一致性校正前后的精度对比图；

图5为开发强度时间序列聚类结果图；

图6为网格化单元建设开发功能关键节点分析；

图7为环渤海地区建设用地时空演变模式挖掘结果。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

请参阅图1-7，本发明提供如下技术方案：一种国土空间结构时空演变模式挖掘方法，本方法采用的实验数据包括：Landsat系列历史影像数据，包括TM/ETM+/OLI等不同类型，辅助分类的数据还包括VIIRS夜间灯光数据、数字高程数据(SRTM)及降水量、气温、积温等气候数据；

其中，VIIRS夜间灯光数据(Nighttime Day/Night Band Composites Version 1)来自可见光红外成像辐射计套件(VIIRS)日/夜带(DNB)的夜间数据的月平均辐射合成图像；

数字高程数据(SRTM)来自美国土地分布式活动档案中心，且已使用开放源数据(ASTER GDEM2，GMTED2010和The National Elevation Dataset)进行了填充处理；

气候数据使用FLDAS数据；

研究区选取中国环渤海地区，包括环绕着渤海全部及黄海的部分沿岸地区所组成的广大区域(如图1)；以京津冀为核心、以辽东半岛和山东半岛为两翼的环渤海经济区域，主要包括北京、天津、河北、山东、辽宁；主要范围为经度113.4°E～125.8°E，纬度34.3°N～43.5°N；全区陆域面积51.8万平方公里，具体步骤如下：

S1、国土空间地域功能的分类分级：

根据环渤海地区整体区域特征，将国土空间划分为生态用地、农业生产用地、建设用地3类一级基本功能，并进一步细分为林地、草地、水体、裸地、耕地、不透水面7类二级基本功能分类。

S2、国土空间地域功能长时序分类：

数据搜集与预处理，利用GEE平台收集1990-2019年覆盖环渤海地区的所有云量少于15％的Landsat TM/ETM+/OLI等不同类型数据，共筛选9242景影像，并完成影像拼接、裁剪处理，同时下载环渤海地区夜间灯光数据(Nighttime Day/Night Band CompositesVersion 1)、数字高程数据(SRTM)，按照以下步骤进行分类处理；

样本点选取与分布优化，借助于Google Earth高分辨率历史影像，选取2020年不同功能类型地物的样本点，并以此年为基准年，切换2019年参考影像，逐点比对，对地物类别发生改变的样本点进行增删修改，从而得到2019年的样本点，以此类推，直到完成全部年份的样本点选取，共选取样点37275个。将其中80％的样本点作为训练样本点，20％作为验证样本点，使用GEE内置随机函数对训练样本点进行20次随机位置分布，取精度最高的结果，作为最终训练样本；

特征向量计算与遴选，根据研究区域的特征，计算归一化植被指数(NDVI)、归一化建筑指数(NDBI)、归一化差水指数(NDWI)特征变量，引入GEE平台可计算的所有纹理特征，与VIIRS夜间灯光数据、数字高程数据(SRTM)构成不同特征变量的组合方式，对比不同组合方式的精度结果，筛选出最优的特征向量，最优特征向量包括：归一化植被指数(NDVI)、归一化建筑指数(NDBI)、归一化差水指数(NDWI)、修正归一化差水指数(MNDWI)，与VIIRS夜间灯光数据、数字高程数据(SRTM)；纹理特征包括差异值(diss)、惯性矩(inertia)、相关度(idm)，最佳的窗口尺寸范围大致在3-5，每幅影像的最佳窗口尺寸不尽相同；

基于随机森林分类器的长时序分类，利用最终训练样本、最优特征变量组合对随机森林分类器进行训练，完成训练后，逐年对1990-2020年的影像进行分类，分类结果如图3所示；

经过计算发现，2020年环渤海地区耕地、林地、不透水面的面积比重分别为47.83％、23.18％、9.77％，是国土空间地域功能类型最多的三类，这符合广大沿海平原地区的主要特征，1990-2020年环渤海地区变化最为显著的土地类型是耕地、裸地与林地。其中，耕地向建设用地转变是土地利用变化的主要类型，共有3.95万平方千米的耕地转变为了建设用地，占全部变化面积的14.08％，由耕地转变而来的建设用地主要分布在北京、天津、沈阳、济南、石家庄等中心城市，淄博、潍坊、临沂等区域中心城市以及大连、青岛、烟台等沿海港口城市周边；

初始精度评价，利用验证样本，依次对每一年国土空间地域功能分类结果进行初始精度评价，计算总体精度和kappa系数，将多年平均总体精度和kappa系数作为初始精度评价指标，经过计算，1990-2020年间共29个时间节点(1999年和2012年无全覆盖影像)，分类时的平均总体精度为93.58％，kappa系数为91.71，2004年以后的影像分类结果均大于93.58％，且总体上越靠近2020年的年份影像分类结果越高。

S3、基于哈希表的像素级时空立方体构建：

以像素为基本单元，根据像素级时间序列排序类型，对研究区域616189127个像素进行类型标记，哈希表每个唯一的key值对应存储不同序列类型，共记录了389846个不同序列类型。

S4、基于像素级时空立方体的时间一致性精校正：

构建时空立方体的滑动窗口W_m和种子窗口W_s，滑动窗口W_m的初始位置位于紧邻于对应的种子窗口W_s的外侧，对研究区国土空间地域功能长时序分类数据进行时间与空间一致性校正，具体步骤如下：

针对每一个key值存储的时间序列，分别构造两个种子窗口W_s和两个滑动检测窗口W_m，两个种子窗口W_s的初始位置位于1990年和2020年，将滑动时间窗口步长n设置为5；

将滑动窗口W_m置于种子窗口W_s外侧；

计算滑动窗口W_m内的主导地类f_m，滑动窗口W_m内某地类出现频率高于r则称该地类为该滑动窗口W_m的主导地类f_m；若滑动窗口内所有地类出现频率均低于r，则该滑动窗口没有主导地类，不做接下来的处理；r取值设定为0.8；

判断种子窗口W_s内地物类型f_s与紧邻的滑动检测窗口W_m内的主导地类f_m是否一致，若一致，则将滑动窗口W_m内首尾主导地类栅格之间所有地类均置为f_s，并将种子窗口位置移动到滑动窗口内最后出现主导地类f_s的栅格位置，转至步骤b2直至两个滑动窗口相遇，处理完哈希表存储的所有时空立方体；若不一致，则将种子窗口向内移动一位，转至步骤b2直至两个滑动窗口相遇，处理完哈希表存储的所有时空立方体。

S5、基于像素级时空立方体的空间一致性精校正：

将每一个像素级时间序列作为一个处理单元，利用3*3滑动窗口，采用众数滤波方法，对窗口中心像元时间序列分类数据进行整体性纠正，直至对所有像素级时空立方体完成精校正。

经过步骤S4和步骤S5的精校正，校正后的分类平均总体精度达到95.81％，大幅提高2.23％(图4)，这种方法在一些分类精度较低的年份得到了较好的验证。其中有5年的分类精度提高了4％以上，在超过1/3的时间点上，分类精度提高了3％以上；

S6、时空立方体的空间网格化自动分割与功能参数计算：

针对建设用地，对时空立方体的空间网格化分割与功能参数计算，按以下步骤进行计算：

采用10km网格化单元进行分割，将环渤海地区像素级时空立方体分割为5996个网格化单元；

以建设用地占网格化单元的比重为功能参数，即开发强度，逐年计算1990-2020年开发强度数值。

经过计算发现，共分割出5996个网格化单元，在2020年，开发强度在20％以上有930个网格化单元，在5-10％之间有843个网格化单元，在1％以下有4050个网格化单元，占比为67.55％。

S7、网格化单元长时序开发强度的自动聚类：

使用开源R语言跨平台库TSclust包中的diss.EUCL()函数，采用欧氏距离进行相似度量；使用包中的pam()函数，采用K-medoids聚类方法对所有网格化单元开发强度长时序数据进行自动聚类，结果如图5；

聚类结果表明，聚为5类效果最佳，不同聚类间时间序列数值能够有比较好的区分，高值类、中高值类、中值类、中低值类、低值类分别有56个、185个、539个、1245个、3971个网格化单元，占比分别为0.93％、3.09％、8.99％、20.76％、66.23％。

S8、网格化单元长时序功能参数的变化点识别：

时间序列功能参数的平稳性分析，从1991年开始，依次将后一年的开发强度减去前一年的开发强度，并按照以下公式计算，实现时间序列的平稳性分析：

取对数时加1是为了避免当面积差不为正数时取对数产生缺失值；

将变化节点数量k确定为2，通过计算三个时间分段的标准差，将标准差小于0.01的网格化单元的k值设置为0；

时间变化点检测，使用开源R语言changepoint包中的cpt.var()函数，根据子分段与整条时间序列之间的标准差的大小关系，在具有指定数量的变化点的所有可能分段中查找分割成本最低的分段；

在设定固定两个变化点时，该方法会把每条时间序列都划分为三个分段，结果如图6所示，可见，每条子分段与全局具有相同的全局平均值，但标准差会在变化点位置发生变化。此时函数的应用方法为：

cpt.var(data,method＝”SegNeigh”,Q＝3)；

经过计算，3971个网格化单元无明显的变化点，2025个网格化单元有2个变化点，且变化点时间距离在5年以上。

S9、时空演变模式识别与分类：

按以下步骤对所有网格化单元的时空演变模式进行识别与分类：

根据环渤海整个区域不透水面的变化，将2003、2014两个节点作为划分节点，划分为阶段I、II、III三个阶段(图6)，对应前期、中期和后期三个时间段；

网格化单元关键节点确定与归类：依次计算变化点前后分段的功能参数均值，取变化点前后均值绝对值最大的节点，作为关键节点，根据变化关键节点所处的时间位置，将网格化单元归类到整体研究区域对应的时间分段，由于关键节点位于前期I阶段较少，因此把I阶段和II阶段合并，统一称为前期阶段；

将空间自动聚类结果与关键变化节点分段结果进行集成，结合研究区域国土空间结构的总体特征，将所有网格化单元划分为稳定高值类(56个)、前期启动中高值类(108个)、后期启动中高值类(77个)、前期启动中值类(271个)、后期启动中值类(268个)、前期启动中低值类(614个)、后期启动中低值类(631个)、稳定低值类(3971个)8个类型。

S10、不同时空演变模式的特征分析与应用：

结合环渤海地区具体情况，稳定高值类主要指北京、天津、沈阳、大连、石家庄、济南、青岛等主要大城市的老城区，前期启动中高值类、后期启动中高值类是主要大城市以及中等城市扩张的新城新区，前期启动中值类、后期启动中值类主要是等级较小城市与县城、较大规模中心镇及其扩张区，前期启动中低类、后期启动中低类主要指中心镇以及规模较大的农村居民点，稳定低值类主要指人口密度较低的广大农业主产区与生态地区。通过时空演变模式挖掘结果的分析(图7)，可以发现，环渤海地区建设用地扩张呈现出以省会城市和沿海大城市为核心，逐步向沿海和辽东半岛、山东半岛两翼中小城市扩散的时空格局，建设用地扩张规模和速率也呈现出沿中心城市和主要轴线依次降低的空间梯度变化。

本发明通过采用开源地理信息大数据平台，利用多源大数据实现对国土空间地域功能的长时序分类，通过像素级时空立方体构建与网格化分割，在此基础上对网格化单元进行聚类与变化点识别，实现长时序国土空间结构演变模式的识别、分类与挖掘，适用于大范围、长时序、高频率的国土空间地域功能和空间结构识别与分类，尤其是适用于人类活动主导的国土空间结构演变模式挖掘。

需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性地包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (8)

1.一种国土空间结构时空演变模式挖掘方法，其特征在于：本挖掘方法包括以下步骤：S1、国土空间地域功能的分类分级；S2、国土空间地域功能长时序分类；S3、基于哈希表的像素级时空立方体构建；S4、基于像素级时空立方体的时间一致性精校正；S5、基于像素级时空立方体的空间一致性精校正；S6、时空立方体的空间网格化自动分割与地域功能参数计算；S7、网格化单元长时序地域功能参数的自动聚类；S8、长时序功能参数的变化点识别；S9、时空演变模式识别与分类；S10、不同时空演变模式的特征分析与应用，其中：

S1、国土空间地域功能的分类分级：以地域功能生成机理为基点，从人类生产、生活、生态三类基本功能为划分依据，考虑土地利用方式与调查特征，将国土空间划分为建设用地、农业用地、生态用地三个一级分类；

S2、国土空间地域功能长时序分类：借助GEE平台，采用随机森林算法对国土空间地域功能进行长时序分类；

S3、基于哈希表的像素级时空立方体构建：以像素为基本单元，根据像素级时间序列排序类型，对研究区域所有像素进行类型标记，哈希表每个唯一的key值对应存储不同序列类型；

S4、基于像素级时空立方体的时间一致性精校正：构建时空立方体的滑动窗口W_m和种子窗口W_s，滑动窗口W_m的初始位置位于紧邻于对应的种子窗口W_s的外侧，对研究区国土空间地域功能长时序分类数据进行时间一致性校正；

S5、基于像素级时空立方体的空间一致性精校正：将每一个像素级时间序列作为一个处理单元，利用k*k滑动窗口，采用众数滤波方法，对窗口中心像元时间序列分类数据进行整体性纠正，直至对所有像素级时空立方体完成精校正；

S6、时空立方体的空间网格化自动分割与地域功能参数计算：参考长时序地域功能分类数据的空间分辨率r与研究区范围大小S，采用网格化单元进行分割，网格化单元k，应满足k介于100r～1000r之间，k²介于S/10000～S/1000，k取值过程采用取值区间的二分法进行试验；

以不同功能类型比例关系为地域功能参数，逐年计算第t年第m个网格化单元内第q类国土空间的地域功能参数P_t ^m,q：

其中，

表示对应地域功能类型所占面积，s为网格化单元面积；

S7、网格化单元长时序地域功能参数的自动聚类：网格化单元长时序地域功能参数的自动聚类，根据地域功能参数时间序列的特征，采用欧氏距离作为相似度，利用K-medoids聚类方法，对所有网格化单元地域功能参数长时序数据进行自动聚类；

S8、长时序地域功能参数的变化点识别：针对整个研究区域和每一个网格化单元，识别变化点；步骤S8识别变化点包括以下步骤：

c1、时间序列地域功能参数的平稳性分析，按照以下公式，将第t个时间节点(t∈(T₀,T_n))的地域功能参数

与其前一时刻的地域功能参数/>

相减，然后对差值取对数，去除时间序列的趋势性，转变为平稳性时间序列：

取对数时加1是为了避免当面积差不为正数时取对数产生缺失值；

c2、变化节点数量k确定，根据预设的时间演变类型数目，对应设置k个变化节点，分别计算k+1时间序列不同分段的标准差，如果不同分段标准差在阈值以内，k设置为0；

c3、当时间序列中有k个变化点(k∈[1,n-1])，它们的位置集合为τ_0:k＝(τ₀,…,τ_k+1)，其中τ₀＝0，τ_k+1＝n，k个变化点会将序列分为k+1段，其中第i段为

对包含k个变化点的序列/>

进行分割，其最小成本记为/>

则/>

可表示为：

其中

是分割点的成本函数，采用负对数似然的两倍表示；

在子序列

中(t∈(1,n))，最优分割/>

按如下方法被递归地求解：

其中τ_k∈{k,…,t-1}，从最优分割中提取变化点的位置通过对序列的从后向前递归实现，首先找到

中最后一个变化点的位置，如果该位置在时刻τ上，那么就在/>

中寻找下一个变化点，由此重复，直到提取出指定数量的变化点，如果相邻的t、t+1时刻被识别为变化节点，则留取分割成本较低的变化节点作为该序列中唯一的变化点；

S9、时空演变模式识别与分类：对所有网格化单元的时空演变模式进行识别与分类；步骤S9包括有以下步骤：

d1、确定时间分段，以整个研究区域的变化点为划分标准，将时间序列划分为不同阶段；

d2、网格化单元关键节点确定与归类：依次计算变化点前后分段的功能参数均值，取变化点前后均值绝对值最大的节点，作为关键节点，根据变化关键节点所处的时间位置，将网格化单元归类到整体研究区域对应的时间分段；

d3、将自动聚类结果与关键变化节点分段结果进行集成，结合研究区域国土空间结构的总体特征；

S10、不同时空演变模式的特征分析与应用：结合国土空间类型的具体分析，解析不同网格化单元类型所代表的具体现实含义，并分析其在国土空间开发保护格局、规划与决策中所起的作用。

2.根据权利要求1所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法中步骤S2包括有以下步骤：

a1、数据搜集与预处理：利用地理云计算服务平台收集指定时间范围和研究范围的所有云量少于15％的Landsat数据，以及VIIRS夜间灯光数据、SRTM数字高程数据，进行影像拼接与裁剪的处理；

a2、样本点选取与分布优化：借助于Google Earth高分辨率历史影像，选取指定时间范围内每一年的不同功能地物的样本点，将其中80％的样本点作为训练样本点，20％作为验证样本点，使用地理云计算服务平台内置随机函数对训练样本点进行20次随机位置分布，取精度最高的结果，作为最终训练样本；

a3、特征向量计算与遴选：根据研究区域的特征，计算归一化植被指数(NDVI)、归一化建筑指数(NDBI)、归一化差水指数(NDWI)特征变量，引入GEE平台可计算的所有纹理特征，与VIIRS夜间灯光数据、数字高程数据(SRTM)构成不同特征变量的组合方式，对比不同组合方式的精度结果，筛选出最优的特征向量；

a4、基于随机森林分类器的长时序分类：利用最终训练样本、最优特征变量组合对随机森林分类器进行训练，完成训练后，逐年对相应年份的影像进行分类；

a5、初始精度评价：利用验证样本，依次对每一年国土空间地域功能分类结果进行初始精度评价，计算总体精度和kappa系数，将多年平均总体精度和kappa系数作为初始精度评价指标。

3.根据权利要求1所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法中步骤S4包括有以下步骤：

b1、针对每一个key值存储的时间序列，分别构造两个种子窗口W_s和两个滑动检测窗口W_m，两个种子窗口W_s的初始位置位于时间序列的首末端，设置滑动时间窗口步长n，n为大于3的自然数；

b2、将滑动窗口W_m置于种子窗口W_s外侧；

b3、计算滑动窗口W_m内的主导功能地类f_m；

b4、判断种子窗口W_s内地物类型f_s与紧邻的滑动检测窗口W_m内的主导功能地类f_m是否一致，若一致，则将滑动窗口W_m内首尾主导地类栅格之间所有地类均置为f_s，并将种子窗口位置移动到滑动窗口内最后出现主导地类f_s的栅格位置，转至步骤b2直至时序首或尾年份，处理完哈希表存储的所有时空立方体；若不一致，则将种子窗口向内移动一位，转至步骤b2直至时序首或尾年份，处理完哈希表存储的所有时空立方体。

4.根据权利要求1所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法中步骤S7相似度采用如下公式计算：

其中，dist(i,j)表示第i、j网格化单元的相似度，n表示时间序列长度，t表示第t个时间节点，

分别表示第t个时间节点第i、j个网格的功能参数值。

5.根据权利要求1所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法步骤S1考虑遥感影像分类特点进行二级分类，将生态用地进一步划分为林地、草地、水体、裸地，建设用地为不透水面，农业用地为耕地。

6.根据权利要求2所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法中地理云计算服务平台为GEE平台。

7.根据权利要求3所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法步骤b3中滑动窗口W_m内某地类出现频率高于r则称该地类为该滑动窗口W_m的主导功能地类f_m，若滑动窗口内所有地类出现频率均低于r，则该滑动窗口没有主导功能地类，不做接下来的处理，r取值在0.7-1.0。

8.根据权利要求1所述的国土空间结构时空演变模式挖掘方法，其特征在于：本方法步骤d3中将所有网格化单元划分为稳定高值类、前期启动高值类、后期启动高值类、前期启动中值类、后期启动中值类、前期启动低值类、后期启动低值类、稳定低值类的不同类型。

Images