挠力河流域耕地变化及其对流域植被及蒸散发特征的影响

尹 剑,夏芮兹,卫丹琪,丁 乙,姜洪涛,罗心愿

(1.贵州财经大学 西部现代化研究中心,贵州 贵阳 550025;2.贵州财经大学 大数据应用与经济学院,贵州 贵阳 550025;3.东北农业大学 水利与土木工程学院,黑龙江 哈尔滨 430072)

耕地是保障区域粮食安全和维护区域经济社会稳定的重要支撑,以耕地为主的流域的生态水文环境健康是粮食生产可持续发展的重要前提[1]。挠力河流域地处三江平原腹地,是中国东北重要的粮食生产基地。改革开放40多年来,东北地区粮食生产结构在市场驱动下不断调整。不断发生湿地开垦、旱地改水田(旱改水)和水田改旱地(水改旱)的耕地变迁带来的土地利用结构变化。人类活动的频繁干扰和近年来全球气候变化的影响,给区域带来了水资源紧张和生态系统脆弱化等潜在的生态环境风险。

气候变化和人类活动所引起的生态水文效应已成为当前全球变化研究领域的焦点问题[2]。气候变化改变了与水文过程紧密相关的气候要素,从而直接影响地表水文要素[3]。人类活动所导致的流域下垫面土地利用变化,如耕地开垦、种植结构调整、退耕还林、城市化等,会改变地表反射率、植被叶面积指数等影响水循环的下垫面性质从而影响植被蒸散、下渗等产流过程,总体上影响着生态系统和水循环系统的时空分布[4]。因此,下垫面土地利用变化可以较为直接地反映人类活动的强度和影响程度[5],土地利用变化、特别是耕地的变化,常代表人类活动强度,用于分析人类干扰下生态系统和水文过程的响应[6]。

多年来,学者们对挠力河流域的耕地变化及生态水文环境开展了一系列研究。这些研究有注重分析土地覆被结构的,如李娜等[7]研究了水田进程对耕地时空变化的影响;路昌等[8]研究了土地利用变化的梯度效应,以辅助区域耕地、林地的管理;周浩等[9]研究了土地利用变化强度的空间集聚特征。另一些研究预测了气候变化和人类活动驱动下流域土地利用的变化情景。如在RCPs情境下的水土资源平衡[10-11]以及历史趋势和政策约束条件下可能的耕地规模[12-13]。土地利用变化、种植结构改变对区域生态水文产生较大的影响。基于此,一些学者对流域的生态变化、水文特征进行了分析。曾星雨等[14]基于挠力河流域归一化植被指数(NDVI)分析了流域生态系统生产力的变化情况,对各子生态系统的生产能力进行了评估。宋爽等[15-16]借助景观结构方法,研究了多指标驱动下的挠力河流域的生态健康程度。宫兴龙等[17-19]研究了土地结构变化对多种水循环变量的影响,构建相关矩阵探索了土地利用变化下水循环特征。张莹等[20]从微观尺度,检测了不同作物期间耕地需水过程。周浩等[21]围绕水田化特征,分析了流域未来水资源供给以及水土平衡下的种植规划。王梓云等[22]分析了挠力河流域湿地退化特征以及气候变化和人类活动的影响。

上述研究从多个角度对挠力河流域的土地利用、耕地规模进行了研究,但主要还是基于某一种单一地类的特征,如水田、旱地或者湿地,综合所有土地利用类型来分析流域生态水文变化的研究相对缺乏,对生态和水文特征变量研究的结合相对不足,此外研究各类耕地变化对全流域的生态水文变化的贡献有待进一步开展。研究表明归一化植被指数(NDVI)和蒸散发(ET)是反映流域生态水文的重要指标[23-24],是分析流域生态生产力和水资源承载力的关键。基于此,本文首先利用地统计学方法分析了1980—2018年挠力河土地利用特别是耕地变化及其空间分异,然后基于MODIS再分析数据,通过比较研究NDVI和ET的变化趋势,量化了其生态水文影响,最后探讨了气候变化和人类活动下,耕地变化对流域生态水文变量变化的贡献。

1 材料与方法

1.1 研究区概况

