华北平原蒸散发变化及对植被生产力的响应*

王林娜,韩淑敏,李会龙,杨永辉**

(1.中国科学院遗传与发育生物学研究所农业资源研究中心/中国科学院农业水资源重点实验室/河北省节水农业重点实验室石家庄 050022; 2.中国科学院大学 北京 100049)

蒸散发(evapotranspiration,ET)是联系大气和陆面水文过程的关键环节,对流域水循环和水量平衡具有重要影响,然而由于蒸散发过程涉及辐射、传导、扩散、对流和表面大气相互作用等多个过程,是区域水分循环中最难估算的变量,准确评估陆面蒸散发时空格局变化特征,不仅有利于理解其对水循环的影响及作用机制,对解决水资源危机和完善水资源管理体制也具有重要意义。

传统的蒸散发研究主要基于站点特定条件下的观测,难以反映大区域复杂下垫面的蒸散发过程,而遥感手段使得蒸散发估算由站点尺度扩展到区域尺度。随着遥感技术的不断成熟,衍生出众多遥感蒸散发产品,如PML_V2、MOD16A2、GLEAM、SSEBop、NTSG 等,正在成为全球和区域尺度蒸散发和水循环研究的有利手段。2019年发布的Penman-Monteith-Leuning Evapotranspiration V2(PML_V2)产品(2000-2020年),因具有较高时空分辨率(8 d/500 m),在全球10 种植被类型和95 个通量站点进行了验证,在流域水循环与水资源管理领域具有广泛的应用前景。Zhang 等利用PML_V2 模型和哨兵2号10 m 土地利用数据,估算了北京市不同土地覆盖类型的月蒸散量,为政策制定和城市用水规划提供了有价值的科学依据。Li 等基于PML_V2 数据定量评估了植被变化对黄河流域水储量的影响,发现2003年以来黄河流域植被变化加剧了区域水储量的下降,需要采取一定的措施平衡生态与水文的关系。Luan 等基于PML_V2 模型将黄土高原窟野河流域的植被变化和煤炭开采对径流的影响进行了剥离,量化了两种人类活动对径流减少的贡献率,为制定合理的煤矿开采政策提供了科学依据。

受空间分析手段限制,以往蒸散发研究大多以区域尺度定性描述和变化分析为主。近年来,基于栅格尺度的趋势变化和显著性检验为蒸散发的定量分析和敏感区识别提供了新手段,可有效分辨栅格尺度上蒸散发的变化趋势及变化的显著性,筛选植被生产力、土地利用/土地覆盖等对蒸散发影响的关键因子,为区域生态系统管理提供更细致有效的科学指导。已有学者将Theil-Sen Median 趋势分析方法和Mann-Kendall 显著性检验方法用于栅格尺度的蒸散发分析,研究了不同区域蒸散发的时空变化特征,如Ghalami 等在伊朗、Li 等在淮河流域、谷佳贺等在黄河流域的研究等都佐证了该方法的可行性。

华北平原是我国粮食主产区,水资源短缺一直是限制其农业发展的重要因素,区域中有70%的农业用水依靠地下水开采,这也导致该地区地下水位下降严重。近30年来,一系列的节水灌溉措施有效缓解了华北平原水资源匮乏的现状,但依旧无法解决地下水位不断下降的困局。基于高精度遥感产品研究蒸散发的时空演变特征,对于揭示区域的灌溉节水效果、制定完善的水资源管理对策具有重要意义。近年来,有学者先后对华北平原的蒸散量变化及影响因素开展定性研究,如Yang 等基于1 km 空间分辨率的ET-watch 蒸散发和降水数据,估算了2002-2007年海河平原灌溉需水量,并识别出在太行山前一带由于农业集约化种植对水资源可持续利用造成威胁。Gao 等基于模型方法估算了海河流域1960-2002年的实际蒸散量,发现年降水量和潜在蒸散量的降低是导致流域实际蒸散量减少的主要原因。截至目前,基于栅格数据围绕华北平原开展的遥感蒸散发研究仍显不足,需要利用高空间分辨率产品,分辨蒸散发变化的敏感区和主要影响因素。

本文基于高精度PML_V2 遥感产品,结合Theil-Sen和Mann-Kendall 趋势分析方法,研究2000-2019年华北平原栅格尺度蒸散发的时空变化格局,识别蒸散发变化的演变趋势、变化显著性及成因,为完善水资源管理对策提供科学参考。

