近20年黄土高原NDVI/气象要素对实际蒸散发的影响研究

韩宇平，徐 丹，黄会平，贾冬冬，裴 铎，张俊峰

(1.华北水利水电大学水资源学院，河南 郑州 450046；2.河南省黄河流域水资源节约集约利用重点实验室，河南 郑州 450046；3.华北水利水电大学建筑学院，河南 郑州 450046)

近年来，气候及植被覆被变化引起的水文效应显著增强，进而影响到区域内蒸散发，并引起区域储水量变化[1-2]。蒸散发(Evapotranspiration，ET)包括土壤、水面蒸发和植被蒸腾，是地表水文循环过程的重要组成部分，对保持区域能量及水量平衡具有重要意义[3-4]。黄土高原是中国典型的生态脆弱带，绝大部分属于干旱半干旱地区，1999年以来开始实施退耕还林(草)工程，目前植被在很大程度上得到恢复，植被覆盖度大幅增加[5-6]。植被的生长会改变地表粗糙度和地表反射率，影响汇流产生过程及蒸腾过程，使区域ET发生变化，进而对区域生态水文过程产生影响[7-9]。因此，研究变化环境下ET的时空分布规律，分析并探讨其与各影响因子的关系，对于了解区域水文条件、提高水资源利用效率有一定的指导意义[10]。

遥感数据产品具有多时间尺度、精度高及监测范围广的特点，可提供多时间、空间尺度、长时间序列的精准数据[11]。目前遥感已成为研究区域ET、植被覆盖度、区域水量存储等的重要手段[12]。基于Penman-Monteith公式研究发布的全球陆地蒸散发数据集MOD16已通过全球通量塔数据验证，模拟精度可达86%[13]，被广泛应用于各区域的蒸散发研究中。郭晓彤等[14]基于MOD16/ET数据产品，分析了淮河流域2000—2014年蒸散发的时空分布特征及对土地覆被类型变化、气候变化的响应。何韶阳等[15]选择了华北平原PML_V2、SSEBop_V4和MOD16A2 3种高分辨率ET产品，对不同产品进行精度验证和时空对比，以期选择出适合研究区域的数据产品。黄土高原因地表覆被变化剧烈成为研究的热点地区，韩盟伟等[16]基于Penman-Monteith模型计算潜在蒸发量，研究黄土高原1959—2015年ET时空分布规律及与气象因子的相互关系。Ma等[17]利用黄土高原2001—2015年地表反射率(MOD09GA)、地表温度(MOD11A1)数据，基于能量平衡方程反演ET，并进一步探讨了ET年际时空变化及与气象因子、植被变化的关系。周志鹏等[18]基于MOD16A2数据产品，分析了黄土高原全年及四季ET时空分布特征。退耕还林还草工程使黄土高原地区大量坡耕地转换为耕地和林地，地表植被度提高25%左右，随着植被覆盖度提高，区域ET显著增加[19-20]。大量研究已表明植被生长发育对ET有明显影响，但已有研究较少从植被生长的不同阶段来探讨ET的时空变化特征及影响因子，归因分析中也缺乏综合考虑各因素在空间上对ET的影响。因此本文采用MOD16A2 ET、MOD13Q1 NDVI和气象数据，分析黄土高原地区蒸散发在全年、植被生长季和非生长季3个时间尺度的时空变化特征，并在此基础上研究NDVI、气象要素对蒸散发的综合影响。研究成果有助于掌握黄土高原3个不同时间尺度蒸散发的时空规律及影响因子，以期更好地探讨植被生长发育与ET的关系，为黄土高原生态恢复和水资源综合治理、高质量发展提供理论依据。

1 材料与研究方法

1.1 研究区概况

黄土高原(33°43′～41°16′N，100°54′～114°33′E)位于中国西北部、黄河中上游地区，高程为81～4 945 m，总面积约64.9万km2，主要包括山西、陕西、甘肃、青海、宁夏、河南和内蒙古自治区7个省(区)的部分地区。该地区为温带大陆性季风气候，全年降水主要集中在7—9月，年降水量在150～750 mm之间[21]，降水量由东南向西北减少。1999年起，国家开始实行大规模的退耕还林还草工程，黄土高原为首批试点区域之一。

