基于时序过程的沿黄粮食主产区土地利用分异研究

杨崇科 马泉来 周 浩 万小强 曹艳杰

(1.河南省资源环境调查一院, 郑州 450007;2.河南省自然资源科技创新中心(资源环境承载能力评价与监测预警研究), 郑州 450007;3.湖南师范大学地理科学学院, 长沙 410081)

0 引言

探查特定地区土地利用变化规律,能够发现土地利用行为异常、实现土地资源优化配置[1]。自21世纪初以来,我国城镇化和工业化进程持续推进,土地利用变化表现为传统用地类型减少、景观同质化显著等特点[2-3]。其中,对于主要粮食产区而言,由于粮食生产与经济建设的客观矛盾,土地利用变化不仅表现更为突出的建设占用耕地、林草地退化等空间维度差异,而且在全域时序过程上也呈现高度复杂性[4-6]。因此,有必要围绕该类型区开展土地利用变化的内在规律挖掘和多维表达研究,以促进土地资源优化配置和保障粮食安全。

作为全球变化研究的组成部分和重要动因,面向粮食生产需求的土地利用变化研究已引起了学者高度关注。当前,相关研究主要集中在经济发达区或生态脆弱区,多围绕土地利用变化的“形”(空间形态)、“量”(类型与数量)和“理”(演化机理)等,采用如结构比例、转移矩阵、土地利用动态度、热点分析、地统计学分析等模型来简化土地利用过程,分析其基本模式、数量结构和空间格局,为深化土地利用变化研究提供了基础。如宋戈等[7]引入空间形态学中的标准差椭圆,来研究三江平原耕地利用变化的方向性特征;周浩等[8]通过动态度指数来研究挠力河流域相邻年份间土地利用变化剧烈程度;宋小青等[9]采用转移矩阵来描述耕地在不同时间节点上的数量转型特征;丁智强等[10]运用信息图谱论揭示了云南哈尼梯田长时间序列的土地利用时空变化特征及其地形梯度效应。但上述研究多侧重于研究对象整体或内部结构单元的状态特征描述(如结构比例、动态度指数等指标仅描述了研究对象的整体变化速率,无法细致度量单位内部的异质性),在时空过程解析上重格局轻过程、重模拟轻度量,定量分析仍显不足。同时有研究发现,通过定量表达多时空条件下土地利用的“空间格局”和“时序特征”,对土地利用变化过程进行量化,有利于进一步深化土地利用变化研究[11-12]。因此,从时序过程视角出发,对土地利用变化的时空过程进行度量,能够探查土地利用状态分析中难以发现的差异,从而揭示时空过程潜在规律,以深化土地利用变化研究。

黄河中游沿岸作为我国传统农区和新时期高质量发展战略区,在推动国家生态文明建设及商品粮供应安全上具有重要作用[13-15]。长期以来,经济发展建设及持续性生态退耕政策等综合作用,导致该地区土地利用变化显著且时序过程趋于复杂。基于此,本文以位于黄河中游沿岸的河南6县(以下简称“沿黄6县”)为研究区,从时序过程视角构建累积动态度模型以拓展土地利用变化研究思路,并采用时空栅格建模、空间点模式分析等方法,定量揭示1980—2020年该地区土地利用变化过程特征,挖掘潜在分布模式,以加强对当地土地利用配置的认知,为促进土地资源优化配置提供理论支撑。

1 研究区概况与数据来源

1.1 研究区概况

沿黄6县位于黄河中部的干流北岸、河南省境内,辖原阳县、封丘县、延津县、长垣县、滑县和濮阳县6个县级行政区(图1),总面积7 657.15 km2。地势以平原为主,属温带大陆性季风气候区,四季分明,但降水年际和区域分布不均,易形成季节性干旱。该地的农业开发利用条件较为优越,农业生产空间占国土总面积的近80%,并已成为河南省乃至我国最重要的粮食主产区之一。根据第七次人口普查数据,2020年其常住人口已达到473.37万人,对应GDP为1 739.94亿元,土地利用及开发程度较高,受人类活动干扰强烈。因此,在相对较小的时空尺度,该地区富集了本文所需的自然、经济、政策等信息,具有显著的地域特色和代表性。