1 数据与方法

1.1 研究区概况

华北平原东临渤海,西依太行山,北靠燕山,南抵黄河,地势低平,海拔多在50 m 以下,面积约13.8万km。地处暖温带湿润或半湿润气候区,雨热同期,年均温11～12 ℃,年降水量一般在400～600 mm,季节分配不均,春季常有较严重旱情发生,夏季伴有洪涝灾害。主要农作物包括小麦()、玉米()、棉花(spp.)、大豆()等,种植模式主要为冬小麦-夏玉米一年两作。

农业用地和建设用地是研究区最主要的两种土地利用类型,所占比例分别为74.3%和20.9%(2020年),本文以这两种最主要的土地利用为例,并选择典型区域进行蒸散时空变化特征研究。针对农业用地,以石家庄和保定(山前平原区)、衡水(中部平原区)、德州(东南部黄河灌区) 3 类共4 个农业生态类型区为例; 建设用地以北京、天津、石家庄和保定4 个大中型城市为例(图1)。

图1 华北平原和典型研究区位置Fig.1 Locations of the North China Plain and typical areas

1.2 数据来源及处理

本文研究时段为2001-2019年,遥感蒸散发数据集和总初级生产力(GPP)数据来源于中国科学院地理科学与资源研究所张永强研究员及团队开发的全球PML_V2 陆地蒸散发与总初级生产力数据集(PML_V2); 归一化植被指数(NDVI)数据来源于中国科学院资源环境科学与数据中心提供的中国月度1 km NDVI 数据集; 其中NDVI 和GPP 用于表征华北平原的植被生产力; 土地利用数据来源于我国自然资源部制作的GlobelLand 30(表1)。

表1 数据类型、名称和来源及时空分辨率Table 1 Data type,name and source,and temporal and spatial resolution used in the study

不同时间尺度(年、小麦生长季和玉米生长季)的蒸散发数据在GEE(Google Earth Engine)平台由日数据累加完成,其中小麦生长季为每年的10月至次年的5月,玉米生长季为每年的6月至9月。在进行栅格尺度相关分析时,由于PML_V2 蒸散发的分辨率为500 m,而NDVI 的分辨率为1 km,故在Arcgis 10.4 平台将NDVI 的分辨率重采样至500 m。

在分析农业用地和建设用地蒸散量年际变化时,由于蒸散数据和土地利用数据的空间分辨率不同,在空间上叠加计算时会产生差异,为减少误差,在栅格尺度上将PML_V2 数据与GlobelLand30 数据进行交集运算,将500 m 分辨率栅格内农业用地大于80%的设为农业用地像元,建设用地大于60%的设为建设用地像元,以此来反演不同土地利用的蒸散值。

1.3 研究方法

Theil-Sen Median 斜率和Mann-Kendall 检验方法均为稳健的非参数方法,且很少受到异常值的干扰,在生态水文领域广泛用于判断长时间序列变量的变化趋势和显著性检验研究。本研究首先采用Theil-Sen Median 斜率法计算2001-2019年栅格尺度蒸散量的斜率 β。斜率 β表示蒸散量的变化趋势,当＞0 时,蒸散量呈上升趋势; 当=0 时,蒸散量变化趋势不明显; 当＜0 时,蒸散量呈下降趋势,计算公式为:

式中: Median 为中值函数,1＜＜＜,x、x分别表示第、年的蒸散量。

Mann-Kendall 方法用于研究时段内栅格尺度蒸散量的显著性检验,计算公式为:

其中:

式中:为正态分布函数,均值为0,方差var()=(-1)(2+5)/18。当＞0 时,标准的正态统计变量通过下式计算:

式中:为统计值,大于0 时呈增加趋势,小于0 时呈减少趋势。的绝对值大于1.645、1.960 和2.576时,分别表示通过了置信度90%、95%和99%的显著性检验。

以上公式利用MATLAB 语言在像元尺度上实现,输入项为2001-2019年蒸散量值,输出项为研究时段内的斜率值和显著性检验值。

利用MATLAB 语言,在像元尺度上分析蒸散量与影响因子NDVI 和GPP 的相关程度,采用最常用的皮尔逊相关系数法,计算公式为:

显著性采用检验,检验相关的可信程度,当＜0.01 时为极显著性水平,当＜0.05 时为显著性水平。

2 结果与分析

2.1 华北平原蒸散量的时空分布特征

从图2 可知,华北平原2001-2019年蒸散量为542～665 mm,年际间呈震荡波动态势,平均为588.1 mm。趋势分析结果表明,研究时段内年尺度蒸散量总体呈不显著上升趋势,上升速率为0.65 mm·a;2008年和2015年蒸散量距平变率相对较大。小麦季蒸散量为182.2～258.1 mm,平均235.2 mm,总体呈现不显著下降趋势。玉米季蒸散量为329～415 mm,平均377.8 mm,总体呈现不显著上升趋势。玉米季蒸散量占全年的64.2%,小麦季蒸散量占全年的35.8%。

图2 2001—2019年华北平原全年和不同作物生长季的蒸散量(ET)变化特征Fig.2 Evapotranspiration(ET) of year and crop growth seasons in the North China Plain from 2001 to 2019

华北平原多年平均蒸散量和小麦生长季的蒸散量在空间上的分布格局极为相似(图3a,b),主要表现为: 在东南部黄河灌区及保定、石家庄山前平原带出现高值区,在邢台和邯郸东部等传统棉花种植区出现低值区,低值中心出现在渤海湾一带及北京、天津、石家庄等大中型城市,此外水域蒸散量大都大于1500 mm; 玉米季蒸散量(图3c)在大部分地区无明显的空间差异,主要分布在350～450 mm,城镇地区蒸散量值较低,高值区与全年蒸散量高值区的分布相吻合,主要分布在东南部黄河灌区一带,表明多年平均蒸散量高值主要由玉米季贡献。

图3 华北平原全年及不同作物生长季蒸散量及其显著性变化趋势的空间分布Fig.3 Spatial distribution of evapotranspiration(ET) and its’ significant change trends in the North China Plain for annual and different crop growing seasons

不同时间尺度蒸散量的变化趋势和显著性检验的空间分布特征如图3d-f 所示。年均蒸散量的检验结果为: 约有75.6%面积呈增加趋势,其中34.6%呈显著性增加,主要分布在东南部的黄灌区和中部的衡水和沧州等地; 24.4%的区域呈下降趋势,主要分布在太行山前平原一带,仅有2.8%区域达显著性下降。小麦季蒸散量有59.3%的区域呈现下降趋势,16.3%区域显著减少,主要在太行山前平原呈条带状分布; 有40.7%的区域呈现上升趋势,其中14.9%区域呈显著上升,主要分布在黄河灌区一带。玉米季蒸散量绝大部分区域呈现上升趋势,占比达89.9%,其中呈显著上升趋势占比为46.4%,主要分布在石家庄、保定太行山前平原区、华北平原的北部及整个东南部黄河灌区。

2.2 不同土地利用类型蒸散量与趋势变化

基于华北平原农业用地和城市建设用地占比的绝对优势,以这两种土地利用类型为例,阐明蒸散量的年际变化特征及趋势变化。从图4 可知,农业用地年际蒸散量为562～693 mm,多年平均值为590 mm,其中小麦季为214 mm(占36%),玉米季为376 mm(64%)。4 个典型农业区中,德州的蒸散量最高,保定和衡水较低,与图3a 的整个华北平原蒸散量的空间分布特征相一致。另外,几个典型区域的蒸散量年际变化趋势均在2015年出现拐点,2015年前均呈波动增加趋势,之后呈波动减少趋势(图4a)。华北平原农业用地蒸散量的趋势变化以增加为主,有85.5%面积呈上升趋势,其中有42.3%达显著上升,主要分布在东南部黄灌区一带,而蒸散量减少区域很少,零星分布于山前平原城市周边地带(表2)。

图4 2001—2019年华北平原及典型区不同土地利用类型蒸散量(ET)变化Fig.4 Evapotranspiration(ET) under different land use types in the North China Plain and typical areas from 2001 to 2019

表2 2001—2019年华北平原不同土地利用类型蒸散变化显著性区域面积比例及主要分布区Table 2 Significant area proportions and main distribution areas of evapotranspiration change under different land use types in the North China Plain from 2001 to 2019