1.2 数据来源

ET(MOD16A2)和NDVI (MOD13Q1)数据来源于美国国家航空航天局NASA网站(https://ladsweb.modaps.eosdis.nasa.gov/ search/)，时间序列为2001—2020年，时间分辨率分别为8、16 d，空间分辨率分别为500、250 m。气象数据(降水、温度、风速)来源于中国气象科学数据共享网(http://cdc.nmic.cn/home.do)，时间序列为2001—2020年。Landsat-8影像数据来源于美国地质勘探局(USGS)(https://earthexplorer.usgs.gov/)，时间序列为2013—2020年，空间分辨率为30 m。

首先利用MRT(MODIS Reprojection Tool)软件对下载的数据进行拼接转换等预处理工作，再利用ArcGIS软件裁剪出黄土高原图，ET图像中的无效值赋值为0。黄土高原55个气象站点的降水、温度和风速等数据，利用创建的泰森多边形(Voronoi diagram)求平均值，采用克里金(Kriging)插值得到空间分布。研究区位置及气象站点的空间分布见图1。

1.3 研究方法

1.3.1空间回归分析

黄土高原植被生长发育特点为每年4月返青发芽，10月份枝叶枯黄凋萎。因此，按照植被生长时段规律划分为生长季(4—10月)、非生长季(1—3、11—12月)。利用slope趋势分析[10]，对每个栅格的ET数值进行线性回归分析，计算公式如下：

(1)

式中Eslope——统计时间尺度内的变化速率，Eslope>0表示20年间ET为增加趋势，Eslope<0表示ET为减少趋势，Eslope=0表示没有变化；N——总年数，本文为20；Ei——ET第i年的值。

1.3.2Morlet小波分析

Morlet小波分析[22]被广泛应用于水文过程中，是一种揭示周期性变化规律的研究方法。数学表达见式(2)：

(2)

1.3.3主成分分析

主成分分析(PCA)法[23]是一种把各变量通过正交变换重新组合成一组线性不相关变量的多元统计方法。通过降维技术对多个变量进行处理，在保留原信息的基础上更加客观地筛选出主导的影响因子，根据其特征值及贡献率确定主成分。本文基于PCA分析方法提取各主成分以反映不同因子的综合效应对ET时空变化的影响。在p=0.05 显著性水平下，通过方差最大化正交旋转，提取前2个主成分。并根据主成分初始因子载荷量，计算各线性组合中系数，确定各因子在综合得分模型中的系数并进行归一化处理[24]，得到各因子指标权重。

1.4 MODIS数据验证

选取2013—2020年行列号为127/35、130/35及126/33的Landsat-8影像数据，每年选择2幅图像，每幅影像时间间隔约为半年，基于SEBAL模型，结合洛川、榆中、五寨3个地面气象站点的风速、气温等气象数据，反演各个气象站点相应时间的ET值并扩展至8 d，利用Pearson相关系数法对MODIS ET数据进行精度验证(图2)。结果表明，在95%置信区间相关系数为0.88，证明MODIS ET数据在黄土高原是适用的，可以进行实际蒸散发的相关研究。

图2 SEBAL模型反演ET与MODIS ET关系

2 结果与分析

2.1 ET的时空变化特征

2.1.1ET时间序列特征

a)年际线性变化。2001—2020年黄土高原地区ET年际变化趋势见图3。黄土高原全年、生长季和非生长季ET总体均呈增长趋势，增速分别为6.780 mm/a(p<0.01)、6.256 mm/a(p<0.01)和0.524 mm/a(p<0.05)。3个时段ET数值差异较大，但变化趋势基本一致，其中生长季与全年的波动变化具有高度一致性。ET年均值为279.5 mm，最大值和最小值分别出现在2018年(347.8 mm)和2001年(191.1 mm)；ET生长季均值为225.6 mm，占全年的80.71%，最大值和最小值分别出现在2018年(286.2 mm)和2001年(135.2 mm)；非生长季ET均值为57.8 mm，占全年的19.29%，最大值和最小值分别出现在2017年(71.9 mm)和2010年(45.4 mm)。由此可知，黄土高原ET增速由快到慢依次为全年>生长季>非生长季，且非生长季与全年变化波动的一致性小于生长季与全年变化波动的一致性，其原因可能是全年的蒸散发主要集中在生长季，因此全年ET主要是受生长季的影响较大，非生长季影响相对较小。