研究区(图1)位于三江平原腹地的挠力河流域(45°43′～ 47°35′N,131°31′ ～ 134°10′E),流域总面积约为26 484 km2,平均海拔60 m。属大陆性季风气候,半湿润半干旱。1月和7月的平均气温分别为-21.6、21.6 ℃。年平均降雨量可达565 mm,平均蒸发速率542.4 mm。由于地势低,地表径流少,广泛发育的洪泛平原有利于湿地的形成。流域地形呈现西南高、东北低态势,水系自西南流向东北。地貌类型主要由山地与平原组成,山地分布于流域西南部和南部,约占流域面积的38.31%,其余区域为平原。挠力河流域主要包括5类土壤,分别为暗棕壤、黑土、白浆土、草甸土、沼泽土。流域内植物资源丰富,植物群落为中温带针阔叶混交林。尤其是1980年以来,流域农业开发活动频繁,目前已建成7个现代化农场,总人口达125万人,其中农业人口占65.4%。挠力河流域是三江平原重点商品粮食基地之一,也是中国黑龙江省三江平原大面积沼泽湿地的典型区域,其粮食产量约占三江平原粮食总产量的30%。自20世纪50年代以来,流域湿地面积大规模萎缩,导致湿地生态功能下降。近年来降水量比20 a前减少了180 mm左右,减少幅度是俄罗斯远东地区的2倍。

图1 研究区概况

1.2 数据源

本研究使用的数据包括:①挠力河流域相关地形、水系等栅格、矢量图;②土地利用数据来源于中国多时期土地利用遥感监测数据集(CNLUCC)的30 m分辨率网格土地利用数据;并以Google Earth 及 Google Earth Engine (GEE)提供的数据为辅助数据源,耕地信息数据的解译精度验证则通过Google Earth 软件布控数据采样网格点实现(解译准确率均大于85%);③国家气象信息中心提供的多年来流域的气象资料;④研究时段内的空间分辨率为250 m、时间分辨率为16 d的NDVI数据,提取自MODIS的MOD13Q1产品,8 d的实际ET数据,提取自MODIS的MOD16产品;上述NDVI和ET数据已经由美国地质调查局(USGS)实现了GEE平台的共享,数据处理过程均通过GEE云平台完成。

1.3 研究方法

1.3.1耕地变化数据采集与分析

将土地利用保存为Grid文件,然后导入计算1980、1990、2000、2010、2020年的耕地(水田、旱地)面积。通过叠加两幅农田图,得到多年的农田结构变化。此外,运用经济学中的复利公式计算耕地年变化率,见式(1):

(1)

式中k——耕地面积年变化率;Ua、Ub——首、末年的耕地面积;T——研究年限。

1.3.2NDVI和ET变化率的计算

拟利用NDVI变化率评价研究区的生态状况变化。首先,对250 m MODIS的MOD13Q1产品进行重采样、投影变换、矢量裁剪等预处理后,计算每个像元生长季的平均NDVI值。流域的生长期设定在5—10月,与主要植被生长期相匹配。其次,采用Mann-Kendall(MK)趋势分析法[30-31]计算多年各像元平均NDVI的年变化率,见式(2)—(5):

(2)

(3)

(4)

(5)

式中xj、xi——第j年和第i年各像元NDVI的平均值;F——二值符号函数;n——研究周期的年数;β——NDVI的年变化率,β>0表示NDVI呈上升趋势,β<0表示NDVI呈下降趋势;Kendall-tau 相关系数(τ)是用于衡量2个变量之间的非参数相关性的指标,τ介于-1和1之间,当τ=1时,表示2个变量完全一致的排序,当τ=-1时,表示2个变量完全相反的排序,当τ=0时,表示2个变量之间没有任何相关性。

Kendall-tau相关系数的显著性可以通过Z值和p值来判断。p值表示在零假设下(即2个变量之间不存在相关性),观察到的Kendall-tau 相关系数或更极端的结果出现的概率。p值越小,表示Kendall-tau 相关系数越显著。Z表示NDVI变化的显著性水平,Z为正值表示上升趋势,负值表示减少趋势,Z的绝对值在大于等于1.645,1.96,2.576时表示分别通过了置信度90%,95%,99%的显著性检验。p值和Z值之间的关系可以通过正态分布的累积分布函数来描述。p值越小,表示观察到的结果越不可能是由于偶然因素导致的,因此假设越不可信。一般来说,p值小于0.05被认为趋势是显著的。当x表示为蒸散发变量时,可以计算蒸散发的变化趋势。