图1 研究区地理位置示意图Fig.1 Location of research areas

1.2 数据来源

土地利用原始信息源均来自美国陆地资源卫星Landsat多光谱遥感影像,数据获取自美国地质勘探局(USGS)(http:∥earthexplorer. usgs.gov/),过程中选用不同的卫星数字产品,在影像的几何纠正、图像增强等预处理基础上,采用不同波段的标准假彩色融合以实现人机交互式目视解译,并结合后期的室内数据纠正,完成1980、1990、2000、2010、2020年共5期土地利用信息识别,具体划分类型包括耕地(包括旱地和水田)、林地、草地、水域、建设用地(城镇建设用地、农村居民点和其他建设用地)和未利用地6大类;解译精度验证方面,由于早期(1980年和1990年)图件数据等资料较少,主要通过比对中国科学院资源环境科学数据中心(http:∥www. resdc.cn)“中国土地利用现状遥感监测数据”(30 m)以及历史文本资料进行,中期(2000年和2010年)则通过Google Earth软件布控数据的采样网格验证点来实现(解译准确率均大于85%),而对于现状年(2020年)的验证,通过对实地GPS信息样点的比对验证和记录(考察时间为2020年8月7—13日),完成土地利用解译数据的精度验证(解译准确率89.20%);基础地理信息数据包括行政区划数据(县级和乡镇级)、黄河水系分布以及基础道路交通数据,均来自于当地自然资源管理部门。研究中各图件经投影变换统一转为Albers双标准纬线等积投影。

2 研究方法

通过栅格化处理相邻年期的土地利用变化过程,以揭示土地动态演化在空间维度上的异质性规律,继而通过累加研究期不同时段的空间过程,实现土地利用变化的时序过程探查。

2.1 土地利用动态度

2.1.1动态度计算

动态度客观描述了土地利用变化状态的总体特征,其中单一动态度D为某地类在研究期的面积变化速率,综合动态度C则为所有地类在研究期的整体面积变化快慢情况[8],计算公式为

(1)

(2)

式中Ua——对应地类在研究初期面积,km2

Ub——该地类在研究末期面积,km2

T——研究时段长度,a

CUi——研究初期i种地类面积,km2

ΔCUi-j——i种地类变为j种地类面积,km2

2.1.2动态度栅格化处理

由于动态度无法细致度量土地利用变化的内部空间过程,因此借鉴景观生态学中能反映景观格局特征的网格采样法,并考虑平均斑块面积(平均斑块面积的2～5倍能较好地反映采样区周围景观格局的异质性特征[16]),确定网格采样尺寸为2 km×2 km(研究区边缘的不规则地块,面积大于0.5个规则网格的单独作为一个样区,面积小于0.5个规则网格的并入相邻的样区,最终共有样区2 100个),来网格采样C;继而将对应网格的C视作该网格中心点值,并运用克里金插值法(Kriging),实现沿黄6县综合动态度的栅格化处理(分辨率为30 m)。

2.2 CDDLU模型构建

为实现土地利用变化的三维时空过程向二维平面压缩,对所有相邻时段动态度差值的绝对值进行求和处理,即可获得土地利用累积动态度(Cumulative dynamic degree of land use,CDDLU),以深化土地利用变化研究。CDDLU反映不同时序过程所对应的土地利用变化快慢程度,CDDLU越大则土地利用变化时序过程越强烈,采用K-means聚类法对CDDLU进行等级类型评价,以保证每一类型的差异最小,进而挖掘土地利用变化的时序过程模式。公式为

CDDLU=∑|Ct-Ct+1|