图3 2001—2020年黄土高原不同时间段ET时间序列变化

b)年际周期性变化。图4为2001—2020年黄土高原全年、生长季和非生长季3个时段的小波系数实部和小波方差。小波系数正值代表ET数值增大，负值代表ET数值减小。由图4a、4c、4e可知，3个时段的系数实部等值线在4～8、13～16 a时间尺度上比较密集。由图4b、4d、4f可知3个时段均出现2个峰值，全年及生长季峰值分别出现在6、14 a的时间尺度，其中14 a出现的峰值大，说明全年及生长季时段14 a左右的周期震荡最强，为ET变化的第一主周期，6 a为ET变化的第二主周期；非生长季2个峰值分别出现为5、14 a，最大峰值为14 a，即非生长季14 a为ET变化的第一主周期，6 a为ET变化的第二主周期。不同时段小波方差大小为全年>生长季>非生长季，即在年尺度上黄土高原ET周期变化特征最明显。根据14 a强主周期特征，2020年处于ET值减小阶段，预测2020年之后会进入增加阶段。

图4 黄土高原不同时段小波系数实部图和小波方差

图4 黄土高原不同时段小波系数实部图和小波方差

c)突变检测。区域年均ET的M-K检验见图5a，在2001—2020年，年均ET持续增加(UF>0)；ET突变点发生在2010年附近(99%的置信区间，u0.01=±2.32)，表明2010年后，ET增速加快。图5b为年均NDVI突变检验，研究时段内，NDVI持续增加(UF>0)。NDVI的持续增长对ET产生正相关作用，ET在2010年附近发生突变，原因可能是退耕还林还林工程实施以后，树木生长需要一定时间，约10年左右树木成林，因此2010年以后成年树木对黄土高原ET影响更大，加快ET的增长速度。

图5 突变检验

2.1.2空间变化特征

图6为2001—2020年黄土高原3个不同时间段的ET空间分布。全年及生长季均由东南向西北递减，高值位于黄土高原东南部的陕西省部分区域；东部和南部的山区、灌溉区等地次之；低值主要位于黄土高原西北部。非生长季高值主要分布于南部秦岭山地和西部贺兰山部分区域。

图6 2001—2020年黄土高原不同时段ET空间分布

由表1可知，年尺度上，ET值在200～300 mm之间的面积占黄土高原总面积的31.48%，100～200、300～400 mm的比例分别为23.29%、23.83%，极小区间(0～100 mm)和极大区间(600～800 mm)所占比例仅为4.16%、0.21%。生长季100～200 mm之间的面积占研究区总面积的比例最大为33.59%，200～300、300～400 mm的区间占比分别为30.90%、19.87%，极小区间(0～100 mm)和极大区间(600～700 mm)所占比例仅为9.99%、0.01%。全年和生长季ET在空间上有明显的分异特征，但其极值占比很小，90%左右的区域ET为100～500 mm。非生长季ET数值整体偏小，94.63%的区域小于100 mm，100～200 mm的区域仅占5.37%。

表1 黄土高原不同时段ET多年均值占比

为进一步研究ET时间变化趋势的空间分布特征，对每个像元进行线性回归分析，得到2001—2020年3个时间段ET变化斜率(图7)。由图可知，全年及生长季ET变化空间分布情况较为一致，东南部及东部增长幅度较大；非生长季西部地区增长幅度较大。3个时段呈增加趋势的面积不同，全年及生长季增加趋势面积占比分别为93.34%和93.48%，非生长季为75.33%，非生长季变化幅度远小于生长季和全年。

