CN105184048A

CN105184048A - 一种基于多源数据的逐年土地利用空间数据的制备方法

Info

Publication number: CN105184048A
Application number: CN201510436380.8A
Authority: CN
Inventors: 崔耀平; 朱连奇; 周志民
Original assignee: Henan University
Current assignee: Guangdong Hongtu planning and Design Co., Ltd
Priority date: 2015-07-23
Filing date: 2015-07-23
Publication date: 2015-12-23
Anticipated expiration: 2035-07-23
Also published as: CN105184048B

Abstract

本发明公开了一种基于多源数据的逐年土地利用空间数据的制备方法。本发明使用两套数据：一是利用30m空间分辨率的Landsat？TM/ETM+高分辨率对地观测卫星影像，经过人工解译获取的土地利用变化数据（A-LUCC）；二是对地观测卫星(MODIS？500m空间分辨率)的土地利用数据(M-LUCC)。首先对两套不同源数据进行属性调整和数据融合，以A-LUCC的时间变化速率为约束值，参考M-LUCC的逐年变化率；其次，通过一种类归一化的线性差值计算方法获得逐年LUCC的时间变化；并以A-LUCC的空间变化为准，参考M-LUCC的各土地利用变化位置，完成逐年LUCC的空间定位。本发明克服了土地利用变化研究中，高精度的土地利用数据无法逐年提供的问题，本方法不仅可操作性较强，同时在一定尺度范围内有精度保证。

Description

一种基于多源数据的逐年土地利用空间数据的制备方法

技术领域

本发明涉及土地利用变化逐年数据的制备。主要运用非同源土地利用数据的时间和空间值的精度互补的思想，以A-LUCC空间分布为准，以同分析窗口下M-LUCC的时间变化为准，通过非同源数据的有效结合实现逐年的土地利用变化数据的制备，并且可以在一定的分析窗口和尺度单元下保证精度，从而为土地资源规划管理和发展评估提供科学的参考资料和决策依据。属于地球科学领域。

背景技术

土地利用/土地覆盖变化（LandUse/CoverChange,简称LUCC）时空格局变化可延伸为时间变化和空间变化这两个方面。时间变化研究的主线为：不同时期土地利用变化类型的面积动态指数和转移矩阵，以其为评价指标来量化和分级土地利用的变化速度、变化幅度和转移比例等。空间变化的研究主线为：以地理信息科学技术为研究工具，依据土地利用动态信息，划分区域土地利用类型的空间分布格局，进而明确土地利用类型随时间演化的动态转化过程。

对LUCC空间和时间过程变化的研究主要集中在以下两种类型：第一种是模型模拟研究，DeAlmeida等将元胞自动机模型应用到巴西城市Bauru的土地利用变化和扩展模拟。黎夏、邓祥征等人分别提出的城市CA模型和区域土地利用类型变化模型；第二种是年际变化研究，如刘纪远等以1980年代末、1995年、2000年和2005为时间截面，以5年为一研究阶段，揭示了20世纪80年代末以来中国土地利用的时空变化过程。

基于时间序列的土地利用变化宏观研究已经日臻成熟，受限于数据，对LUCC“年际变化”的长时间序列的过程研究却少之又少。而从时间尺度和空间尺度同时入手，以统计学理论和空间分析技术为工具，探索土地利用时空演化规律，克服土地利用变化研究中，高精度的土地利用数据无法逐年提供的问题，并在一定尺度范围内保证精度，已逐步成为土地利用变化研究的前沿思维。

发明内容

本发明旨在提出一种逐年土地利用空间数据的制备方法，以便可以实现长时间序列的逐年土地利用变化分析。该方法采用不同源的两种土地利用数据，结合其各自提供的有效信息，通过数据融合方式，实现逐年土地利用空间数据的制备。该方法可以在一定程度上替代模型模拟，是一种新的、简便易行的、可快速用于实践研究的方法。

为实现上述目的，本发明提供一种逐年土地利用空间数据的制备方法，其步骤依次为：

