毛乌素沙地蒸散发时空分布及影响因素分析

杨亮彦，黎雅楠*，范鸿建，王雅婷

（1.陕西省土地工程建设集团有限责任公司，西安 710075；2.陕西地建土地工程技术研究院有限责任公司，西安 710021；3.山东省城乡规划设计研究院有限公司，济南 250013）

蒸散发（evapotranspiration，ET）主要包括水面和土壤蒸发、植被蒸腾[1]，在地表水循环和能量平衡中具有重要的作用。全球大约60%的地表降水通过蒸散发的方式返回大气，在干旱区高达90%[2]。同时蒸散发直接影响土壤湿度、地表温度[3]、地下水深度和植被生长状况[4]，是气候变化[5]、水文过程[6]及农作物估产[7]等的重要参数。随着生态环境退化、水资源短缺等问题的出现，明晰ET 的时空变化和影响因素具有重要的研究意义。

由于地表异质性，传统的地表蒸散发量获取方法以气象站提供的水面蒸发量为主，但该数据不能有效地反映区域尺度的地表蒸散发量，无法满足区域水资源合理开发利用的需要[8]。遥感技术的发展为地表蒸散发量的获取提供了新的方法，利用遥感数据可估算出净辐射通量、土壤热通量和显热通量，结合能量平衡公式得到潜热通量。国内外学者在此基础上相继提出了SEBAL（surface energy balance algorithm for land）、SEBS（surface energy balance system）、SEBI（surface energy balance index）等[9]ET 模型，以此获取地表实际蒸散发量。基于区域ET反演模型，众多学者及团队研发了地表蒸散发遥感数据产品，Miralles等[10]基于Priestley-Taylor 方程研发GLEAM 产品（0.25°，Daily）；NASA 基于MODIS 卫星遥感数据，利用彭曼（Penman Monteith，P-M）公式及多种阻力算法制作全球陆地蒸散发量产品MOD16（500 m，8 d）[11]；美国地 质调查 局（United States Geological Survey，USGS）利用陆面能量平衡模型获取实际蒸散量的SSEBop产品（1 km，month）[12]；Chen 等[13]基 于SEBS 模型获得了ET-EB产品；NASA 基于陆面模型和数据同化技术生成全球陆面通量数据GLDAS（0.25°，month）[14]。其 中，MOD16 产品因时空分辨率高、监测精度高、数据容易获取等优势，被广泛应用于典型地表ET时空分布研究。张特等[15]利用P-M 公式在澴河流域验证了MOD16 数据的精度，发现均相对误差（mean relative error，MRE）为7.5%，可满足区域ET 时空变化研究需求。Ma 等[16]基于MODIS 数据探究了2000—2015 年中国黄土高原ET 变化，结果表明，ET 在2000—2010 年有明显增加，但2011—2015年没有明显变化趋势。郭晓彤等[17]基于MOD16数据采用复相关关系探究了淮河流域ET 时空变化的驱动力，发现在淮河流域，人类活动是ET 时空变化的主要驱动力。

毛乌素沙地是中国西北典型的农牧交错区，生态环境十分脆弱。由于该地区降水稀少、气候干燥、太阳辐射强、蒸发强烈等自然因素，毛乌素沙地长期存在水资源短缺问题，土地荒漠化进程迅速，对西部和京津冀地区的生态环境造成严重危害，因此毛乌素沙地一直是中国干旱的热点研究区域。但前人的研究多集中在毛乌素沙地植被覆盖变化、荒漠化、土地利用变化等方面，在蒸散发方面以ET 反演模型估算局部地区日蒸散发量为主，鲜有针对毛乌素沙地降水量和实际ET时空变化的宏观研究。本研究基于MOD16 ET 数据集，利用趋势分析法明晰毛乌素沙地2001—2019年ET 时空变化特征，结合气象数据探究ET 与水热条件的驱动力关系，并分析土地利用∕覆被（land use∕land cover，LULC）和植被指数对ET时空变化的影响，为落实毛乌素沙地水资源管理及开发利用制度和生态环境保护政策提供参考。

1 材料与方法

1.1 研究区概况

毛乌素沙地地处中国西北部、黄土高原腹地，横亘榆林市北部、鄂尔多斯市南部及盐池县和灵武 市（37.45° —39.37°N，107.67° —110.67°E）（图1）。毛乌素沙地处于干旱与半干旱过渡区，是我国典型的农牧交错带，其西北部以沙地和草地为主，东南部以耕地和草地为主。研究区以温带大陆性气候为主，干旱少雨且分布不均，降水多集中在6—8月，年降水量250～440 mm，由西向东南方向递增。毛乌素沙地属于西北典型的生态脆弱区，是我国沙漠化最严重的地区之一，也是我国防风固沙重要的生态功能区之一[18-19]。