(3)

式中t——年份时点由小到大的序列数,本文1980、1990、2000、2010、2020年的t值依次为1、2、3、4、5

Ct——对应年份时点栅格化处理后的综合动态度

2.3 时序过程模式探索方法

研究区属黄河干流北岸豫北平原的核心粮食主产区,农业生产空间规模大,其境内黄河以及国省道、高速公路、铁路等基础道路地理要素对该地区土地利用格局形成影响显著,因此本文选取黄河以及主干道路等线性地理要素,来分析土地利用变化时序过程模式的近邻缓冲尺度效应。同时,在具体模式效应探索过程中,空间点模式分析中的最近邻指数(NNI)[12]通过计算最邻近点对应的平均距离,来比较该点的观测模式与随机模式之间的相似性,进而判定点的分布模式状态,在空间尺度分析上具有较好的探索性和应用性。因此,本文选取该指数来研究单一空间尺度下不同土地利用变化类型的临近缓冲特征,公式为

(4)

式中di——第i个网格与其最近邻网格之间距离,取1、2、3、5、10 km

n——动态度网格采样数目

A——研究区面积

NNI=1,类型点处于随机分布状态;NNI<1,类型点处于集聚分布状态;NNI>1,对应类型点处于均匀分布状态。为了更好地反映实测平均距离和预期平均距离偏离程度,用正态分布检验得出Z值及其置信水平,Z值为负值且越小,则要素分布越趋于集聚分布;Z值为正值且越大,则要素显著偏向均匀分布;位于二者间则随机分布。

由于点分布模式的空间尺度性效应,使得点在不同缓冲距离上,其所表现的聚类效果存在差异,因此,采用多距离空间聚类的Ripley’s k指数法[17],通过比较某一空间尺度的观测K值和期望K值的相对大小确定类型点的分布模式,来描述多空间尺度所对应的土地利用变化时序过程特征,以最大限度挖掘空间点信息,公式为

(5)

其中

(6)

式中r——空间尺度

uij——点i和点j之间的距离

为了更好地解释实际空间格局,提出了L函数代替K函数,由此可得

(7)

式中d——距离

Ki,j——权重,如果没有边校正,当i和j之间的距离小于d时,Ki,j=1,反之Ki,j=0

在L(d)变换下,预期值K等于距离。当处于置信区间内,如果特定空间尺度的K观测值大于K预期值,则该类型点分布的聚类程度更高;如果K观测值小于K预期值,类型点分布的离散程度更高。

3 结果与分析

3.1 1980年以来沿黄6县土地利用现状

利用ArcGIS 10.5对1980年以来土地利用现状数据进行汇总统计。由表1和图2可知,40年来,沿黄6县土地利用结构以旱地、水田、建设用地和水域用地为主,林地、草地以及未利用地面积占比较低,土地利用发生了明显变化。其中,建设用地面积持续增加,且尤以城镇建设用地最为迅速,空间上主要集中于长垣县、延津县、原阳县等建成区周边地带。水田、旱地和林地面积均呈先增加后减少特点,草地、水域用地和未利用地面积持续减少,主要是由于中期黄河沿岸耕地后备资源的持续性开发至殆尽以及经济快速发展导致。具体表现为:

表1 1980—2020年各土地利用类型面积和动态度Tab.1 Statistics on area and dynamics of various land use types from 1980 to 2020

图2 1980—2020年研究区土地利用现状Fig.2 Land use status of study area from 1980 to 2020