图7 2001—2020年黄土高原不同时间段ET变化斜率

2.2 NDVI及气象要素时空变化

2.2.1NDVI时空变化特征

由图8知，全年、生长季和非生长季的NDVI整体呈增长趋势，增速分别为0.005 /a(p<0.01)、0.007 /a(p<0.01)和0.003 /a(p<0.01)，可见退耕还林(草)政策实施以后，黄土高原地区生态环境改善，植被覆盖率逐年上升。20年间，不同时段NDVI均值为0.305(全年)、0.387(生长季)和0.189(非生长季)，NDVI最大值均在2020年，其值分别为0.463(生长季)、0.219(非生长季)和0.362(全年)，NDVI最小值均在2001年，其值分别为0.305(生长季)、0.153(非生长季)和0.242(全年)。黄土高原不同时段NDVI的平均值为生长季>全年>非生长季，且3个时段NDVI时间序列变化基本一致。

图8 2001—2020年黄土高原不同时间段NDVI时间序列变化

受水热条件的影响，3个时段多年平均NDVI空间分布范围基本一致，由东南向西北递减，东南部值最高，其次为东部，西北部值最低(图9)。西南地区植被类型主要是针阔叶林、寒温性灌林和草丛，而西北地区受水热条件限制，植被类型多为草原、草甸和高山稀疏植物，因此西北部NDVI值最低[25]。全年、生长季和非生长季NDVI空间分布规律与其对应时段ET空间分布基本一致。

图9 2001—2020年黄土高原不同时间段NDVI空间分布

2.2.2气象要素时空变化特征

20年间黄土高原年平均降水量为444.99 mm，生长季平均降水量为406.30 mm，非生长季仅为39.12 mm。由图10a可知，3个时段的降水量均呈增加趋势，全年、生长季、非生长季降水量增速分别为3.820 mm/a(p<0.05)、3.295 mm/a(p>0.05)、0.410 mm/a(p>0.05)。图10b表明全年、生长季和非生长季的平均气温分别为9.14 ℃、16.50 ℃、 -1.24 ℃，增长速率分别为0.015、0.010、0.014 ℃/a，但均不显著(p>0.05)。图10c表明全年、生长季和非生长季平均风速增速分别为2.14、2.16、2.13 m/s。20年间黄土高原全年、生长季和非生长季平均风速增速分别为0.011、0.010、0.013 m/(s·a)，增长速率较低，但平均风速年际间波动较大。

图10 黄土高原不同时段降水、温度、风速变化

由图11可知，3个统计时段内不同气象要素空间分布规律一致，降水和温度由东南向西北递减，与ET分布规律基本一致；风速由东向西递减，东北部风速最高。

图11 黄土高原不同时段降水(mm)、温度(℃)、风速(m/s)空间分布

2.3 ET影响因子分析

2.3.1ET影响因子的主成分分析

选择NDVI、气温、降水、风速作为影响因子，运用SPSS 23.0 提取主成分特征值大于等于1.0的主成分(表2)。生长季、非生长季和全年的前2个主成分累计贡献率分别为79.11%、76.19%和83.60%，这表明前2个主成分基本能够反映这3个不同时间尺度NDVI与气象因子对ET的贡献度。

表2 特征值、贡献率、累计贡献率

由表3可知生长季第1主成分特征值为1.958，贡献率为48.95%，NDVI、降水和风速荷载较大，均呈正相关；第2主成分特征值为1.207，贡献率为30.16%，温度在第2主成分中有较大荷载，且呈正相关。非生长季中，风速、NDVI和温度在第1成分中占比较大，特征值为1.874，贡献率为46.86%，第2成分特征值为1.173，贡献率为29.33%，降水有较高荷载。在年时间尺度上，第1主成分特征值和贡献率分别为2.161和54.04%，温度和风速有较大荷载，第2主成分降水和NDVI为主要荷载，其特征值和贡献率分别为1.183和29.56%。

表3 旋转成分矩阵