图1 毛乌素沙地区位Fig.1 Distribution of land use types in Mu Us Sandy Land

1.2 数据来源及处理

本研究采用的数据主要包括ET、气象、LULC和植被指数，其中ET 数据为MOD16 产品，由NASA 基于P-M 模型计算的数据集，空间分辨率为500 m，时间分辨率为8 d。选用2000—2019 年MOD16A3GF 数据集（https:∕∕search.earthdata.nasa.gov∕），由MOD16 产品全年统计合成的L4 级产品，时间分辨率为年。本研究利用MODIS Reprojection Tool（MRT）工具对HDF 类型数据进行格式转换、投影转换、裁剪、镶嵌等预处理。归一化植被指数（normalized difference vegetation index，NDVI）为MOD13数据（https:∕∕ladsweb.modaps.eosdis.nasa.gov），其空间分辨率为500 m，时间分辨率为16 d，每月可获取2 期影像，取二者平均值作为当月植被指数，并利用最大值合成法获取每年NDVI数据。

气象数据来源于中国气象数据网（https:∕∕data.cma.cn∕），研究区内共有8个国家气象站点，分别为横山站（编号53740）、靖边站（编号53735）、定边站（编号53725）、盐池站（编号53723）、神木站（编号53651）、榆林站（编号53646）、鄂托克旗站（编号53529）和乌审旗站（编号53644）（图1）。气象数据包含地表温度、气温、降水量等8种日值数据，本研究通过统计分析获取每个站点的年平均气温和降水量数据，并计算毛乌素沙地7个气象站点的年均气温和降水量。LULC 栅格数据为30 m 空间分辨率 的GlobeLand30（http:∕∕www.globallandcover.com∕），由中国研制，该数据采用WGS-84 坐标系，包括耕地、林地、草地、灌木地、建设用地、未利用地等10 个一级LULC。目前已发布2000、2010 和2020版。

1.3 评价指标

1.3.1 ET 变化趋势 趋势分析可模拟每个像元多年ET 的变化趋势，分析ET 在研究时段的时空变化规律。计算整个研究区ET 长时间序列的像元回归斜率K，若斜率为正则表示ET 随时间变化呈上升趋势，为负则表示ET 随时间呈下降趋势，且斜率的绝对值越大，说明ET 变化趋势越明显。其计算公式如下。

式中，K为ET 变化趋势斜率；n为监测时间段的年数；ETi表示第i年的ET。

1.3.2 相对变化率 相对变化率是相对值对时间的导数，主要用于探究研究时间段内因素的波动变化状况。

式中，r为相对变化率为研究时段ET 多年平均值。

1.4 数据分析

用分析软件EXCEL 和SPSS 进行数据统计和相关性分析。相关系数r取值范围-1～1，r为正表示参数为正相关关系，r为负表示参数呈负相关关系，绝对值越大相关程度越高。

2 结果与分析

2.1 毛乌素沙地ET时空分布

2.1.1 毛乌素沙地ET 空间变化特征 图2 为毛乌素沙地2000—2019年ET逐年空间变化趋势，可以看出，毛乌素沙地ET分布具有较强的空间差异性，整体呈西北低、东南高的空间规律，并且不同行政区域内ET分层明显，陕西省明显高于内蒙古自治区，这一特征可能与LULC和地表植被覆盖度有关。毛乌素沙地2000—2019年ET变化明显，除研究区西部的灵武市农业区ET较高外，20年内研究区ET在不断上升，其中在2000—2009年，ET在300 mm·a-1以下，2010年后高于300 mm·a-1的区域明显增多，主要分布在陕西省境内的定边县、靖边县、神木市、榆阳区和横山区。

2.1.2 毛乌素沙地ET 时间变化特征 为进一步分析不同时期ET的分布规律，统计了研究区2000、2005、2010、2015 和2019 年ET 分布区 间及占 比（表1），可以看出，2000和2005年，研究区ET主要集中在250 mm·a-1以下，分别占研究区的99.0%和89.5%，比例略微下降；到了2010 逐渐增加至250～300 和300～350 mm·a-1，分别占研究区的55.2%和23.4%；2015 年，研究区ET 出现回落，主要集中在200～300 mm·a-1之间，占研究区的74.6%；2019 年毛乌素沙地ET主要区间恢复到250～350 mm·a-1，其中300～350 mm·a-1首次成为ET 分布最多的区间。综上所述，毛乌素沙地年ET有上升趋势，且存在地域性差异，毛乌素沙地ET变化受行政区划生态环境政策和自然因素共同影响。

表1 毛乌素沙地ET分布区间及占比Table 1 Distribution range and proportion of ET in the Mu Us Sandy Land

