1. 引言
近年来,随着全球气候变化以及不断加剧的人类活动,导致土地资源退化成为对社会经济系统可持续发展和人类生存环境构成严重威胁的重大全球性环境问题之一 [1],并受到国际组织与国家的广泛关注 [2]。联合国可持续发展目标跨机构专家组(IAEG-SDGs)依据UNCCD关于土地退化定义 [3],确定土地利用/覆盖变化趋势、土地生产力变化趋势和碳存量(土壤碳)变化趋势为识别土地退化与否的核心指标 [4]。土地生产力是土地重要功能之一,既能反映土地资源质量状况又能表征土地的生产能力,其变化能有效反映土地退化状态特征。在SDG15.3.1元数据表 [5] 以及《SDG15.3.1良好实践指南(2017)》 [6] 中指出植被指数(直接观测得到)及衍生数据被认为是土地生产力分析中独立、稳健的土地生产力数据,是土地生产力的常用替代指标数据 [7] [8] [9]。同时,大量研究表明NDVI能够有效反映与评估植被生长状况,且受气候影响显著。因此,通过分析NDVI时空变化规律及其自然驱动因素,探讨土地生产力时空演变趋势及对气候的响应,对评价区域土地资源质量、预测未来发展与土地退化治理具有现实意义。
中巴经济走廊(China-Pakistan Economic Corridor, CPEC)是连接中国与阿拉伯海的重要陆上走廊,是“一带一路”倡议的旗舰项目。CPEC巴基斯坦段主要位于干旱和半干旱地区,生态环境非常脆弱,地形复杂,气候多样 [10] [11]。此外,近年来随着气候变化以及人口迅速膨胀,不合理的人类活动加剧,如森林砍伐和过度放牧,过度垦荒,加之工业化和城市化导致生态环境恶化、土地退化、森林退化、土壤退化和土壤侵蚀等日益严重 [12] - [17],造成了巨大的经济和环境损失 [14] [15]。Ullah等 [14] 分析巴基斯坦马尔丹区土地退化程度发现:土壤侵蚀、过度放牧、人口快速增长、农业扩张和土地管理不善是引起该区域土地退化加剧的主要原因。Muhammad等 [15] 发现土地退化对农业可持续发展和粮食安全产生潜在威胁。Niazi [16] 认为不平衡的土地所有权以及耕地不均衡分布导致过度使用土地,最后造成土地生产力下降。目前针对土地生产力的研究,缺乏对直接观测的时间序列遥感数据的有效应用,以遥感技术为主的对地观测技术以其全覆盖、可重复、客观性强的优点,获取的时间序列遥感数据更有效评估土地生产力时空变化并分析其变化的驱动因素。
基于此,本文利用NDVI数据以及降水、气温、水汽压等气象数据,采用Theil-Sen Median (Sen)趋势分析、Mann-Kendall (MK)显著性检验、Hurst指数以及相关分析等方法,逐像元分析,研究巴基斯坦地区2000~2020年土地生产力退化情况、预测未来变化趋势,以及巴基斯坦地区土地生产力对气候的响应。以期为巴基斯坦区域生态环境的恢复与建设、可持续土地资源管理、荒漠化与土地退化防治等提供参考。
2. 研究区
巴基斯坦位于南亚次大陆(23˚30'~36˚45'N和60˚53'~75˚31'E之间) (图1),国土面积约88 × 104 km2 (含巴控克什米尔地区),下辖4个省和4个特区,东北、北、西、东部分别于与中国、阿富汗、伊朗和印度接壤,
注:该图基于自然资源部标准地图服务系统下载的审图号为:GS (2016)第1665号的标准地图制作,底图无修改。
Figure 1. The geographical location of study area
图1. 研究区示意
南濒阿拉伯海。巴基斯坦地势总体上由西北向东南逐渐降低,地形复杂,北部和西部以高原、山地为主,中东部为平原,南部以沙漠为主,高山和丘陵占全国面积的一半以上。巴基斯坦气候,从南到北大致分划为热带沙漠气候、热带季风气候和高山山地气候,境内大部分地区处于干旱和半干旱状态,全年高温,整体降雨少、干湿两季分明且空间分异显著,主要分布在北部山区以及中东部地区。巴基斯坦是一个典型的农业国,超过62%的人口(2017年总人口超过2亿)直接或间接以农业为生,其境内最重要的粮食产地是印度河平原。近二十年来,巴基斯坦城镇化和工业化快速推进,经济发展与粮食安全、生态安全之间冲突日益突出,土地荒漠化现状严峻。
3. 数据与方法
3.1. 数据
本文根据SDG15.3.1中元数据表 [5] 及UNCCD《SDG15.3.1良好实践指南(2017)》 [6] 中关于土地生产力指标评估数据的说明以及结合巴基斯坦具体实际情况,选择NDVI数据土地生产力评估数据,降水、平均气温(气温)、最低气温、潜在蒸散、水汽压以及最高气温等气象数据,作为本文研究基础数据。
NDVI数据源于美国国家航天局(NASA)的MODIS-NDVI产品数据(MOD13Q1, https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD13Q1),时间分辨率为16 d,空间分辨率为250 m,时间序列为2000年1月~2020年12月,由GEE平台集成处理,采用最大合成法(MVC) [18] 获取月尺度数据,进一步求平均获取年尺度数据。降水、水汽压气温等气象数据选用CRU Ts v.4.05 [19] (https://crudata.uea.ac.uk/cru/data/hrg/),空间分辨率为0.5˚ × 0.5˚,时间序列为2000年1月~2020年12月,对月尺度数据采用双线性插值法重采样与NDVI数据匹配,后进一步通过求均值、最大值、最小值等合成所需年尺度气象数据。
3.2. 方法
本文利用NDVI数据作为土地生产力评估数据 [20] [21],采用Sen趋势分析与MK显著性检验相结合逐像元分析,研究巴基斯坦二十年来土地生产力的时空变化特征,在此基础上结合Hurst指数探讨土地生产力变化趋势的持续性特征。同时,利用相关分析法分析气候数据与土地生产力数据之间相关程度,以此探讨巴基斯坦土地生产力变化对气候的响应特征。
3.2.1. Sen趋势分析与MK显著性检验
Sen与MK相结合的趋势分析与显著性检验方法,由于数据不需要服从一定分布、能有效抵抗误差数据、对离群数据不敏感、以及具有坚实的统计学理论基础,其结果科学性和可信度高 [22],被广泛应用于NDVI时间序列分析中 [8] [23] [24]。借助MATLAB软件,采用Sen方法对2000~2020年NDVI进行逐像元变化趋势分析,并采用MK方法对变化趋势进行显著性检验,将二者结果叠加,综合分析巴基斯坦土地生产力时空变化情况。
利用Sen趋势分析方法计算NDVI变化趋势 [8] [25],公式为:
(1)
式中,
为NDVI Sen变化趋势,
和
分别为j年和i年NDVI的值。当
时,NDVI增加,反之则减少。
利用MK非参数据检验方法判断NDVI变化趋势显著性 [26] [27],计算公式为:
设定:
(2)
定义Z统计量为:
(3)
其中:
(4)
(5)
(6)
式中,统计量Z取值范围为
,
和
分别为研究区像元j年和i年的NDVI值,f为符号函数,n表示NDVI序列长度。采用双边检验,在给定显著性水平
下,当
时,表示研究序列数据在
水平上变化趋势具有显著性,反之则趋势不显著。
参照相关文献 [8] [28] 等的NDVI趋势分类,由于
的情况较少,选取
为退化区,介于
之间为稳定区,
为改善区。选取
(
),即
时NDVI变化趋势显著。以此将土地生产力变化水平定义为5级:明显改善(
,
,)轻微改善(
,
),稳定(
,
),轻微退化(
,
),严重退化(
,
)。
3.2.2. Hurst指数
Hurst指数是定量描述时间序列长度依赖性的有效方法,被广泛应用于气候学、植被变化等研究中 [8] [28] [29] [30] [31]。Hurst指数计算方法较多,本文采用较为常用的Hurst指数计算方法,R/S分析法 [31] [32],公式如下:
对于时间序列
,
,定义该时间序列:
1) 均值序列
,
(7)
2) 累积离差
,
(8)
3) 极差
,
(9)
4) 标准差
,
(10)
对于比值
,如
,则时间序列存在Hurst现象,H为Hurst指数,而H值由
利用最小二乘方法拟合得到。根据的H值可以判断NDVI时间序列变化趋势的持续性。一共分为三种形式:
,表明时间序列具有持续性,未来与过去变化趋势一致,且H越接近1,持续性越强;
,表明时间序列为随机序列;
,表明时间序列具有反持续性,未来与过去变化趋势相反,且H越接近1,反持续性越强。
3.2.3. 相关分析法
本文采用应用最为广泛的Person相关系数来衡量各气候因子与NDVI之间的相关程度,两个变量相关系数计算如下:
(11)
式中,
为x与y的相关系数,n为时间序列长度,
是时间;
代表NDVI的观测值,
代表气候因子的观测值,
和
为两个变量的平均值,
为相关性系数,
时,二者相关性为正,
时,两个变量呈负相关,
越大,相关性越强。本文基于MATLAB软件,计算时间序列中NDVI与各气候因子的相关系数。根据相关系数显著性检验表可知,在0.05和0.01的显著性水平下,相关系数临界值分别为0.433和0.549。基于此,本研究将相关性分为六种情况:负相关显著(
),负相关较显著(
),负相关不显著(
),正相关不显著(
),正相关较显著(
),正相关显著(
)。在此基础之上,统计显著与较显著(
)的栅格占研究区总面积的百分比,作为NDVI与气候因子的相关度。
4. 结果与分析
4.1. 2000~2020年巴基斯坦NDVI与气象因子时间变化特征
以年为时间尺度,基于年平均NDVI,分析21年来NDVI变化情况(图2)。从图2看出,巴基斯坦年均NDVI值呈波动上升趋势增速为0.022/10a (P < 0.001, R2 = 0.7937),波动范围为0.15~0.25,仅在2007年和2018年出现较大波动,2000年NDVI值为0.157,为21年来最低值,2020年NDVI值为0.217,为21年来最大值。2000~2020年全境年平均NDVI增长率为37.9%。降水、水汽压、潜在蒸散等气象数据在年尺度上基本成正玄变化趋势,年降水的波动程度最大,变化范围为185.30 mm (2018年)~391.98 mm (2020年),21年平均值为293.07 mm,整体呈上升趋势,增速为3.15 mm/a。潜在蒸散与降水变化趋势相似,变化范围为49.92~51.81 mm,增速为2.77 mm/10a。六项气象因子中仅水汽压呈微弱下降趋势,变化范围为12.8~14.5 hPa,增速为−0.087 hPa/10a。
4.2. NDVI时空变化特征
利用2000~2020年平均NDVI数据,计算21年平均值得到NDVI空间分布(图3),从中可以看出,巴基斯坦高NDVI值地区主要分布于印度河流域以及中东部平原,北部山区向中部平原过渡地区为巴基斯坦主要森林覆盖地区,NDVI值最高;印度河流域及中东部平原为农业集中区,农业依赖灌溉,受干旱影响较小,NDVI较高。西部高原高海拔地区土地贫瘠,主要为裸土、稀疏植被等,北部山区主要为稀疏植被、冰川/永久积雪等,以及南部沙漠地区NDVI值较低。NDVI值<0.1的地区主要为裸土、沙漠和冰川/永久积雪。