1）数据获取：

两种不同源数据：一是利用较高分辨率对地观测卫星影像30m空间分辨率的LandsatTM/ETM+，通过人工解译制备的土地利用变化数据集（A-LUCC）；二是500m空间分辨率MODIS的土地覆盖数据（M-LUCC）。

2）数据预处理：首先对多期A-LUCC数据进行地理空间分析算法中的矢量转栅格，并基于最近邻算法进行栅格数据重采样，即将空间分辨率重采样为500m，以保证和M-LUCC的空间分辨率一致。然后对A-LUCC和M-LUCC的土地利用类型，采用大类的分类标准，分为农田、林地、草地、湿地（含水域）、建设用地和其他未利用地等六大类（简称农、林、草、湿、城、荒）。

3）预设分析窗口

本发明通过均匀网格化来判定最小的分析窗口，即把整个研究区分成1×1,4×4,…等亚区作为分析窗口，在亚区内逐一判读2个研究时段内均有的土地利用转类，而满足此条件的最小网格亚区即为所需的最优分析窗口，分析窗口的上限是整个研究区。在取值区间内分析窗口越小精度越高。

4）时间变化的判定

对不同土地利用类型的时间变化速率，以A-LUCC的时间段变化值为限定值，参考M-LUCC的年际变化值来计算，这里称之为逐年类归一化线性拟合处理，具体推导过程如下：

A-LUCC数据中某土地利用转类i的j+n年（j=1,2,…;n≥1）变化总量ATR _{ni_sum}可用如下公式表示：

(1)

A-LUCC数据中某土地利用转类i的逐年变化量之和MTR _{ni_sum}并非一定等于ATR _{ni_sum}，此值的表达式为：

(2)

式中，MTR为M-LUCC对应转类i的对应年份j+n的变化速率；(MTR _j+n -MTR _j+n-1)、(MTR _j+1 -MTR _j)、…为逐年际间各土地利用转类的面积变化差值。则，如果MTR _{ni_sum}≠ATR _{ni_sum}时，需要采用如下方法处理：

使M-LUCC数据各土地利用类型变化的量等于两期A-LUCC数据各土地利用类型变化的量，即：

(3)

式中，j+m，j+n和j分别为研究时间段的任意年份、末年和首年对应的土地利用类型i变化值，其中1<m<n。则有此式可以推导出逐年LUCC类归一化的土地利用类型i对应的面积值为：

(4)

式中，TR _j+m为LUCC数据中某土地利用转类i的面积；AT和MTR分别为A-LUCC和M-LUCC中某一转类i的面积。

5）空间变化的判定

以A-LUCC空间分布为准，同分析窗口下M-LUCC逐年的同转类斑块分布位置作为参考，在变化速率值的限定下，利用条件语句判定空间转换位置。即：

如果该斑块的动态分布范围与在A-LUCC连续两期研究时段内(t和t+1)的同类斑块(Grid)动态分布范围一致，则为正确的动态转类空间定位，即任意栅格位置对应满足：Grid _t -Grid _t+1 =0；否则，空间范围不一致时，则定位该转类斑块到同窗口下A-LUCC的最邻近同转类斑块的所属空间，超出A-LUCC对应转类斑块空间的依照原M-LUCC的分布空间，进行同分析窗口下的质点所在空间趋势位置分为八个方向(东、南、西、北、东北、西北、东南、西南)，再随机分布定位：即，当质点偏向某一方向，就在某方向所在位置随机布置空间转类栅格类型。

其中，质心是描述地理目标空间分布最有用的单一量算量，它是保持目标均匀分布的平衡点。这里引入质心是旨在通过几何中心的基础上的加权运算来追踪土地类型的变化。任意多边形的质点求算可以通过公式表示：

任意边n边形顶点A _i (x _i ,y _i ),(i=1,2,…,n)依照逆时针排列，其面积为：

(5)

则，对应任意边n边形的质点可表达为：

(6)

6）精度验证和长时间序列时空分析