华北平原城市用地的年际蒸散量为481～614 mm,多年平均值为526 mm(图4b)。4 个典型城市蒸散量2008年和2015年出现峰值,2015年后迅速下降,年际变化整体呈下降趋势,表明城市化过程中由于建筑用地增多和耕地的减少,导致不透水面增多,蒸散量减少。不同城市多年平均蒸散量大小为: 天津(467 mm)＜北京(487 mm)＜石家庄(511 mm)＜保定(557 mm),这一定程度上表明城市化水平越高的城市,建筑用地面积占比偏高,蒸散量偏低。在变化趋势上(表2),城市用地年均蒸散量呈下降与上升趋势的面积占比分别为50.9%和49.1%,下降趋势略高于上升趋势,其中呈显著上升和下降变化的面积占比均在10%左右。值得注意的是,在北京、天津主城区蒸散量呈现上升趋势,而在城市外围区呈下降趋势,而其他大部分城市内部区和外围区均呈下降趋势。

2.3 蒸散量与植被生产力的相关性分析

2001-2019年华北平原蒸散量与植被生产力(NDVI 和GPP)的逐像元相关性分析如图5 所示。蒸散量与NDVI 相关系数介于-0.92～0.96,整个华北平原平均相关系数为0.21,绝大部分地区蒸散量与NDVI 呈正相关,面积占比达76.5%(图5a),其中呈显著正相关的面积占比为23.6%,主要分布在太行山前平原的石家庄、保定一带及黄河灌区的聊城、德州一带(图5b),这些区域是冬小麦-夏玉米主要种植区;呈负相关的区域占比为23.5%,主要分布在邢台和邯郸东部及京津唐地区,其中仅有1.5%的区域蒸散量与NDVI 呈显著负相关,主要分布在北京和天津一带。

图5 2001-2019年华北平原蒸散量(ET)与归一化植被指数(NDVI)、总初级生产力(GPP)的相关系数及显著性相关的区域Fig.5 Distributions of correlation coefficients and significance between evapotranspiration(ET),gross primary productivity(GPP)and normalized difference vegetation index(NDVI) in the North China Plain from 2001 to 2019

华北平原蒸散量与GPP 相关系数介于-0.91～1.00,整个华北平原平均相关系数为0.54,有87.6%的区域蒸散量与GPP 呈正相关,几乎覆盖了整个区域,相关性更强的地区主要分布在东南部黄灌区一带(图5c); 达到显著正相关的区域占到整个华北平原的61.5%(图5d),反映了区域内作物产量增加对蒸散发的正向反馈,尤其在种植冬小麦-夏玉米一年两作的黄河灌区及太行山前平原区地带,这种正向反馈尤为突出。

典型农业区耕地的蒸散量与植被生产力的变化趋势如图6 所示。2001-2019年4 个典型农业区蒸散量、GPP、NDVI 的年际波动基本趋于一致。石家庄、保定和衡水3 个区域的蒸散量、GPP 和NDVI 总体增长趋势不明显,未达显著水平,但在2015年后有不同程度的下降趋势(石家庄和衡水NDVI 除外); 导致蒸散量与植被生产力下降的原因可能与河北省逐步推广实施的地下水压采政策有关,2015年后冬小麦种植面积逐步减少。在整个研究时段,德州的蒸散量、GPP 和NDVI 均达到显著上升(＜0.05),表明该区域的蒸散量、植被生产力增加较为明显。

图6 2001—2019年华北平原4 个典型区域蒸散量(ET)、总初级生产力(GPP)、归一化植被指数(NDVI)的年际变化Fig.6 Annual evapotranspiration(ET),gross primary productivity(GPP) and normalized difference vegetation index(NDVI) in four typical agricultural areas of the North China Plain from 2001 to 2019

2001-2019年4 个农业区耕地的蒸散量与GPP、NDVI 的相关系数均为正,反映了植被生产力对蒸散发的正向反馈(图7)。蒸散量与GPP 的相关性在德州地区呈极显著相关(0.61,＜0.01),在衡水地区呈显著相关(=0.23,＜0.05)。蒸散量与NDVI 的相关性只在德州表现出显著相关(=0.32,＜0.01)。总体上,蒸散量受GPP 的影响更大,尤其是黄灌区的德州地区和低平原的衡水地区,表明蒸散量增加与农田增产有关; 石家庄和保定山前平原地带蒸散量与GPP 和NDVI 均未达显著相关水平,表明该区蒸散发受植被生产力等多因素共同影响。