Figure 3. Spatial distribution of annual average NDVI in Pakistan
图3. 巴基斯坦多年平均NDVI空间分布
由图4、表1可以看出,巴基斯坦大部分地区土地生产力得到改善(71.81%),其中明显改善与轻微改善的区域分别占51.62%和15.19%,明显改善区主要分布于北部山区向中部平原过渡地区、中东部平原以及印度河流域地区,轻微改善地区多分布于西部高原与印度河平原过渡地区以及南部沿海附近区域。土地生产力退化区占3.82%,其中严重退化与轻微退化的区域分别占1.12%和2.70%,严重退化区主要分布于南部印度河入海口以及印度河流域,轻微退化区主要位于别部山区以及南部印度入海口附近地区。北部山区向中部平原过渡地区以及印度河流域地区土地生产力出现一定退化现象。基本稳定区域为24.37%,主要位于西部荒漠以及北部山区,该区域土地贫瘠,土地生产力极低。
巴基斯坦NDVI的Hurst指数阈值为0.10~0.99,均值为0.51,研究区H < 0.5的面积占比为46.71%,H > 0.5的面积占比为53.29%,近一半区域NDVI变化趋势具有反持续性(H < 0.5),表NDVI变化趋势具有波动性,结果与NDVI年际变化(图2)呈波动上升趋势相符。将NDVI空间变化与Hurst指数结果叠加,判断土地生产力变化趋势的持续性(图3),持续退化区域占1.97%主要沿着印度河流域呈零星分布;由退化到改善区域占1.84%,主要位于北部山区向中部平原过渡地区;基本稳定区占18.35%,主要位于西部高原和北部山区,西部高原与北部山区土地贫瘠土地生产力相对稳定;由改善到退化区占34.34%,广泛分布于北部山区向中东部平原过度以及印度河流域地区,该区域为巴基斯坦境内土地利用以及人类社会经济活动热点区域,过度使用耕地、过度放牧、过度垦荒等对土地生产力影响强烈;持续改善区面积占比最大为37.47%,主要分布于西部高原与印度河流域平原过渡地区;变化趋势不确定区占10.53%,主要位于西部高原和北部山区等裸地、山地等土地贫瘠、生态环境恶劣地区。