利用研究时段内的A-LUCC数据，随机选取研究区中某一年的数据，对应检验各类型的空间分布和面积是否一致，进而给出精度。或者更进一步，随机选取某年或某区域内人机交互解译的A-LUCC数据，以此对比检验同一年份的新制备数据，给出对应精度。

土地利用长时间序列时空分析是指根据制备出来的数据结果，看其空间变化情况，主要统计该研究时间段的各土地利用类型的数量。

本发明方法的优势和特点：

1）可间接获取连续时间序列的土地利用数据，对不连续的遥感数据进行精度分析和质量评价，从而提高土地利用变化模拟的效果。

2）将规则网格方法运用到土地利用转类的判读中，使分析对象更加地精细化，也使逐年信息的扩展更加地合理和准确。

3）通过中等空间分辨率和较高分辨率的遥感数据之间建立定量关系，互取其优点，为其他多源数据的融合处理提供了技术参照。

附图说明

图1是制备数据的分析窗口判读示意图

图2是制备逐年土地利用空间数据的方法流程图

具体实施方式

下面结合发明内容和附图1和2对本发明的原理及具体实施方式作进一步阐述。包括如下步骤：

1.数据获取

利用较高分辨率对地观测卫星影像30m空间分辨率的LandsatTM/ETM+，通过人工解译制备的土地利用变化数据集（A-LUCC）；利用500m空间分辨率MODIS的土地覆盖数据（M-LUCC）。

2.数据预处理

(1)A-LUCC人机交互解译。首先对多期A-LUCC数据进行地理空间分析算法中的矢量转栅格，并基于最近邻算法进行栅格数据重采样，即将空间分辨率重采样为500m，以保证和M-LUCC的空间分辨率一致。利用屏幕数字化方法以研究区内对相同的视图比例人机交互解译放大的较高空间分辨率影像，同时开展野外考察，确定分类标准，并更正分类结果，以满足本研究的需要。

(2)对M-LUCC土地覆盖数据对应A-LUCC的解译分类标准，最终统一将两套数据的土地利用类型分为农、林、草、湿、城、荒六类。A-LUCC的解译结果即为农、林、草、湿、城、荒六类，因此不涉及栅格属性的调整；而对M-LUCC，其土地覆盖类型（表1）需要，根据表2统一属性为：农、林、草、湿、城、荒六类。

表1五种土地覆盖分类方案

代码	IGBP	UMD	LAI/fPAR	NPP	PFT
						0	水体	水体	水体	水体	水体
1	常绿针叶林	常绿针叶林	草地/谷类作物	常绿针叶林	常绿针叶林
						2	常绿阔叶林	常绿阔叶林	灌木	常绿阔叶林	常绿阔叶林
3	落叶针叶林	落叶针叶林	阔叶作物	落叶针叶林	落叶针叶林
						4	落叶阔叶林	落叶阔叶林	草原	落叶阔叶林	落叶阔叶林
5	混交林	混交林	常绿阔叶林	一年生阔叶林	灌丛
						6	郁闭灌丛	郁闭灌丛	落叶阔叶林	一年生草地	草地
7	稀疏灌丛	稀疏灌丛	常绿针叶林	荒漠	谷类作物
						8	多树草原	多树草原	落叶针叶林	城市用地	阔叶作物
9	稀树草原	稀树草原	荒漠		城镇建设用地
						10	草地	草地	城市用地		冰雪
11	永久性湿地				裸地/低植被覆盖地
						12	农田	农田
13	城镇建设用地	城镇建设用地
						14	农田/自然植被镶嵌
15	冰雪
						16	裸地/低植被覆盖地	裸地/低植被覆盖地

注：标准1:IGBP,国际地圈生物圈计划-全球植被分类方案；标准2:UMD,马里兰大学方案；标准3:LAI/fPAR,基于叶面积指数/光合有效辐射吸收比例的MODIS方案；标准4:NPP,基于净初级生产量的MODIS方案；标准5:PFT,植物功能类型方案。