图7 华北平原4 个典型农业区蒸散量(ET)与归一化植被指数(NDVI)、总初级生产力(GPP)的相关性Fig.7 Correlations between evapotranspiration(ET) and normalized difference vegetation index(NDVI) / gross primary productivity(GPP) in four typical agricultural areas of the North China Plain

3 讨论

3.1 PML_V2 适用性

本研究采用的PML_V2 为8 d/500 m 高精度遥感产品,在全球和区域尺度都得到了广泛应用。何韶阳等为选择更适用于华北平原的高分辨率蒸散量产品,对PML_V2、SSEBOP、MOD16A2 等3种高精度遥感产品在华北地区进行精度验证,通过与5 个涡度相关测量数据对比,发现PMLV2 产品在华北地区精度最高。研究时段内(2001-2019年)华北平原年蒸散量的波动范围为542～665 mm,该结果显著低于张喜英、Liu 等在中国科学院栾城农业生态系统试验站分别利用长期定位试验取得的年蒸散量为777～834 mm(1980-2017年)和850 mm的观测结果,但与Yang 等基于类似的遥感蒸散发模型ETwatch 估算的海河平原2002-2007年蒸散量波动范围527～679 mm 相近。这可能与遥感蒸散发产品单像元内(0.5 km×0.5 km)土地覆盖类型较复杂、灌溉农田占比不高有关,降低了单像元内的年蒸散量。

研究时段内小麦生长季节华北平原的蒸散量多年平均值为220.5 mm,其中山前平原和东南部黄灌区蒸散量为200～300 mm,北部和滨海地区较低,不足200 mm,而豫北地区较高,为300～400 mm。该结果与莫兴国等基于PM 法和双源-PT 法发展的简化遥感蒸散模型获取的遥感蒸散发数据,分析得到的有关山东及河北中东部农业区蒸散量为250 mm左右,滨海区不足200 mm 的研究结论一致,也与吴喜芳等、方蓓婧在华北平原估算的结果相似。玉米季大部分区域蒸散在300～450 mm,这也与莫兴国等得出的华北平原玉米季大部分区域蒸散在360～420 mm 结果相似。

3.2 华北平原蒸散量时空分布格局与成因

华北平原年均蒸散量整体以年际波动为主,增加趋势不明显。其中小麦季蒸散量特别是2015年后呈波动下降趋势(图2),一定程度上与2014年中央一号文件提出的“先期在华北地区河北省开展地下水超采综合治理试点工作”有关,农业逐步实施的高效节水灌溉工程、水源置换工程、种植结构调整、地下水压采、冬小麦季节性休耕等一系列措施,缓解了华北平原蒸散量的上升趋势。小麦生长季有59.3%的区域蒸散量呈下降趋势,主要分布在河北省境内,反映了该省实施冬小麦休耕、压采等政策的效果。

为进一步明确玉米季蒸散量的变化趋势,图8 分析了玉米季各像元尺度2001-2019年时间序列GPP与蒸散量的相关性及结果的显著性,两者相关系数达到显著性的区域与玉米季蒸散量显著升高的区域(图3f)有较高的一致性,表明玉米季蒸散发增加与GPP 或秋季作物生产力和产量增加有关。

图8 2001-2019年华北平原玉米季蒸散量(ET)与总初级生产力(GPP)的相关性及显著性Fig.8 Correlation and significance between evapotranspiration(ET) and gross primary productivity(GPP) for corn growing season in the North China Plain from 2001 to 2019

华北平原蒸散量的空间分布格局及趋势变化有明显的异质性(图3)。在北部地区和山前平原城市群一带蒸散量较低,在500 mm 左右,这与周琳估算的2003-2012年北京城市蒸散量为517 mm 接近,由于城市用地下垫面多为不透水层,蓄水能力差,蒸腾作用弱,年蒸散量低; 但值得注意的是北京、天津大城市内部蒸散量明显增加,这与2014年国家提出的建设“海绵城市”政策有关,城市绿地公园、透水路面、绿化屋顶等增多,取代了传统大城市中心的不透水地面,蒸散量随之增加。邢台、邯郸东部出现一个低值中心,原因是该区域为河北省主要产棉区,受一年一季耕作方式影响,蒸散量较低,趋势变化不明显。沿渤海地区蒸散量一般在500 mm 以下,该区域遍布盐碱地及盐碱土,农业灌溉较少,导致蒸散量也较低。2013年以来“渤海粮仓”工程的实施,大量低洼盐碱开发为耕地,对农业增产、增收效果显著,导致低平原地区蒸散量呈增加趋势。在华北平原东南部的黄河灌区一带,蒸散量呈现高值区,一般在650～750 mm; 近年受经济效益驱动和劳动力成本上涨的影响,冬小麦-夏玉米的种植模式逐渐取代一年一季棉花种植,复种指数增加也会造成区域蒸散发增加,因此黄灌区蒸散量呈显著增加趋势。