1.3.3耕地变化的生态水文影响

研究采用了基于NDVI变化幅度的农田变化生态影响的评价指标——生态影响指数。该指数定义为耕地变化引起的NDVI变化对研究区域平均NDVI变化的贡献。同时,基于年ET变化幅度的农田变化耗水影响的评价指标来评估耗水影响,该指数定义为耕地变化引起的ET变化对研究区平均ET变化的贡献。研究区的耕地变化包括3种类型:生态恢复(Ecological Restoration,ER),在研究区主要以退耕还林为主;新开垦农田(New Farmland Reclamation,FR),即研究期间新开垦的农田;耕地结构变化(Farmland Change,FC),旱地水田之间的交替变化。

对于生态恢复(ER)对区域NDVI变化的贡献估算见式(6):

(6)

式中Cer——ER对区域平均NDVI变化的贡献率;ker——退耕还林或草地栅格的平均NDVI变化率;kh——整个流域NDVI的平均变化率;Aer、Ah——ER、流域的总面积;kf——成熟森林NDVI的平均变化率,大致代表2000—2020年气候变化的影响,ker和kf之间的差距被认为是ER的贡献率。

对于FR,使用类似的式(7)来估计贡献,其中贡献率等于kfr与kn之差;kfr为新开垦农田栅格NDVI的平均变化率;kn为未改变为耕地栅格的用地类型引起的NDVI变化;Cfr为FR对区域NDVI变化的总贡献率;Afr是FR的总面积。

(7)

对于FC,研究区可能存在的有旱改水和水改旱情况,对区域平均NDVI变化(Cfc)的贡献率由式(8)估算:

(8)

式中kfc——转换后耕地类型的NDVI平均变化率;ko——转换前耕地类型的NDVI平均变化率;Afc——对应的FC类型的面积。

耕地变化对流域平均NDVI变化的总贡献率(Ch)为3种耕地变化类型对NDVI变化的贡献率之和,见式(9):

Ch=Cer+Cfr+Cfc

(9)

基于蒸散发的耗水程度以分析年蒸散发为主,同样针对全流域、旱改水、水改旱、开垦等兴趣区进行。公式类似于耕地变化的生态影响计算过程。在此略过表达。

2 结果

2.1 耕地变化及其空间分异

挠力河流域经过多年开垦,土地利用发生了较大变化,图2显示了1980—2020年40 a的土地利用变化情况,由于2010、2020年土地利用情况类似,因此只显示2020年的土地利用情况。

图2 研究区多年土地利用变化

表1所示,1980—2020年,流域耕地总面积从7 273.22 km2显著增加至15 313.62 km2,增加了110.55%,年均增加1.88%;增长最快的时段为1980—1990年,增加了35.07%,年均增加4.41%。水田面积的增加明显高于旱地的增加。40 a间,水田增加了542.14%,年均增加4.76%,其中1980—1990年增加最快,增加了57.14%,年均增加8.84%;旱地增加了52.79%,年均增加1.07%,同样也是1980—1990年增加最多,为30.26%,年均增加3.67%。旱地在2010—2020年减少了9.41%。2000—2010年,耕地面积变化不明显。

表1 1980—2020年耕地面积的变化情况

2.2 统计尺度的土地利用转移情况

为了探求各类土地利用的转换情况,研究选取20 a间隔进行了土地利用的转移探索。图3借助比例弦图,表示了各类土地利用状况以及相互之间的转换关系。方向关系为前年期向后年期变换,外围的圆弧表示后一年期的土地类型在流域中的面积占比,联系条带表示土地利用类型之间相互转换的关系。互相有联系条带表示存在转换关系,宽度越大表示对应的转换量越大,当联系两端宽度近似时,说明彼此之间转换面积近似,但具体的区域位置不同。

图3 各类土地利用转移比例

1980—2000年的土地利用变化可以结合表1和图3a得出。1980年水田面积相对较小,2000年水田的面积发生了较大的变化,1980—2000年增加了242.38%,旱地也增加了60.98%。1980—2000年有大量的湿地和草地转换为旱地和水田,约有45.33%的湿地消失,主要转换为耕地,占湿地比例的43.20%。草地减少了80.86%,其中73.74%的草地转移为耕地。林地主要转换为耕地,占原有林地面积的13.49%,主要转换为旱地。总体来说耕地的增加来自于湿地、林地和草地等自然生态用地的开垦,湿地、草地、林地大面积减少。