从图3 可以看出，2000—2019 年研究区年均ET 在179.7～321.4 mm·a-1，多年平 均ET 为258.8 mm·a-1。研究时段内，高于和低于平均水平的年份分别有9 和11 个，其中高于平均水平的年份集中分布在2010—2019 年，低于平均水平的年份分布在2000—2009 年，距平相对变化率最大的年份为2000 和2016 年，ET 分别为179.7 和321.4 mm·a-1。以上数据表明，2000—2019 年毛乌素沙地ET 在波动中呈明显增加的趋势，且2010年之后波动范围变小。

图3 2000—2019年毛乌素沙地ET年际变化及距平相对变化率Fig.3 Interannual change and relative change rate of ET in the Mu Us Sandy Land from 2000 to 2019

2.2 毛乌素沙地ET时空变化趋势

根据毛乌素沙地ET 年际变化将研究时段分为2 个阶段：2000—2009 和2010—2019 年，并将ET 的趋势变化率分为7 个等级（图4，表2）。可以看出，毛乌素沙地各阶段ET 存在一定的差异性。在2000—2009年，毛乌素沙地ET平均变化速率为4.60 mm·a-1，呈快速上升趋势，ET增长速率集中在2～10 mm·a-1，占研究区的96.86%，其 中K＞5.00 mm·a-1的区域主要分布在研究区的周边，下降的区域仅占0.04%，主要位于灵武市的农业区；2010—2018年，毛乌素沙地ET的平均变化速率为5.71 mm·a-1，ET 继续保持增加趋势的同时，速率有所加快，与2000—2009年相似，ET的增长速率集中在2～10 mm·a-1，但K＞5 mm·a-1区域成为毛乌素沙地最主要区域，主要分布在研究区的中部和东部。整体分析，2000—2019 年毛乌素沙地ET 以上升为主，平均变化速率为6.87 mm·a-1，其中2～5 mm·a-1的区域占研究区的14.89%，分布在西北部；5～10 mm·a-1的区域占研究区的79.75%，主要分布在中部、南部和东部区域；K＞10 mm·a-1的区域占研究区的5.25%，散落在研究区的东部和南部；基本不变和下降的区域仅占0.11%。

表2 2000—2019年毛乌素沙地ET的趋势变化率不同等级的比例Table 2 Proportion of different grades of ET trend change rate in the Mu Us Sandy Land from 2000 to 2019

图4 2000—2019年毛乌素沙地ET变化速率空间分布Fig.4 Spatial distribution of ET change rate in the Mu Us Sandy Land from 2000 to 2019

2.3 ET影响因素分析

2.3.1 气候变化对毛乌素沙地ET 时空变化的影响 气候变化是影响区域水热条件空间分布的重要因素，ET 的变化与区域水热变化规律密切相关[20]。在区域研究中，年均气温表征区域热力条件，降水量表征水分条件，研究气温和降水量的年际变化有助于探究区域ET 的影响因素。图5 为2000—2019 年毛乌素沙地ET 与气温、降水量年际变化。研究区年均降水为337.8 mm，降水年际波动较大，上升趋势不明显，其中乌审旗站点2001—2010 年年降水量为383.1 mm，2011—2019年年降水量为430.1 mm，具有明显上升趋势；毛乌素沙地气温年际变化相对较小，年均气温在9.3 ℃左右。研究区ET 与降水量的年际变化趋势基本一致，呈波动上升趋势，除2019 年外，二者的相关系数为0.74；与气温变化趋势无明显相关性。以上数据表明，毛乌素沙地ET 受水热条件的影响，其中降水量为主要驱动因素。

图5 2000—2019年毛乌素沙地ET与气温、降水量年际变化Fig.5 Interannual changes of ET，temperature and precipitation in the Mu Us Sandy Land from 2000 to 2019

2.3.2 LULC 对毛乌素沙地ET 时空变化的影响毛乌素沙地LULC 以草地和耕地为主，其次为裸地、灌木地、林地、建设用地和水体。以2010 年为例（图6），草地、耕地、裸地、灌木地、建设用地、林地和水体及湿地占比依次为73.05%、14.63%、9.59%、1.73%、0.50%、0.26%和0.26%，草地遍布整个研究区，耕地分布在东部和南部，未利用地主要分布在研究区中部，呈条带状分布，其他LULC 散落分布在研究区内。由于MOD16 数据缺乏水体ET数据，本研究不再分析水体ET变化。

图6 毛乌素沙地土地利用∕覆被分布Fig.6 Distribution of LULC in the Mu Us Sandy Land

