基于时序MODIS的黄河上游2002—2018年耕地时空变化特征分析*

韩春雷,沈彦俊,武兰珍,郭 英,陈晓璐,4

(1.河北师范大学资源与环境科学学院 石家庄 050024;2.中国科学院遗传与发育生物学研究所农业资源研究中心 石家庄 050022;3.甘肃农业大学水利水电工程学院 兰州 730070;4.青海师范大学地理科学学院 西宁 810008)

耕地是社会生产发展的最基础资源,其不仅仅是粮食的来源、人民生活的基础物质保障[1],还具有空间承载、生态服务等多种价值。耕地的分布及其变化对区域水资源的利用和生态环境的变化具有重要影响,不合理的耕地利用会导致水土流失、耕地质量下降、水资源供需紧张、地下水位下降等一系列问题[2-3]。准确监测耕地分布的变化是区域水土资源优化管理的重要前提,对于构建和保持区域或流域合理的山水林田湖草国土利用空间格局,实现区域或流域生态与经济协调发展具有重要的意义。

黄河流域是我国重要的生态屏障区,水资源严重匮乏,缺水和水资源时空分布不均给全流域的高质量发展带来挑战。黄河上游产出全流域56.77%的水资源[4],而农业用水占全流域40%以上[5],农业产出不高,水资源利用效率低下等问题突出。黄河上游处于干旱半干旱地区,降水量少,不能满足耕作需求,大量灌溉用水也加剧了水土流失问题。同时,农业开发与生态环境保护和水资源高效利用矛盾突出,因灌溉导致的次生盐碱化已成为影响耕地生产力可持续维持的突出问题。因此,研究过去几十年黄河上游农业用地的时空变化对于合理管理水土资源,制定科学合理的水资源利用规划具有重要意义。

遥感和地理信息技术在耕地变化的监测上已得到广泛的应用,对于提升区域水土资源管理起到重要的支撑作用。国外学者Fritz 等[6]通过收集国家与全球尺度的耕地数据,制作了以2005年为基准的1 km 全球IIASA-IFPRI 耕地百分比地图,且用不同方式验证有较好的精度;Clark 等[7]提出了一种土地利用和土地覆盖制图方法,通过结合MOD13Q1 数据,用该方法在阿根廷、玻利维亚和巴拉圭的干查科生态区进行试验,结果的精确度较好并指出大豆(Glycine max)和种植牧场扩张是造成林地毁坏的因素;Arowolo 等[8]利用遥感数据、气象数据、统计数据等探讨了尼日利亚耕地变化的驱动因素;Belgiu 等[9]基于高分辨率Sentinel-2 数据,使用TWDTW(timeweighted dynamic time warping)方法在罗马尼亚、意大利和美国分别进行基于像素和基于对象的各种作物类型分类,并与随机森林算法结果进行比对,最终得到基于对象的效果要更好。国内也有苏锐清等[10]利用Landsat 影像分析了潮白河区域2001−2017年耕地利用时空变化;王凤娇等[11]利用遥感与统计数据探讨了黄土高原耕地变化与粮食安全之间的关系;赵明明[12]利用无人机影像和高分一号影像对黑龙江海伦市的农田耕地进行提取,对比分析发现无人机影像更适于进行农田精准管理区划分。在黄河流域,也有耕地变化轨迹[13]、土地利用[14-16]、生态功能区耕地治理[17]等相关研究,而这些研究中大多使用5年一期的土地利用数据或几个时间节点的数据,虽然空间分辨率精度相对较好,但还是缺少长期的动态监测,在时间的连续性上有所欠缺。因此进行连续耕地变化监测能够更好地反映耕地长时间的时空差异,提供更为细致全面的耕地利用决策支持。

本文采用决策树分类方法,以黄河上游龙羊峡以下地区为研究区域,借助遥感技术提取了黄河上游2002−2018年整体的耕地分布,分析了耕地的时空变化特征,为黄河上游生态环境保护、高质量发展、土地合理规划、水资源合理规划等方面提供依据。

1 材料与方法

1.1 研究区概况