对于水田来说,除了保留了1980年自身的92.93%面积外,其主要来源为草地、湿地和旱地的转换,分别占到2000年当期水田比例的21.77%、27.76%和21.89%。旱地的主要转换来源为林地、草地和湿地,分别占其总面积的12.27%、15.41%和17.70%。1980—2000年存在一定比例的旱改水,其中10.03%的旱地变为水田;水改旱相对较小,占水田面积的6.91%,只占到旱地总面积的0.57%,总体来说水田-旱地之间的转换不明显。

2000—2020年土地利用类型变化情况见图3b。这20 a间,湿地、草地和林地进一步退化,转换为水田和旱地。有53.46%的湿地转化为耕地,其中18.48%的湿地转换为旱地,35.27%的湿地转换为水田。39.36%的草地转换为耕地,8.56%的林地转换为耕地。旱地和林地存在面积近似相等的互相转化。此外,该时段内,旱地和水田存在较大面积的互转,但旱地转向水田的面积大于水田转向旱地的面积。旱地发生了显著的改水田变化,有31.60%的旱地转换为水田,占2020年水田面积的59.17%,也就是说超过半数的新增水田来自于旱地的转移。而也有47.41%的水田转换为旱地,占到2020年旱地面积的14.22%。由于一定比例的湿地转换为旱地补充了相应的旱改水面积缺失,占2020年旱地面积的12.76%,使得旱地总体面积变化不大。

纵观1980—2020年,较大面积的林地、草地和湿地转变为耕地,是旱地的主要来源。湿地向水田旱地的转换最多,1980—2020年40 a间湿地向旱地的转移显著大于前2个20 a。水田面积增加除了来自于前三者的转移外,旱改水占据较大的比例。在旱改水的过程中,同样进行着水改旱的过程。虽然图3b 显示2000—2020年水改旱较为明显,但是1980—2020年水改旱在图3c中呈现出不明显的情况,可能的原因是在2000—2020年水改旱过程中,转出的水田很大比例上是1980—2020年湿地转换为水田的部分。为了探求这一原因,需要借助网格尺度的土地利用转换分析。

2.3 斑块尺度的土地利用转移情况

为进一步分析旱改水和总体土地利用变化,图4显示了总体1980—2000、2000—2020、1980—2020年的土地利用变化在斑块尺度上的地图,图例表示了靠前年份的土地覆被类型向靠后年份的土地覆被类型的转移情况,若两者一致表示未发生转换。为了便于显示主要的耕地和自然生态用地之间的转移,转移面积较小的情况在途中略去。由图4a可以看出湿地转水田的情况在1980—2000年十分明显,这部分主要集中分布在支流交汇地区,而这一区域在2000年后进行了集中的水田改旱地工程,在图4b中由水田转换为旱地。图4a 和图4b反映出上游山地是林地开垦为农田的集中区。图4c则反映出中下游地区湿地垦荒区域,最终转化为旱地和水田的交错带。2000—2020年发生了显著的水改旱过程,北部支流沿岸大量的旱地转换为水田。

图4 挠力河流域土地利用转移

在当前保障国家粮食安全的前提下,现有农田的保护是十分重要的,进一步将评估近20 a来,旱改水、水改旱以及生态用地退化所产生的生态影响。2000—2020年的土地利用类型转移中,面积排名前十的兴趣区单元分别为:不变林地(28.48%)、不变旱地(23.37%)、旱地→水田(12.23%)、不变水田(5.42%)、水田→旱地(5.27%)、不变湿地(5.00%)、湿地→旱地(4.70%)、林地→旱地(2.54%)、湿地→水田(2.46%)、旱地→林地(1.70%),共占流域总面积的91.16%。基于此,研究在图4的基础上,进行了进一步的分区,获得研究感兴趣的旱地、水田、旱改水、水改旱、湿地转水田、湿地转旱地、林地转旱地,退耕还林这几类占比较大的转换地区,由于研究范围较大,相对来说建设用地、水体、草地的占比较小,研究暂不讨论,重点讨论与耕地相关的区域。