(1)长期以来,作为粮食主产区的沿黄6县,耕地保护目标战略执行较好,整体数量格局稳定,其二级类的旱地、水田面积占比分别由1980年的73.29%、6.50%变为2020年的73.10%、6.30%;同时,沿黄6县属河南省农业人口相对稠密区,境内城市化、工业化进程持续推进,以及郑州市、上海市等发达城市对该地区的人口虹吸作用,使得其城乡发展呈显著的“二元化”特点,并表现在建设用地内部结构变化差异上。40年间,该地区建设用地中城镇建设用地扩张显著,其面积占比由1980年的0.42%上升至2020年的2.29%,与之对应的是农村居民点面积占比仅增加0.15个百分点,二者相对面积扩张比达到5.33∶1;对于生态性用地而言,林地和水域用地面积波动变化,其中林地面积先增后减,而水域用地自2000年来受当地退耕还湿等政策影响,研究期内面积先减后增,草地面积占比由1980年的0.17%降至2020年的0.07%,未利用地则由1980年的0.82%变为0.02%。

(2)从动态度来看,各时段的土地利用变化状态存在较大差异。1980—2020年,沿黄6县水田和旱地面积缓慢下降,动态度分别为-0.08%和-0.01%,面积下降时段主要集中在2000年后;建设用地中城镇建设用地的动态度为所有地类最高(10.97%),其他建设用地为4.90%,二者均表现出快速扩张的特点,而农村居民点面积呈增加—减少—增加态势,平均动态度仅为0.03%;草地面积持续下降且尤以2000—2010年最为突出,年均变化率达到-4.72%,水域用地先减后增,未利用地先增后减,二者平均动态度分别为-0.64%和-2.44%。

在空间上,草地、水域用地以及未利用地变化区均集中分布于原阳县和封丘县的黄河沿岸带,反映了黄河对该地区生态性用地格局形成的重要作用。

3.2 土地利用累积动态度

3.2.1动态度空间格局

计算各采样网格中心点的土地利用动态度并运用Kriging空间插值法,实现动态度的栅格化处理。从图3可看出,经过栅格化处理的动态度,能够较好反映沿黄6县在不同时段的土地利用变化特点。其中,1980—1990年和1990—2000年2个时段的土地利用变化较为平缓,整体动态度偏低。但由于原阳县和封丘县南部以及长垣市东南等地邻近黄河,受黄河沿岸的水域湿地变化影响较大,使得该地区动态度的高值区多集中于此,且1990—2000年的变化强度高于1980—1990年;自21世纪初以来,沿黄6县土地利用变化加剧明显。2000—2010年,南部黄河沿岸以及各主要乡镇周边土地利用开发活动强烈,动态度偏高,其相对高值区(大于3%)主要分布在封丘县的荆隆宫乡、曹岗乡、李庄镇以及长垣县的恼里镇。2010—2020年,土地利用变化强度有所放缓,但高值区更集中分布于各县域的建成区周边地带。

图3 1980—2020年研究区不同时段土地利用动态度分布Fig.3 Dynamic degree of land use in different periods in study area from 1980 to 2020

3.2.2累积动态度时序过程

基于各时段栅格化处理的动态度数据,实现沿黄6县土地利用累积动态度(CDDLU)提取,采用K-means聚类法并依据CDDLU,划分3个土地利用变化过程类型(类型Ⅰ、Ⅱ、Ⅲ)(图4)。研究结果表明,CDDLU能较好表征长时间序列土地利用变化的时空过程信息,1980年来,沿黄6县土地利用变化呈显著的空间异质性和阶段性分布特点,尤其对于黄河沿岸及若干主要县域建成区周边地区而言,土地利用变化频繁且强烈,差异性时序过程特征突出。具体表现为:

(1)CDDLU范围为0.16%～7.30%,其中高值区主要分布在黄河沿岸及其邻近缓冲区,受人类农业开发活动以及黄河湿地退化等因素影响较大,同时各县域建成区周边的土地利用开发活动频繁,土地利用稳定性差,导致对应的CDDLU也偏高。中部、北部等地的CDDLU偏低,主要是由于该地区为传统商品粮种植区,土地利用变化时序过程平缓。从行政区来看,封丘县的CDDLU均值最高,达到1.38%,位于其境内的李庄镇、恼里镇和曹岗乡以及位于延津县的新乡食品纺织工业聚集区的CDDLU均值依次为3.77%、2.87%、2.66%和2.96%,土地利用变化过程强烈。然而,作为河南省县域粮食产量最高的滑县,CDDLU仅为0.64%,土地利用变化时序过程平缓,实地考察资料显示,该地区支柱性产业以农业生产及其次级加工为主,耕地数量保护较好。