黄河干流全长5464 km,流域总面积79.5 万km2。黄河上游是指从源头到内蒙古头道拐之间的流域范围,其中唐乃亥水文站以上部分称为源流区,主要为青藏高原的草甸草原,是主要的牧业区。本文研究区集中在龙羊峡以下的黄河上游地区主要农业耕地分布区,该区位于100.68°～111.52°E,34.12°～41.84°N(图1)。研究区地跨青藏高原、黄土高原两大地形区,地势整体西南高、东北低,起伏较大,甘青地区山脉河谷较多,宁蒙地区多以平原为主,北部地区被巴丹吉林沙漠和乌兰布和沙漠环绕。气候随地形分为高原气候与大陆性气候,降水多集中在每年的6−9月,占年降水的75%以上,多年平均降水量从青海的446 mm,递减到河套平原的170 mm。

由于研究区地形起伏明显、南北跨度较大,使得研究区内温度、光照、降水等条件有一定差异,东西、南北各地区耕地上会出现同一作物有不同的生育期或不同作物有相近的生育期情况,因此将研究区再分割出小区域以减少作物生育期差异对分类的影响。小区域分别为:内蒙段、宁夏段、甘肃段和青海段。

1.2 数据来源

研究所用数据包括遥感数据、作物生育期数据、社会经济数据、实地考察数据。其中遥感数据为MOD13Q1 数据集中NDVI 数据,空间分辨率为250 m,时间分辨率为16 d,研究使用2002−2018年的各年数据,一年共有23 景,第一景为“B1”、第二景为“B2”,依次类推。作物生育期数据包括作物类型、播种日期、收获日期等。社会经济数据包括研

究区内各省、县级的作物面积等。实地考察数据主要为实地获取的作物类型点。表1 展示了各个数据的来源。

表1 研究数据来源Table 1 Research data sources

1.3 研究方法

1.3.1 时间序列谐波分析法(Harmonic Analysis of Time Series,HANTS)

MOD13Q1 数据集使用最大值合成法得到16 d一期的数据,经过该方法处理已经消除了一部分噪声,但是仍然存在一定不足。为进一步提升数据质量,减少数据噪声影响,使用时间序列谐波分析法对MOD13Q1 数据进行平滑去噪。该方法的核心算法是傅里叶变换和最小二乘法拟合,即把时间波谱数据分解成许多不同频率的正弦曲线和余弦曲线,从中选取若干个能够反映时间序列特征的曲线进行叠加,以达到时间序列数据的重建目的[18]。通过调整参数,可以观察曲线变化,且能够在时间序列中任何时刻再现无云图像[19]。

式中:谐波的余项A0等于序列的平均值;Aj为各谐波的振幅;ω=2jπ/N为各谐波的频率;N为序列的长度;θj为各谐波的初相位;m为谐波个数,m=N−1。这些正弦函数叠加构成傅里叶序列[20]。

HANTS 的工作过程是:首先由所有离散数据生成最小平方拟合曲线,然后检查每一个数据值,将它与曲线进行比较;图像上受云影响的点NDVI 值会很低,明显低于拟合曲线的点是受云干扰的点,要剔除并赋0 值;偏离量超过阈值最大的点优先剔除,然后根据剩余的采样点重新生成拟合曲线,再检查每个数据值,再剔除偏离曲线超过阈值的点,反复循环,最后生成光滑曲线[21]。

MOD13Q1 数据集的时间分辨率为16 d,一年则共有23 个波段,将每一个波段所对应的NDVI 值连接即可生成NDVI 时间序列曲线。但是原始的NDVI 时序数据不够平滑,时间序列曲线较为曲折,较难判断作物生育期,利用HANTS 对其进行滤波平滑,得到较为平滑的NDVI 时序曲线,以耕地上的作物为例平滑前后的NDVI 曲线对比如图2所示。结合作物生育期以及平滑后的NDVI 时间序列曲线可以进行特征识别,为耕地提取提供依据。

1.3.2 决策树分类

使用决策树分类方法进行耕地提取。决策树算法的基本思想是通过一些判断条件对原始数据集逐步二分和细化,其中,每一个分叉点代表一个决策判断条件;每个分叉点下有两个叶节点,分别代表满足条件和不满足条件[22]。通过设置判别条件,计算输入的数据是否满足该节点所设置的条件,满足和不满足分别再进行计算,层层递进,直到没有规则进行判别为止,此时计算结束。该方法快捷、高效、易于操作,已经被众多学者应用。地面典型地物由于自身的光谱特性,在NDVI 时间序列曲线上的反映不同,例如,大部分地物的NDVI 值均大于0,但因水体自身对光线的吸收作用,一般NDVI 值小于0;城市、裸地等地物因植被覆盖相对较低,NDVI 值也相对较低,一般小于0.3[23]。依据地物不同NDVI 曲线进行决策树编写,部分决策树如图3所示。