Figure 4. Land productivity change and sustainable spatial distribution in Pakistan from 2000 to 2020
图4. 2000~2020年巴基斯坦土地生产力变化以及持续性空间分布

Table 1. Statistics of NDVI change trend from 2000 to 2020
表1. 2000~2020年NDVI变化趋势统计
4.3. 相关性结果与分析
本文采用气候因子与NDVI相关度作为衡量NDVI对气候因子的响应程度,相关度越大,则该气候因子对NDVI影响越显著。
研究区内土地生产力与降水、水汽压、潜在蒸散等气候因子相关性分析结果(图5、表2),可以看出降水是巴基斯坦地区最显著的影响因子,相关度为53.43%,且降水与土地生产力呈正相关(98.92%),降水对土地生产力有明显的促进作用。西北部地区高原地区干旱严重,降水对土地生产力促进作用十分显著,而位于印度河流域的耕地,由于绝大多数处于灌溉区,土地生产力受降水影响相对较弱,北部山区多为无植被覆盖区,土地生产力极低,降水对其土地生产力促进作用有限。潜在蒸散、水汽压、气温等对土地生产力影响相对较弱。总体来看,降水是影响巴基斯坦地区土地生产力最主要的气候因子,其次为潜在蒸散水汽压、最高气温、最低气温,平均气温影响最弱。