(2)类型Ⅰ(CDDLU为2.97%～7.30%)所对应的土地处于剧烈时序过程变化状态,面积占比为5.44%,主要分布在南部黄河沿岸,西北部也分布着少量该过程类型土地。上文分析显示,该地区土地利用变化以生态性用地(草地、水域用地和未利用地)为主,强烈的生态性用地变化过程,易导致其生态系统脆弱。因此,类型Ⅰ所属地区应尤为关注“黄河流域生态保护和高质量发展”国家级战略契机,以促进生态性用地趋于稳定及良性结构发展;类型Ⅱ(CDDLU为1.04%～2.96%)面积占比为23.75%,时序过程较为明显,土地利用变化频繁,空间上呈环研究区边缘地带分布的特点,同时原阳县、封丘县、长垣县等县域中心建成区周边地带也分布一定面积的该类型土地,后续国土空间管理与生态保护建设需重点关注该类型区对应的土地;类型Ⅲ(CDDLU为0.16%～1.03%)为该地区最主要的土地利用变化时序过程类型,面积占比达到70.81%,广泛分布于研究区中部、南部等地,与耕地分布存在较高重合性,时序过程较为平缓,需重点维持其对应的用地功能和结构面积稳定性。从行政区来看,长垣县、封丘县和原阳县的土地利用变化显著区(类型Ⅰ和类型Ⅱ)依次占对应土地总面积的48.61%、48.34%和46.17%,而滑县的变化显著区面积占比仅为9.52%,各县域单元时序过程的类型差异明显。

3.3 土地利用变化聚集性模式

将CDDLU赋值给网格中心点,继而采用NNI指数,来探索单一空间尺度所对应的土地利用变化过程类型的聚集性特点(表2和图5)。从时序过程的全域NNI指数来看,沿黄6县土地利用动态演化特征突出。其中,类型Ⅰ和类型Ⅱ的时序过程显著,对应样点数分别为160个和530个,在99%置信度下,二者均处于聚集分布状态。类型Ⅲ的时序变化过程平缓,对应样点数为1 384个,分布均匀,同样呈现统计显著性;从河网、公路网基本线性要素的不同缓冲区NNI统计结果来看,在距离线性要素较近的范围内(1 km和2 km),类型Ⅰ、Ⅱ、Ⅲ均处于聚集分布状态,但随着缓冲区范围增加,类型Ⅰ的聚集性程度先增后减,类型Ⅱ的聚集性则持续减弱,类型Ⅲ的聚集性同样持续减弱,并同时导致其分布模式由聚集分布逐渐转变为均匀分布,说明土地利用变化的时序过程受距离道路、河流远近影响显著,但影响类型和范围均存在差异。

表2 3类土地利用变化时序过程的NNI指数Tab.2 NNI of spatial and temporal processes for land use change types

图5 研究区线性要素在不同缓冲区范围内的NNI指数Fig.5 NNI statistical results of linear elements in different buffer ranges in study area