1.3.3 Mann-Kendall 检验

Mann-Kendall 检验是一种非参数统计检验方法,主要分析要素随时间序列的变化趋势。该方法的优势是样本不需要遵循一定的分布状态,不受少数异常值的影响,计算简便、应用范围广泛[24-25]。本文利用Mann-Kendall 检验对构造的MODIS NDVI 时间序列数据进行趋势检验,分析研究区内植被随时间的变化情况。利用Mann-Kendall 趋势分析方法对研究区内植被变化情况进行分析,设置α=0.05,临界值为±1.96,其中临界值大于0 但小于1.96 为不显著增减,大于1.96 为显著增加;临界值小于0 但大于−1.96为不显著减少,小于−1.96 为显著减少。

2 结果与分析

2.1 耕地提取结果及精度验证

利用决策树分类方法进行耕地提取,结果如图4所示。本文利用外出考察所取得的实地验证点对提取结果进行混淆矩阵验证;为了进一步验证提取结果,将提取结果与研究区内的县级统计数据进行相关分析。利用手机照片所含的GPS 信息在外出考察中取点4190 个,对取得的点进行筛选整理,除去其他用途、重复、无效的点后,可用耕地点有808 个。对提取结果进行混淆矩阵验证,其结果如表2所示,各地区耕地提取结果精度在75%以上。对提取结果和县级统计数据进行相关分析,其结果如图5所示,耕地提取结果与统计数据的R2达0.85,说明结果能够反映研究区耕地时空分布。

表2 耕地提取混淆矩阵验证结果Table 2 Result of cultivated land confusion matrix verification

2.2 研究区耕地时空分布

研究区各年份耕地数量及变化情况如图6所示。2002年到2018年耕地面积呈增加状态,从276.16 万hm2增加到364.37 万hm2,共增加88.21 万hm2,研究区整体增加速率为5.18 万hm2·a−1;宁夏和内蒙地区耕地增长最为迅速,2018年宁夏段耕地总面积达76.61 万hm2,共增加64%,增加速率为1.87万hm2·a−1;内蒙段耕地面积也大幅增加,2018年达145.47 万hm2,占总面积的44%,增加速率为2.46 万hm2·a−1;甘肃段耕地面积同样呈增加趋势,到2018年共增加18.89 万hm2,增加速率为1.05 万hm2·a−1;青海段耕地则为减少趋势,共减少5.36 万hm2,减少速率为0.20 万hm2·a−1。

在空间分布上,青海段地处青藏高原边缘,地形起伏大,耕地主要分布在河谷地区,沿黄河、湟水、大通河等谷地与两侧缓坡分布;甘肃段地处青藏高原和黄土高原的过渡地带,地势复杂多样,耕地主要分布在较为平坦的河谷、山前冲积地区,除沿黄河、洮河、祖厉河、庄浪河等河流分布外,还有景泰县等引黄灌溉地区;宁蒙地区地形较为和缓,耕地主要分布在河套平原地区。宁夏南部和中部因灌溉工程建设,也有耕地分布。

研究区内的耕地变化情况如图7所示,青海段整体上耕地减少,减少主要集中于城市中心地区和坡地,如西宁市区、海东城市沿线、贵德县城等。这些耕地主要是由社会经济发展、城市扩张引起。大通回族土族自治县、互助土族自治县等地区山中坡耕地减少而西宁市区、海东市、同仁县等缓坡地区耕地增加。甘肃段整体上耕地增加,定西市、会宁县地区最为明显,临洮县、榆中县、东乡族自治县也有增加,而天祝县南部山区、永登县山前和城市地区、甘南地区耕地主要呈减少趋势。宁蒙地区耕地都呈增加状态,在城市周边耕地减少。此外宁夏段在南部原州区、中部红寺堡地区也有明显增加。

2.3 耕地变化驱动因素分析