表2对LUCC和MODIS进行重分类规则表

新地类	MODIS数据对应地类（标准1）
		1农	12农田14农田/自然植被镶嵌
2林	1常绿针2常绿阔3落叶针4落叶阔5混交林6郁闭灌丛7稀疏灌丛8多树草原9稀树草原
		3草	10草地
4湿	0水体11永久性湿地15冰雪
		5城	13城镇建设用地
6荒	16裸地/低植被覆盖地

3.规则网格

网格尺度即分析窗口大小的选择非常重要，它直接影响了分析窗口的大小和信息提取的精度。最佳的网格尺度是既能避免对象多边形太破碎，又能保证A-LUCC和M-LUCC在多个时间段均有相同的土地利用转换类型(图1)。本发明通过每次设定不同的网格宽度和高度值，如(1×1),(4×4),…等来把整个研究区划分成i×i,i+1×i+1,…等规则亚区，如果第n+1个亚区内两套数据在多个（至少两个）时间段有相同的土地利用转类而第n个亚区内没有，那么n+1×n+1就是所需的最小分析窗口。

4.时间变化速率的判定

时间变化速率的判定在研究区不同的分析窗口内逐一进行。对不同土地利用类型的时间变化速率，以A-LUCC的时间段变化值为限定值，参考M-LUCC年际变化值，采用对应数据差值的方法来计算，这里称之为逐年类归一化线性拟合处理，具体推导过程如下：

(1)

(2)

(3)

(4)

5.变化栅格的空间定位

其中，质心点坐标（X _cm ,Y _cm）的计算公式为：

(5)

式中，m为多边形第i个顶点的质量；x _cm,y _cm分别为对应多边形顶点i的坐标。

由于在平面二维几何中，平面几何多边形图形可以认为是质量均匀的物体，因此多边形的质心和重心保持在同一点上，且质量为面积和密度之积：，因此，上式可以表示如下：

(6)

式中，S为多边形第i个顶点的面积。基于此，求算质心问题转化为求算多边形的面积。

任意三角形的三个顶点坐标A(x ₁ ,y ₁ ),B(x ₂ ,y ₂ ),C(x ₃ ,y ₃ )，按照逆时针方向排列，则三角形ABC面积为：

(7)

对应三角形的质心表达为：

(8)

任意四角形的四个顶点坐标A(x ₁ ,y ₁ ),B(x ₂ ,y ₂ ),C(x ₃ ,y ₃ ),D(x ₄ ,y ₄ )，则其面积可以通过求算两个三角形面积之和来获取，最终表示如下：

将上式带入公式(6)，则对应四角形的质心为：

(9)

任意边n边形顶点A _i (x _i ,y _i ),(i=1,2,…,n)依照逆时针排列，则其面积可表达为：

(10)

对应任意边n边形的质点可表达为：

(11)

以A-LUCC空间分布为准，同分析窗口下M-LUCC逐年的同转类斑块分布位置作为参考，在变化速率值的限定下，利用条件语句判定空间转换位置：如果该斑块的质心和A-LUCC连续两期研究时段内的同类斑块的质心在规则网格中心点的同侧，那么该转类斑块保持不动；如果在规则网格中心点的两侧，则将该转类斑块移动到同窗口下A-LUCC的最邻近同转类斑块，超出部分仍依照原M-LUCC的空间分布进行同分析窗口下的随机分布定位。

6.结果的精度检验和时间序列分析

利用研究时段内的A-LUCC数据，随机选取研究区中某一年的数据，对应检验各类型的空间分布和面积是否一致，进而给出精度。或者更进一步，随机选取研究区中的某一年的人机交互解译A-LUCC数据，对比检验同一年份的新制备数据，或者在研究区内随机选几个小区域找到逐年的较高分辨率影像进行人机交互解译，再对比解译的分类和本方法制备的分类：主要分析各类型的空间分布和面积是否一致，给出对应精度。