图7 所示为2000—2019 年毛乌素沙地不同LULC 面积及ET 变化。2000—2019 年间，建设用地和耕地面积增加，草地、裸地面积减少，未利用地面积变化不大。6 种LULC 的ET 均呈显著增加的趋势，与2000年相比，2019年耕地、林地、草地、灌木地、建设用地和未利用地ET 分别增加70.33%、81.43%、64.52%、57.50%、43.51% 和58.28%。不同LULC的平均ET从大到小依次为耕地＞林地＞灌木地＞建设用地＞草地＞裸地，耕地和建设用地面积的增加，裸地、草地面积的减少显著增加了毛乌素沙地ET，与研究区ET 增加的趋势相符。综上所述，LULC 是影响毛乌素沙地ET 变化的重要因素。但同时发现2019 年各LULC 的ET普遍高于2000年，说明其他因素对毛乌素沙地ET变化也具有重要影响。

图7 2000—2019年毛乌素沙地不同土地利用∕覆被面积及ET变化Fig.7 Changes of different LULC areas and ET in the Mu Us Sandy Land from 2000 to 2019

2.3.3 NDVI 对蒸散发时空变化的影响 为了更好地探究NDVI 对ET 的影响，本研究将毛乌素沙地多年平均NDVI 和ET 分为6 级（图8）。毛乌素沙地多年平均NDVI 为0.33，8.56%的区域NDVI≤0.20，主要分布在西北部和腹部；仅有1.25%的区域NDVI＞0.60，分布在灵武市农业区和靖边、定边境 内；0.20≤NDVI＜0.60 的区域 占研究区的90.19%，遍布整个研究区（表3）。毛乌素沙地NDVI整体呈西北到东南逐渐升高的趋势，与研究区ET空间分布较为一致，表明ET与NDVI具有较好的相关性。为了定量化分析二者的关系，图9为NDVI 与ET 年际变化曲线，二者变化曲线具有高度的一致性。图10 为ET 与NDVI 多年均值的全区域散点图（红色代表点密度大），二者决定系数为0.628 8，进一步证明了ET 与NDVI 的相关性。整体来说，ET 的变化情况与植被生长状态具有高度的相关性，植被的生长趋势是影响ET的重要因素。

图8 2000—2019年毛乌素沙地多年平均NDVI和ET空间分布Fig.8 Spatial distribution of multi-year average NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

图9 2000—2019年毛乌素沙地NDVI与ET年际变化Fig.9 Interannual change of NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

图10 2019年ET与NDVI多年均值的全区域热点图Fig.10 Region-wide heat map of the multi-year average of ET and NDVI in 2019

表3 2000—2019年毛乌素沙地多年平均NDVI和ET不同区间占比Table 3 Proportion of different intervals of the multi-year average NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

3 讨论

ET 作为干旱区水循环的重要途径，明晰其时空变化规律对毛乌素沙地水资源合理开发及可持续发展具有重要的作用。毛乌素沙地ET 在时空变化上存在一定的规律性，在空间分布上具有明显的分异性，呈西北向东南逐渐增加的趋势，其中陕西行政区内的ET明显高于内蒙古自治区的，与刘静等[18]研究结果一致。ET与降水量年际变化趋势较为一致，与气温无明显相关性，表明在干旱半干旱区，降水量是ET的主要驱动因素，这与毛乌素沙地缺水的自然特征相符。在不考虑水体的情况下，耕地和林地是毛乌素沙地ET 最大的2 种LULC，与王卓月等[8]的研究结果一致。多年平均ET与NDVI空间分布一致，年际变化曲线相似，表明毛乌素沙地植被覆盖及长势是ET的主要因素，与谷佳贺在黄河流域的研究结果相符[1]。

2000—2019 年毛乌素沙地ET 增加明显，其增加速率略高于降水量的变化速率，考虑其原因是耕地、林地面积的增加导致ET 的增加，过度的植树造林存在恶化当地生态环境的可能性，因此如何平衡植树造林工程对研究区的影响是未来需要考虑的问题。本研究发现，在不同LULC 中，灌木地和草地的ET明显低于林地和耕地的，因此增加灌木地和草地面积可有效减少水分流失，调整林木类型是改善植树造林政策合理性的重要选择。

本文针对2000—2019 年毛乌素沙地时空特征及其与水热条件、LULC 和NDVI 的关系进行了深入研究，但研究仍存在一定的局限性：MOD16数据空间分辨率为500 m，且水体部分的ET 缺失，较粗的空间分辨率和缺失值的问题导致部分研究信息无法获取，研究缺乏完整性；由于研究区内气象站点数量较少，且气象站周边缺乏ET 数据，因此在时空立方体和新兴热点分析[21]等研究方法上受到一定的局限性；在探究相关驱动机制时，不同行政区域的土地政策不同，间接影响研究区ET的变化，本文仅从宏观角度对ET进行分析，并未区分行政区域解读相关土地利用政策，这将是日后工作重点研究的方向。