植被变化与降水情况息息相关,本文利用研究区内及周边46 个气象站点的降水数据进行克里金插值,得到2002−2018年研究区降水量分布图(图8),研究区内降水分布呈南多北少格局,北部地区因受气候与沙漠影响降水最少。南部甘南地区处于青藏高原东北部,山脉纵横,受地形阻挡影响水汽充足。研究区降水呈增加状态,2002−2018年平均降水量为333 mm,2002年平均降水量为310 mm,2018年平均降水量为441 mm,2018年增加量尤为显著。从图9可知,在研究时段内,大部分地区的NDVI 呈增加趋势,其中西宁−兰州−定西一带以及临夏州、白银市等地区植被呈显著增加趋势;研究区内主要山区、丘陵和宁蒙灌区的部分地区植被呈不显著增加趋势;而河湖等地植被呈减少趋势。对降水和NDVI 进一步分析(图10),研究区内降水和NDVI 都呈增加趋势,降水是影响植被生长的重要自然条件,研究区年平均降水量呈增加趋势,NDVI 也随之升高。研究区降水主要集中在6−9月,且秋收作物种植广泛,降水量增多能够为农作物生长提供有力支持,进而成为耕地扩张的动力。韩海青等[26]、张露洋[27]的研究中也指出降水是耕地变化的动力。

农田水利技术水平对耕地变化有很大影响,宁蒙灌区地势和缓,灌区内大面积采用渠道输水和地面灌溉,其灌排工程基础条件、灌溉用水习惯等在黄河流域引黄灌溉有较高代表性[28]。自1998年以来,宁蒙灌区在国家政策支持下进行了较为完备的节水改造与水利工程建设,渠道调控和灌区供水保障能力显著提升。河套灌区在1998−2015年间进行大规模节水改造与建设,河套灌区拥有渠、沟10万余条,各级灌排渠道6 万余公里,各类设施物10 万余座[29],灌溉水利用系数从2005年的0.300 提高到2019年的0.557。在宁夏灌区部分,有干渠、支干渠15 条,总长1540 km;排水干沟32 条,总长790 km;小型电力排灌站570 座,排灌机井506 多眼;各大干渠总引水能力750 m3·s−1,年引水量67 亿m3[30],宁夏灌溉水利用系数由2003年的0.350 提高到2016年的0.514。节水与水利工程建设成效明显,宁蒙灌区农田水利完善使得农业有更多灌溉水可用,这一因素促进了耕地面积的扩展。

从2002年开始,我国全面实行退耕还林政策,山区耕地逐渐减少。青海省从2000年开始试点退耕还林,2002年全面实行,到2016年退耕还林还草19.33万hm2;而甘肃的甘南州自2000年开始退耕还林到2017年末,退耕2.61 万hm2。图7 也反映出青海段(大通回族土族自治县、互助土族自治县、湟中县、湟源县)、甘肃段(甘南地区、天祝藏族自治县)等山区耕地呈减少趋势,这与政策要求相一致,说明政策的影响很显著。另外研究区内沿河流分布的耕地可以利用地表水进行农业灌溉,而像景泰、原州区、红寺堡区等地表水缺乏的地区难以进行灌溉,但灌溉工程建设(如景电提灌工程、固海扬黄灌溉工程等)对这些地区的农业耕作起到了决定性作用。如宁夏红寺堡灌区在扬水工程建设以前并没有耕地分布,在1998年工程建设后至2016年已经成为解决4万多hm2、30 万人口用水的关键。还有如甘肃景泰县内地表水资源量很少,能用于灌溉的更少,而景电提灌工程的建设为景泰县、古浪县、民勤县和阿拉善左旗提供了有效的灌溉用水,总灌溉面积达7.2 万hm2。

在经济因素上,收集2002−2018年间研究区内县级统计数据,将同一年份的数据合并得到一个总值,而后使用逐步回归分析方法进行影响因素分析,从而得到影响耕地变化的驱动因素。在参考统计资料与前人研究[31]的基础上选取以下指标:X1总人口(万人)、X2城市人口(万人)、X3乡村人口(万人)、X4城市化率(%)、X5地区生产总值(亿元)、X6第一产业生产总值(亿元)、X7第一产业结构比重(%)、X8地区人均生产总值(元)、X9农民可支配收入(元)、X10粮食产量(万t)。每年一个总量指标,共计10 个。利用逐步回归分析方法对各个指标与耕地变化进行分析,结果中只有农民可支配收入被保留,其余指标经过分析后皆被排除,因此农民可支配收入是影响耕地面积变化的主要因素。经青海地区实地调研,青海省有相当一部分耕地弃耕、撂荒,原因就是当地居民因为经营耕地所得到的收入少于外出打工,致使劳动力外流,进而导致耕地变化。

3 结论与讨论