2.4 耕地变化对区域NDVI值的影响

为了保障粮食安全,在现有耕地的基础上,短期内不会发生较大的土地利用性质的改变,而调整区域种植结构是改善区域生态安全和水资源安全的有效手段。基于此,重点研究耕地结构变化显著的2000—2020年,探求耕地结构变化是否给区域生态带来影响。

首先研究NDVI的变化情况。图5表达了全流域的生长季各月(5—10月)、生长季平均、全年平均和最大值NDVI年际变化。从各月NDVI变化的情况来看,除了10月变化不显著外,全流域NDVI在生长季均呈现微弱的上升趋势,由于大量的湿地转换为水田和旱地,湿地的植被密度不及耕地,因此在生长季,总体的生长季NDVI得到提升。8月是各类植被生长最为茂盛的时期,不论哪种植被在此段时间叶面积较大,因此上升趋势相对其他月份变化不大。6、7月是变化最大且最显著的月份,可能的原因是不断地湿地开垦,增加耕地,农作物的分布密度远远大于湿地,在生长季提升了全流域的NDVI。整体生长季NDVI来看,也有显著的上升趋势,同时,全年的NDVI极大值亦有上升趋势,全年总体平均值没有明显的变化。

图5 全流域生长季各月逐年NDVI变化

针对不同土地利用变化引起的NDVI变化,图6分析了各类土地利用兴趣区域中的NDVI均值,在20 a间的变化情况。图中的林地、旱地、水田和湿地均值未发生土地覆被类型改变的区域,下同。图6主要表达了各月林地、湿地2种土地类型不变区域的NDVI,不变旱地、不变水田区域的NDVI,以及林地转旱地、林地转水田、水田转旱地、旱地转水田等集中用地类型转换,但最终是耕地的土地利用类型的NDVI。

图6 各类兴趣区生长季各月逐年NDVI变化

研究发现,除了8月以外,林地在各月一直是NDVI的高值区。5、6月湿地的NDVI总体高于其他耕种区域,7—9月湿地的NDVI整体低于耕地,10月湿地耕地高于耕地。林地改旱地地区的NDVI值小于林地,且在作物生长初青期(5、6月)逐渐与林地的NDVI值差距增加;作物生长期(7、8月)与林地NDVI的差值逐渐减少,成熟期(9、10月)又呈现出差距增加,这说明林地改旱地是一个不断增加的过程,直至2020年完成所有兴趣区的土地利用转换。各类耕地兴趣区的NDVI相差不大,整体来看,水田NDVI在所有耕地中较大,旱地改水田的NDVI逐渐与水田接近,其他几种向旱地转换的耕地兴趣区的NDVI逐渐与旱地接近。

2.5 耕地变化对区域ET的影响

进一步研究年蒸散发量在各类土地利用变化兴趣区中的变化过程情况(图7)。各类用地均呈现一定的蒸散发上升趋势,且趋势均通过显著性检验。林地、湿地等自然生态用地上升趋势较小,显著性和相关性也较小,上升趋势分别为2.56、2.93 mm/a。各类耕地均有显著的上升趋势和较高的上升率,这表明20 a来耕地开垦导致流域蒸散发耗水的增加,这一问题值得重视。总体来看,水田及各类转向水田的用地的蒸散发大于旱地和其他用地转向旱地的兴趣区,而旱地和转旱地兴趣区的蒸散发耗水增加速率高于水田。在林地和湿地等自然生态用地增加率低的前提下,可以推断,耕地蒸散发受气候和人类活动影响增加较大,在未改变种植方式的前提下,气候变化更容易影响土地类型为耕地的耗水。从流域整体的蒸散发变化情况来看(图7k、7l),挠力河流域蒸散发增加显著性和增加量都大于自然生态用地(林地和湿地),呈现显著的增加趋势,蒸散发增加量为4.492 8 mm/a,林地和湿地的蒸散发增加量只有2.75 mm/a。

图7 挠力河流域各类用地年蒸散发量变化趋势

2.6 耕地变化对区域ET的影响

进一步研究各类土地利用变化兴趣区对总体NDVI和ET变化的影响。将几类主要土地转换类型进行汇总,计算2000—2020年生长季年平均NDVI相对变化率和全年ET的相对变化率(表2)。由于林地开垦变化没有通过显著性检验,因此暂时不显示相关数据。