采用Ripley’s k指数描述多空间尺度(1、2、3、5、10 km)对应的时序过程类型的聚集性特征,以研究沿黄6县土地利用变化空间模式的尺度性效应。从图6可以看出,在分析尺度1 km下,3个时序变化过程类型的观测K值均小于预期K值,表示该空间尺度下各时序过程类型均处于离散状态,但随着空间尺度递增,其差值转为正数且绝对值也越来越大,聚集性逐渐提高。同时,在不同空间尺度下,类型Ⅰ、Ⅱ、Ⅲ的离散和聚集性程度也存在差异。类型Ⅰ聚集性的递增速度最快,且至分析尺度9 km,其观测K值与预期K值的差值达到最大(10 073.30),而对应的类型Ⅱ为4 543.53。类型Ⅲ在分析尺度1 km的离散程度最高,但随着空间尺度增大,聚集性递增程度在3个类型中最低,其在9 km处的观测K值与预期K值差值仅为2 879.82。上述分析结果表明,沿黄6县土地利用变化程度受空间尺度影响较大,但变化越强烈的地区,其受空间尺度影响程度越弱。因此,沿黄6县在后续土地利用开发管理中,应在重点考虑线状要素(黄河和主干道路)布局的前提下,根据线状要素近邻程度特点,制定差别化的利用管理措施,以促进土地利用优化布局。

图6 研究区类型Ⅰ、Ⅱ、Ⅲ在各缓冲区下的Ripley’s k指数统计结果Fig.6 Statistical results of Ripley’s k indices for land use change types (Ⅰ, Ⅱ and Ⅲ) under various buffers in study area

4 讨论

土地利用变化为自然、社会和经济等因素驱动下的复杂性过程。本文从动态度着手,通过对单一状态过程的差值累加来实现土地利用变化的时序过程提取,为土地利用变化的时空模式挖掘提供了新思路。但由于地理事件或地理现象的空间统计模式存在尺度依赖,本文得出的结论只是在特定采样尺度下(2 km×2 km)的特定结论,是该尺度下土地利用变化的时序过程规律的表征,但采样方式和网格单元的取值方法不同,或会对土地利用变化的时序过程模式结论产生影响,采样方式取决于其对时空模式结论的影响最小,因此,这种影响的规律性需要进一步定量研究[18-21];其次,CDDLU模型是本文提出的概化土地利用变化时序过程核心手段,但由于不同类型划分结果会对土地利用变化的时空模式的度量产生影响,如何依据CDDLU来科学划分不同过程类型,为影响时空模式挖掘的关键,后续需全面分析不同分级结果对土地利用变化模式度量的影响;再次,河南沿黄6县位于黄河以北,是我国重要的粮食主产区,其土地利用格局形成受黄河、基础交通道路信息影响较大,为保证时空模式挖掘的简洁性,本文选取包括黄河和道路因子在内的线性要素进行分析,重点讨论了上述线性要素的不同邻近性对土地利用变化时序过程的影响,但这种随着线性要素邻近性的时空模式还受到其他因素影响,后续需进一步对其他要素进行论证以提高模式的精确性。

5 结论

(1)河南沿黄6县土地利用变化呈显著的空间异质性与阶段性分布特点。耕地面积缓慢下降,城镇建设用地扩张明显,水域用地和林地面积波动变化,同时在前期黄河沿岸的农业土地开发活动尤为强烈,至后期,原阳县、封丘县和长垣县等县域建成区周边地带土地利用变化强烈。

(2)CDDLU能较好表征长时间序列土地利用变化的时空过程信息。类型Ⅰ、Ⅱ均处于聚集分布状态,面积占比分别为5.44%、23.75%,类型Ⅲ均匀分布;变化时序过程受道路、河流距离影响较大,但影响类型和范围存在差异。距离线性要素较近范围的类型Ⅰ、Ⅱ、Ⅲ均处于聚集分布状态,但随着缓冲区范围增加,类型Ⅰ的聚集性程度先增后减,类型Ⅱ、Ⅲ则持续减弱;变化时序过程受空间尺度影响较大,变化越强烈地区,受空间尺度影响的程度越弱。分析尺度1 km下,各类型均处于离散状态,但随着空间尺度递增,聚集性逐渐提高;沿黄6县土地利用开发管理,应适度考虑黄河和主干道路等线性要素近邻特点,制定差别化的利用管理措施,以促进当地土地利用优化布局。