3.1 讨论

本文使用遥感MOD13Q1 数据结合作物物候等信息提取黄河上游至河口镇的耕地分布信息。利用MODIS 数据进行耕地面积提取已经被学者广泛应用,且在区域尺度上耕地提取效果较好[32-33]。本文研究区范围较大,同样适宜使用MODIS 数据。MOD13Q1数据时间分辨率是16 d,全年连续,而Landsat 卫星数据受到一定限制,一些年份内时间连续性较差。另一面,MOD13Q1 数据的图幅较大且为处理后的产品数据,另外还有专门用于处理MODIS 数据的工具,相对于其他遥感影像数据在处理应用上更为便捷。但是MOD13Q1 数据的空间分辨率为250 m,相较于SPOT、高分系列卫星数据等高分辨率数据有一定的混合像元问题,尤其是在像青海和甘肃部分地块较为破碎的地区,耕地提取的效果可能相对较差。本文使用决策树分类进行耕地提取,受主观判断影响较大,但研究区范围较大对其内部分别进行决策树规则建立要优于统一使用单一规则。其他学者如史飞飞等[34]使用15 m 分辨率的高分遥感影像提取耕地面积的分布特征,结果精度为85.7%,但是高分影像适用于局部地区,而本研究范围和时间尺度跨度较大,使用高分影像难度较大。汤敏[35]使用随机森林算法与面向对象分类方法结合SPOT6 影像进行耕地提取,结果显示面向对象分类方法精度较高,但是面向对象分类进行分割时所需参数较多且不容易确定,需要进行多次尝试,在便捷性上稍弱于决策树方法。

耕地提取结果显示2018年在宁夏和甘肃南部地区耕地面积有明显增加。图8 显示两区域2018年降水量较往年明显增加,对甘肃、宁夏南部区域进行裁剪,分别统计各年份的平均降水量和NDVI 均值,结果如图11所示,可以明显看出2018年降水和NDVI 值要高于历年统计值。当年降水量较往年偏高,植被生长条件变好,使得当年MOD13Q1 数据的NDVI 值也较往年高,从而使得利用NDVI 提取耕地的面积增多,误差增加。另一方面,甘肃省在2007年开始进行引洮一期工程,2015年8月正式运行,为会宁、安定等7 个县(区)城乡用水提供了水资源保障,发展灌溉1.27 万hm2,受益人口达225.35 万人[36];引洮二期工程是一期工程的延伸,于2015年10月正式开工,主要解决定西、白银等8 县267 万人的用水,计划发展灌溉面积1.95 万hm2[37]。引洮工程实施使得甘肃南部地区耕地扩张。

耕地的变化还与城市发展扩张有关。耕地变化结果显示在城市的集中分布区由于城市发展占用耕地而使耕地呈减少状态。在对社会经济因素的分析上,选用指标和方法存在一定的局限性,一个是指标选取依赖于统计资料,可能会出现一些偏差情况。另一个由于收集资料有限,在指标选取上有一定的取舍。

3.2 结论

耕地是社会发展的基础资源,耕地变化会影响粮食安全,明确耕地变化对保持粮食稳定、细化用水等方面有重要影响。本文基于MOD13Q1 数据集,利用决策树分类方法并结合作物物候数据提取了黄河上游至河口镇流域的耕地面积,获得了研究区耕地时空变化情况,并对耕地变化因素进行了分析。研究结果如下:

1)研究区耕地主要沿河流、河谷分布,青甘地区地块较为破碎,多呈斑块分布;宁蒙地区地势较为平坦,多呈集中分布。城市周边和部分山区耕地减少,甘肃南部、宁夏南部、内蒙沿黄地区耕地增加较多。

2)2002−2018年研究区耕地增加88.21 万hm2,增加速率为5.18 万hm2·a−1;青海段减少5.36 万hm2,减少速率为0.20 万hm2·a−1;甘肃段增加18.89 万hm2,增加速率为1.05 万hm2·a−1;宁夏段增加29.93 万hm2,增加速率为1.87 万hm2·a−1;内蒙段增加44.74 万hm2,增加速率为2.46 万hm2·a−1。

3)影响耕地变化的因素有降水增加、节水改造、灌溉工程建设、政策因素、经济收入等,水资源政策和地方社会环境变化是研究区耕地面积变化的主要驱动因素,另外自然条件变化也趋向于促进植被生长,对耕地扩张有一定促进作用。