3.3 植被生产力对蒸散量的影响

蒸散量的时空变化受气候变化和人类活动等多个因素影响,其中植被生产力增加是影响陆地蒸散发的主要驱动因素,植被变绿解释了全球84%的蒸散发变化。Jin 等在黄河流域的研究也表明一半区域内的蒸散量变化与植被变绿相关,探索植被变化与蒸散发的关系成为揭示区域蒸散发和水循环变化的关键科学问题。本文通过研究华北平原ET与NDVI、GPP 的空间相关性,发现分别有76.5%和96.0%的面积呈正相关关系,其中23.6%和68.1%的区域内相关性达显著水平,需要说明的是,虽然呈显著正相关,但两者的相关系数并不高,故一种因子不能全部解释影响蒸散的因素,但在一定程度上仍能表明以灌溉农业为主的华北平原,与产量相关的GPP更能反映作物种植强度和植被生产力对蒸散发的影响。

针对华北平原3 类农业类型区蒸散发与植被生产力关系的研究表明,以德州为代表的黄河灌区蒸散量受植被生产力的影响显著,蒸散量与NDVI、GPP 均显著相关,而两者相关性越高则代表两者的变化趋势趋于一致,表明在东南部黄灌区一带,与GPP 相关的复种指数增加和农田增产等使植被冠层吸收更多的辐射,植被蒸腾速率和蒸散发量也随之提高。以衡水为代表的中部低平原区蒸散量受GPP 影响显著,与GPP 的相关系数为0.48,达到显著相关; 这与中低产田改造、冬小麦及蔬菜、水果种植占比增加有较大关系。而石家庄、保定为代表的山前平原区,是传统的一年两作高产区,蒸散量与植被生产力的关系均不显著,与NDVI、GPP 的相关系数仅在0.30 左右,与冬小麦休耕、大面积发展苗圃、药材种植等减慢了蒸散发和GPP 等的增长速率有关。

4 结论

本研究基于高精度蒸散发产品PML_V2 探究了2001-2019年整个华北平原、主要土地利用类型及4 个典型农业区的蒸散量的时空分布特征及变化趋势,解析了蒸散量时空差异的成因以及植被生产力对蒸散量的影响,主要结论如下:

1)华北平原年均蒸散量为588.1 mm,整体呈波动增加趋势,且有明显的季节性差异,具体表现为玉米季(377.8 mm)＞小麦季(235.2 mm)。空间尺度上,华北平原75.6%的区域年蒸散量呈增加趋势,东南部黄灌区增加趋势最显著; 玉米季89.9%的区域蒸散量呈增加趋势,但小麦季59.3%区域蒸散量呈下降趋势,这与近年来逐步实施的地下水压采和小麦休耕等政策有一定关系。

2)植被生产力对蒸散发的影响显著,ET 与NDVI、GPP 都呈正相关关系,但ET 与GPP 的相关性更高,说明在耕地占绝对优势的华北平原GPP 更能反映作物种植强度和复种指数等对蒸散发的影响; 不同农业类型区植被生产力对蒸散量的影响程度不同,以德州为代表的黄灌区蒸散量受植被生产力影响更大。

3)不同土地利用类型蒸散量空间分布特征和变化趋势差异显著。农业用地内蒸散量均值为590 mm,农业用地蒸散量总体呈上升趋势; 建设用地蒸散量总体呈下降趋势,但在北京、天津大城市中心呈上升趋势,这可能与海绵城市建设有一定关系。

感谢Google Earth Engine 平台提供的PML_V2数据,感谢中国科学院资源环境科学与数据中心提供的NDVI 数据,感谢国家自然资源部提供的30 m土地利用数据。