土地利用长时间序列时空分析是指根据制备出来的数据结果，看其空间变化情况，并统计该研究时间段的各土地利用类型的数量和占总面积的百分比，从而进一步得出土地利用类型转换概率矩阵和转换面积矩阵。

本发明包括权利要求中所引用的主旨的所有修改及等效形式，这在适用的法律中是允许的。此外，上述要素的所有相关领域技术人员可能进行的步骤变更，及任何组合也被本发明所包含，除非另外指出或在上下文中明显差异。

Claims

1.一种基于多源数据的逐年土地利用空间数据的制备方法，其特征在于，包括以下步骤：

1）数据获取：

两种不同源数据：一是利用较高分辨率对地观测卫星影像30m空间分辨率的LandsatTM/ETM+，通过人工解译制备的土地利用变化数据集（A-LUCC）；二是500m空间分辨率MODIS的土地覆盖数据（M-LUCC）；

2）数据预处理：首先对多期A-LUCC数据进行地理空间分析算法中的矢量转栅格，并基于最近邻算法进行栅格数据重采样，即将空间分辨率重采样为500m，以保证和M-LUCC的空间分辨率一致；

然后，对A-LUCC和M-LUCC的土地利用类型，采用大类的分类标准，分为农田、林地、草地、湿地（含水域）、建设用地和其他未利用地等六大类（简称农、林、草、湿、城、荒）；

3）预设分析窗口：本发明通过均匀网格化来设定分析窗口，即把整个研究区分成1×1,4×4,…等亚区作为分析窗口，而判断最小或最优精度的分析窗口大小的依据为：在该窗口内同时具有A-LUCC和M-LUCC的土地利用类型；

4）时间变化的判定：

对不同土地利用类型的时间变化速率TR _ij，以LUCC的时间段变化值为限定值，参考M-LUCC年际变化值，采用对应数据差值的逐年类归一化线性拟合处理：

式中，MTR为M-LUCC对应转类i的对应年份j+n的变化速率；(MTR _j+n -MTR _j+n-1)、(MTR _j+1 -MTR _j)、…为逐年际间各土地利用转类的面积变化差值；则，如果MTR _{ni_sum}≠ATR _{ni_sum}时，需要采用如下方法处理：

式中，j+m，j+n和j分别为研究时间段的任意年份、末年和首年对应的土地利用类型i变化值，其中1<m<n；则有此式可以推导出逐年LUCC类归一化的土地利用类型i对应的面积值为：

式中，TR _j+m为LUCC数据中某土地利用转类i的面积；ATR和MTR分别为A-LUCC和M-LUCC中某一转类i的面积；

5）空间变化的判定：以A-LUCC空间分布为准，同分析窗口下M-LUCC逐年的同转类栅格分布位置作为参考，在变化速率值的限定下判定空间转换位置；即：

如果该斑块的动态分布范围与在A-LUCC连续两期研究时段内(t和t+1)的同类斑块(Grid)动态分布范围一致，则为正确的动态转类空间定位，即任意栅格位置对应满足：Grid _t -Grid _t+1 =0；否则，空间范围不一致时，则定位该转类斑块到同窗口下A-LUCC的最邻近同转类斑块的所属空间，超出A-LUCC对应转类斑块空间的依照原M-LUCC的分布空间，进行同分析窗口下的质点所在空间趋势位置分为八个方向(东、南、西、北、东北、西北、东南、西南)，再随机分布定位：即，当质点偏向某一方向，就在某方向所在位置随机布置空间转类栅格类型；

其中，边形的质点求算可以通过以下公式表示：

设，任意边n边形顶点A _i (x _i ,y _i ),(i=1,2,…,n)依照逆时针排列，则，其对应任意边n边形的质点可表达为：

6）精度验证：利用研究时段内的A-LUCC数据，随机选取研究区中某一年的数据，对应检验各类型的空间分布和面积是否一致；或随机选取某年或某区域内人机交互解译的A-LUCC数据，以此对比检验同一年份的新制备数据，给出对应精度。