各因子指标权重计算结果见表4。生长季时段，风速、气温和NDVI所占权重较大，降水权重比例较低，其指标权重大小为风速>气温>NDVI>>降水，说明在生长季时段对ET影响有较大贡献的是风速、气温和NDVI。非生长季中，NDVI、降水和气温权重占比较大，风速权重比例相对较低，其指标权重大小为NDVI>降水>气温>风速。在年尺度上，NDVI、降水和风速权重比例较高，气温权重比例较低，其权重大小为NDVI>降水>风速>>气温，即在全年中，黄土高原NDVI和降水对ET的贡献较大，其次为风速，气温对黄土高原ET影响并不明显。

表4 主成分因子赋权计算

2.3.2各主成分因子空间分布

利用ArcGIS10.2软件Multivariate中的Principal Component功能，取2个主成分分析各自空间分布(图12)。生长季第1主成分由西向东、由北向南逐渐增大，南部和东北部最大，即NDVI、降水和风速对东部、南部影响大，对西部、北部影响较小。第2主成分为温度，对东南部影响较大，对南部影响较小。非生长季第1主成分从西北、西部向东南部逐渐升高，即风速、NDVI和温度对东南部地区影响较大，对西北、西部地区影响较小。第2主成分为降水，高值呈散点状分布在东北及南部和西南小部分地区。全年尺度第1主成分温度和风速由西北向东南逐渐升高，即温度和风速对ET影响为东南高、西北低。第2主成分由东南、西北部向外逐渐减小，即降水和NDVI对东南和西北地区影响较大，对东部、西部及东北部地区影响较小。

图12 ET各主成分空间分布

3 讨论

黄土高原自1999年实施退耕还林(草)工程以来，植被覆盖度明显提高，伴随着全球气候变化，区域的实际蒸散发量及水文条件都发生了明显改变。本文研究结果表明实际蒸散发以6.256 mm/a的速度显著增加，年均ET值自东南向西北递减，这与韩盟伟等[16]、周志鹏等[18]研究结果基本一致。不同时段NDVI及气象要素对黄土高原ET的影响程度不同，但均呈正作用。在干旱半干旱的黄土高原地区，区域蒸散发受植被密度和降水影响较显著[2,26]，区域NDVI、降水的空间分布规律与ET高度一致，且在年尺度上黄土高原NDVI和降水对ET的贡献较大，表明NDVI、降水对ET的影响较大。植被大规模恢复及气候变化加剧了蒸散发，增加了黄土高原水资源需求量，进而影响到水资源承载力[27]。影响ET的因子除NDVI、气象要素外，还包括海拔[28]等自然因素以及土地利用变化[8]、修建水工建筑物[29]等人类活动的影响。后续研究将考虑其他因素对实际蒸散发的影响，并进一步分析实际蒸散发量变化对区域水文环境及陆地水储量的影响。

4 结论

基于MOD16A2 ET、MOD13Q1 NDVI和气象数据，采用线性回归分析、Morlet小波分析和主成分分析等方法，探讨了黄土高原2001—2020年ET时空演变特征及NDVI和气象因子的影响，结论如下。

a)黄土高原全年、生长季和非生长季ET均增加，增速为全年(6.256 mm/a)>生长季(0.524 mm/a)>非生长季(6.780 mm/a)；3个时段均存在2个主周期，第一主周期为14 a；ET值全年及生长季由东南向西北递减，90%以上区域为100～500 mm，非生长季ET小于200 mm。

b)全年、生长季和非生长季NDVI均显著增加，增速为生长季>全年>非生长季；不同时段年均NDVI由东南向西北递减，与ET空间分布基本一致。降水量在全年和生长季显著增加，非生长季增加不明显，温度和风速在3个时段内均无明显上升趋势。

c)生长季时段对ET影响有较大贡献的是风速、气温和NDVI，降水对其影响最小；非生长季中，NDVI、降水和气温对ET贡献较大，风速权重比例即贡献率相对较低；全年中黄土高原NDVI和降水对ET的贡献较大，其次为风速，气温对黄土高原ET影响并不明显。