表2 2000—2020年的生长季NDVI和全年ET变化情况 %

从表2可以看出,ET的增加率整体高于生长季NDVI的增加率,耕地在蒸散发上均有显著增加,且旱地大于水田。旱地的NDVI增加率也是最大的。湿地的NDVI减少,可能与近年来三江平原湿地退化有关。

进一步对各类兴趣区进行汇总,得出人类活动干扰耕地对区域NDVI和ET的贡献(表3)。对于NDVI来说,新开耕地对总体NDVI的提高产生负向贡献,即开垦农田影响了流域NDVI提升,对生态健康造成负向影响;退耕还林对总体NDVI提升产生正向贡献;旱改水同样产生正向贡献。总体耕地变化对NDVI提升的贡献率之和为-18.56%,可以证明当前人类活动强度对区域生态健康提升有一定的阻碍。对于ET,旱改水的ET增加贡献有9.13%,加速了流域蒸发耗水;水改旱贡献为负,说明通过水田改旱地工程,对流域水资源起到一定的保护作用;湿地开垦贡献较大,说明湿地在转换为耕地的过程中,增加了耗水。

表3 2000—2020年的各类兴趣区对生长季NDVI和全年ET的贡献情况

3 讨论与结论

3.1 讨论

通过对结果的分析,可以发现,耕地面积的增加主要来自于林地和湿地的开垦。首先,尽管耕地的增加有助于提高粮食产量,但是这也会导致生态环境的恶化。因此,需要采取措施保护湿地和林地等自然生态用地。其次,流域内近20 a发生了较为显著的水田和旱地的互转,在面积上旱改水大于水改旱。这说明农民在种植作物时更多地选择水稻,这可能与水稻的经济效益有关。但是,水田和旱地之间的转换也会对流域生态环境造成影响,需要进一步探究如何平衡经济效益和生态保护。第三,耕地区域NDVI的增加主要集中在作物生长季节,相反在非生长季对NDVI变化贡献低。这表明耕地对植被覆盖的影响主要集中在作物生长季节。因此,在种植作物时需要采取措施促进作物生长,同时也要注意非生长季节的土壤保水和养分保持。最后,各类耕地变化对流域NDVI和年蒸散发的贡献率不尽相同。退耕还林对保护生态和减少耗水有贡献;湿地开垦虽然增加了生长季NDVI但也增加了流域耗水;水改旱有利于流域保水,而旱改水工程增加蒸发耗水,可能带来水资源安全风险。因此,需要根据不同类型的土地利用变化采取相应的生态保护和水资源管理措施。

3.2 结论

研究基于云平台遥感数据对中国粮食重要产区的挠力河流域进行了耕地干扰下的生态水文关键变量NDVI和ET的分析,获得了近年来该地区耕地变化对流域植被生态和耗水量的影响。

40 a间,大量的湿地、草地转换为旱地,部分林地也被开垦为旱地,水田的增加主要来源于湿地的开垦和2000—2020年的旱改水,流域的耕地面积在40 a间增加了110.55%,其中水田面积增加542.14%,旱地增加了52.79%。2000—2020年旱改水和水改旱同时发生,且集聚区不同,旱改水发生在流域北部,水改旱发生在流域中部。

结合遥感云平台得出,2000年以来,全流域及各类与耕地相关的土地类型的生长季NDVI和年ET均有不同程度的增加趋势。在不变林地NDVI增加的同时,湿地NDVI减少,表明近年来流域湿地生态存在一定风险。流域NDVI在全年无明显变化。各类耕地变化对流域NDVI和ET变化的贡献率不尽相同,总体上退耕还林对保护生态和减少耗水有贡献;湿地开垦虽然增加了生长季NDVI但也增加了流域耗水;水改旱有利于流域保水,而旱改水工程增加蒸发耗水,可能带来水资源安全风险。未来将进一步探究耕地种植结构的细部变化,分析影响生态水文的详细原因和驱动机制。如对耕地利用变化的内在机制进行深入探究,以便更好地理解耕地利用变化对流域生态水文的影响;进一步探索如何平衡经济效益和生态保护,以促进可持续发展。