Figure 5. Spatial distribution of correlation between NDVI and climate factors in Pakistan
图5. 巴基斯坦NDVI与气候因子的相关性空间分布

Table 2. Correlation between NDVI and climate factors in Pakistan/%
表2. 巴基斯坦NDVI与气候因子相关度/%
5. 结论与讨论
5.1. 讨论
CPEC巴基斯坦段地形复杂、气候多变,全年高温,降雨少、干湿两季分明且空间分异显著,土地荒漠化和土地退化现象严重,研究NDVI的时空变化特征,分析土地生产力变化情况,有助于促进区域生态环境保护和巴基斯坦荒漠化与土地退化防止。本文基于MODIS-NDVI数据与CRU再分析气象数据,结合Sen趋势分析 + MK非参数据统计的显著性检验、Hurst指数、相关分析等方法,逐像元分析,探讨了巴基斯坦地区土地生产力时空变化特征、变化趋势持续性特征以及自然因素作用下巴基斯坦土地生产力的响应特征。研究结果能够在一定程度上反应巴基斯坦土地生产力变化情况以及对气候因素的响应特征,有助于巴基斯坦土地荒漠化防治及区域生态环境保护。但必须注意到巴基斯坦土地生产力退化及土地荒漠化等,不仅仅由于自然因素驱动,由于其生态脆弱,人类活动对土地生产力退化有严重影响。巴基斯坦四分之三的土地贫瘠,且是一个典型的农业国,超过62%的人口(2017年总人口超过2亿)直接或间接以农业为生 [33]。农业严重依赖灌溉,灌溉导致土地内涝和盐碱化,同时耕地分布不均,导致过度使用土地情况严重,造成土地生产力下降 [16]。人类活动以及气候变化对巴基斯坦土地生产力退化影响、不同耕作或植被对地区的适应性值得更进一步研究。
5.2. 结论
1) 从时间上看,2000~2020年巴基斯坦年均NDVI值呈波动上升趋势,波动范围为0.157~0.217,增速为0.022/10a,2018年之后NDVI呈现明显的快速增长趋势。从空间上看,巴基斯坦NDVI主要呈现印度河流域及中东部平原地区高,西部高原、北部山地以及西南沙漠地区低的分布特征。
2) 巴基斯坦2000~2020年土地生产力改善地区远大于退化地区,其中明显改善与轻微改善的区域分别占51.62%和15.19%,严重退化与轻微退化的区域分别占1.12%和2.70%,基本稳定区域为24.37%。巴基斯坦NDVI的Hurst指数,近一半区域NDVI变化趋势具有反持续性(H < 0.5),表明NDVI变化趋势具有波动性。土地生产力持续退化区域占1.97%、由退化到改善区域占1.84%、基本稳定区占18.35%、由改善到退化区占34.34%、持续改善区面积占比最大为37.47%、变化趋势不确定区占10.53%。
3) 降水是巴基斯坦地区最显著的影响因子,对土地生产力有明显的促进作用。西北部地区高原地区干旱严重,降水对土地生产力促进作用十分显著,而位于印度河流域的耕地,由于绝大多数处于灌溉区,土地生产力受降水影响相对较弱。其次为潜在蒸散水汽压、最高气温、最低气温,平均气温影响最弱。