格尔木河流域山前平原区蒸散量的分布特征

2019-10-14朱晓倩金晓媚张绪财

水文地质工程地质 2019年5期

朱晓倩，金晓媚，张绪财，张京

(中国地质大学(北京)水资源与环境学院，北京 100083)

蒸散(Evapotranspiration，ET)是土壤蒸发和植物蒸腾的总称。蒸散量是自然水循环系统中十分重要的组成部分，是干旱内陆地区地表水和地下水的主要排泄方式，干旱区地表降雨量的90%以上是通过蒸散的形式返回到大气[1]。格尔木河是柴达木盆地第二大河，不仅为格尔木地区社会经济发展提供水资源保障，同时也是山前平原区地下水主要补给源和涵养植被生态的主要水源。河流自昆仑山流出后，经山前戈壁砾石区时大量下渗对地下水进行补给，后穿过细土平原区至洪积扇前缘，地下径流受岩性影响溢出形成泉集河，最终流经盐沼平原汇入终端达布逊湖排泄。由于格尔木河冲洪积扇下部地下水埋深浅，则会受强烈蒸发作用，使地下水位变化受到影响[2]，加之近年来对地下水的开采，径流汇入的终端达布逊湖面积逐渐缩小。因此对格尔木河流域长时间序列的蒸散量时空分布及其影响因素的分析就显得尤为必要。虽然前人对格尔木河流域蒸散发有过研究，但缺乏长时间序列蒸散量的时间与空间变化，无论时间尺度还是空间尺度均难准确反映格尔木河流域的年际与年内变化趋势。

就蒸散量估算的研究方法而言，传统的蒸散量估算方法多为基于Penman公式[3]及Penman-Monteith公式[4]估算潜在蒸发及非饱和下垫面蒸发。后发展的传统实测方法如波文比法、涡度相关法、蒸渗仪法等。波文比法[5-6]简单且精度较高，但其适用范围却因其要求下垫面均匀且无平流影响而受到很大的限制;蒸渗仪法精度最高可达0.01～0.02 mm，但该法所需费用高，无法进行大面积布点测量，还需专业技术人员进行监控。涡度相关法是直接测量乱流通量的方法，较其它实测方法更为精密和可靠，但也有必须设立在平坦、植被单一的下垫面上的缺点，且其维护成本较高[7-8]。而在研究区以往的水资源评价中也不可避免地涉及到潜水蒸发问题，多位学者[9-10]基于潜水蒸发强度与潜水埋深的关系利用阿维里扬诺夫经验公式对格尔木河流域蒸发量进行计算，而受包气带土质、岩性及气候等的影响，计算结果误差较大，同时计算所得结果仅是某时段内的蒸发量。总之，传统估算蒸散量的方法受下垫面、气象因素、植被等因素的影响，其结果存在的误差难以估计，且传统方法存在难以获取大尺度时空范围内参数的问题，这也使得长时间序列下大尺度范围的蒸散量时空分布的研究无法准确地进行。

而随着遥感技术的不断推广，基于遥感手段推出了许多地表蒸散发的估算模型，扩展了传统方法在时空尺度上的应用范围，同时减小了传统方法以点代面引发的误差，使得蒸散量的估算结果更为精准。Su Z[11-12]提出地表能量平衡系统(Surface Energy Balance System, SEBS)，其目的在于将卫星对地观测的可见光、近红外和热红外波段资料与大气模式输出数据结合，从而更连贯地估算大气湍流通量，该模型目前应用较为广泛。之后多位学者[13-18]将该模型应用于不同的空间及时间尺度，并进行了不同程度的验证，都得到了较好的蒸散量估算结果，证明了SEBS模型计算蒸散量的合理可行性。同时，利用全球陆面同化系统(GLDAS)生成的气象数据较气象站观测数据具有更好的时空连续性，由GLDAS气象数据估算得出的蒸散量结果准确性也相对较高。

因此本文以格尔木河流域山前平原区为研究区，应用连续序列的MODIS数据及GLDAS气象数据基于SEBS模型估算研究区2001-2016年间16年的蒸散量，并就蒸散量时空分布特征及其影响因素进行分析。通过准确估算研究区长时间序列的蒸散发时空分布，客观评价格尔木河流域可持续利用量，为研究区水资源合理开发利用及生态环境保护提供参考。

1 研究区概况

研究区位于格尔木河流域的山前平原区，北临达布逊湖，南靠昆仑山，地势南高北低，平均海拔为2 780 m。格尔木流域山前平原区具有典型的高原大陆性气候特征，干旱，多风少雨，蒸发强烈。多年平均降水量为42.7 mm，多年平均蒸发量约为2 600 mm，年均气温4.7 ℃，相对湿度32%。研究区地形由西向东倾斜且南高北低，地貌复杂，地貌类型由南向北分别为山前洪积平原、冲积平原、冲湖积平原以及湖积平原，地下水的赋存及运动依赖于巨厚的第四系松散沉积物，从而形成了研究区孔隙地下水系统，由北向南主要可分为盐沼地区、细土平原区及山前戈壁砾石区。研究区工业以盐化工为主，植被类型主要为枸杞、藜麦、芦苇、红柳等。

格尔木河为青海境内河流，发源于昆仑山北麓，经格尔木市，最终流入东达布逊湖，全长468 km，流域面积1.786×104km2。河水补给来源以高山融雪为主，次为河谷潜水，平均流量21.3 m3/s，年径流量6.41×108m3。

图1 研究区地理位置示意图Fig.1 Geographical location map of the study area

2 研究方法及数据

2.1 研究方法

(1)SEBS模型

SEBS(Surface Energy Balance System)模型主要通过地表能量平衡方程估算地表净辐射通量和土壤热通量，之后计算出感热通量，由地表净辐射通量、土壤热通量和感热通量得出潜热通量及蒸发比，最终估算得出日蒸散量。

任一时刻，地表的净辐射通量与该区域感热通量、潜热通量及土壤热通量处于平衡，由此得任一时刻的表面能量平衡方程为：

Rn=G0+H+λE

(1)

式中：Rn——净辐射通量/(W·m-2)；

G0——土壤热通量/(W·m-2)；

H——感热通量/(W·m-2)；

λE——潜热通量/(W·m-2)；

λ——水的汽化潜热/(J·kg-1)；

E——实际水分蒸散总量/(mm·d-1)。

潜热通量指沿垂直方向，单位时间、单位面积上地表和植被以及大气之间因为水的相变而交换的热量，而在地表和大气的能量交换中，潜热通量所占的比例即为蒸发比。其计算公式为：

(2)

式中：Λ——蒸发比；

Λr——相对蒸发比；

H——实际蒸散量与最大蒸散量之比。

日蒸散量计算公式为：

(3)

式中：Edaily——地表实际日蒸散量/(mm·d-1)；

λ——汽化潜热/(J·kg-1)；

ρw——水的密度/(1 000 kg·m-3)。

(2)一元线性回归分析

回归分析是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。在回归分析中，只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。通过一元线性回归分析可模拟研究区每个栅格的变化趋势，本文通过该方法分析研究区2001—2016年每年7月份蒸散量的变化趋势。其计算公式为：

(4)

式中：n——总年数；

MET·i——第i年的ET值，i为年序号(1、2、3……)。

计算得出的Slope可反映研究区各栅格点ET变化趋势。若Slope>0，则该栅格点ET呈增长趋势；若Slope<0，则该栅格点ET呈减少趋势。

(3)Mann-Kendall显著性检测法[19-21]

Mann-Kendall是一种已经广泛应用于水文数据分析非参数检验方法，该方法能很好地揭示时间序列的变化趋势，但是不需要样本数据服从一定的分布，同时也不受少数异常值的干扰。其检测过程如下：

对于时间序列X(x1、x2、x3……)：

(5)

其中，n为数据样本长度，sgn为符号函数，其定义如下：

(6)

当n<10时，直接使用统计量S进行双边趋势检验。

当n≥10时，统计量大致服从正态分布，使用检验统计量Z进行趋势检验，Z值计算公式为：

(7)

其中，

(8)

式中：m——序列中结(重复出现的数据组)的个数；

ti——结的宽度(第i组重复数据组中的重复数据个数)。

2.2 研究数据

(1)MODIS数据

本文选取2001—2016年中等分辨率的MODIS数据，共计1472景，基于SEBS模型估算地表蒸散量。其中包括空间分辨率分别为500 m、1 000 m，时间分辨率均为8 d的MOD09A1数据(地表反射率产品)和MOD11A2数据(地表比辐射率与地表温度产品)。

MODIS数据均由美国NASA网站获取(https://reverb.echo.nasa.gov/reverb)。为保证研究中数据的一致性，将数据统一重采样为500 m分辨率。

(2)GLDAS气象数据

由于研究区气象站较少，本文采用GLDAS(Global Land Data Assimilation System)数据中提取的比湿/(kg·kg-1)、风速/(m·s-1)、气温/K和气压/Pa四个波段参与SEBS模型蒸散发的计算。选取2001—2016年时间分辨率为1月、空间分辨率为0.25°的GLDAS气象数据，共计192景。GLDAS数据来源于https://disc.gsfc.nasa.gov/datasets。为保证研究中数据的一致性，数据需重采样为500 m分辨率。

(3)DEM数据

为定量地反映研究区地形的空间变化特征，本文选取了空间分辨率90 m的SRTM-DEM数据，并于ENVI软件中进行空间投影、重采样等预处理，使其与MODIS、GLDAS数据空间分辨率一致。数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http://www.gscloud.cn)。

(4)气象数据

为研究气象因素对蒸散量的影响，本文选取了研究区内格尔木气象站气温、降水、相对湿度、蒸发等气象数据。数据来源于 “中国气象数据网(http://data.cma.cn)”的中国地面气候资料日值数据集。

3 结果与分析

3.1 蒸散量时空分布特征

本文基于SEBS模型结合GLDAS气象数据得到格尔木河流域山前平原区2001-2016年的日蒸散量，乘以相应月份的天数得到月蒸散量，将各月蒸散量相加即为年蒸散量。通过研究区2016年蒸散量空间分布图(图2)可知：研究区北部为大片盐池，蒸散量最高，最高可达780 mm；研究区内蒸散量与河流分布特征具有较好的一致性，自南向北的河流流经区域蒸散量整体较大；中南部格尔木市东西两侧为研究区植被主要发育带，市西侧为格尔木西农场，以耕地为主，农作物类型主要为枸杞、藜麦、小麦等，市东侧除农田主要发育芦苇、红柳、沙柳、芨芨草等，蒸散量较高，约200～300 mm；中北部裸土地区和南部荒漠地区蒸散量较小，约100 mm。

图2 研究区2016年蒸散量空间分布图Fig.2 Annual ET distribution of the study area in 2016

对研究区2001—2016年蒸散量年际变化及多年平均蒸散量年内变化进行统计(图3)，结果表明：研究区2001—2004年蒸散量略有升高，到2005年有大幅度增长，增长到82.94 mm，总蒸散量约3.79×108m3，2005—2008年蒸散量基本保持不变，浮动于90 mm左右，2008—2016年蒸散量有所波动，但2001—2016年蒸散量总体呈上升趋势。其中2015年蒸散量最大，为185 mm，总蒸散量约8.45×108m3；2001年最小，为25 mm，总蒸散量约1.16×108m3。研究区多年平均年内夏季即6、7、8三个月蒸散量最大，约20 mm；冬季蒸散量小，11、12、1、2月蒸散量值在0.5～3 mm之间。通过与前人对格尔木河流域蒸散量估算结果[22～24]对比表明，本文计算所得蒸散量空间分布规律与前人计算结果吻合，这也进一步证明了SEBS模型估算蒸散量结果的准确可靠性；但由于研究区范围不同，蒸散量绝对值上有差异；同时，前人计算区域蒸散发时只采用格尔木一个气象站的气象值代表全区，计算结果有一定误差。本文采用GLDAS全区气象数据计算区域蒸散发，所得结果相对更为可靠。

图3 研究区2001—2016蒸散量年际变化规律图Fig.3 Annual variation of actual evapotranspiration in the study area from 2001 to 2016

3.2 基于线性回归和Mann-Kendall的蒸散量变化趋势分析

图4 研究区连续时间序列蒸散量空间变化趋势(a)及显著性变化分类(b)图Fig.4 Spatial variation trend(a) and classification of significant changes(b) of continuous time series evapotranspiration in the study area

为提高分类结果的精确性，将两种方法结合得到研究区16年间蒸散量空间变化趋势图(图4)，即当Slope<0且|Z|>1.96,序列呈显著下降趋势；当Slope<0且|Z|≤1.96,序列呈下降但不显著趋势；当Slope≥0且Slope≤10，序列呈基本稳定趋势；当Slope>10且|Z|>1.96,序列呈显著增长趋势；当Slope>10且|Z|≤1.96,序列呈增长但不显著趋势。最终将研究区连续时间序列蒸散量空间变化趋势分为五类：显著增长、轻微增长、基本稳定、轻微降低、显著降低。结果表明：线性回归分析法和Mann-Kendall显著性检测法得出的结果具有良好的一致性。研究区大部分地区蒸散量呈稳定不变状态，其面积占研究区总面积的68.58%；蒸散量显著增长的地区主要集中在研究区中南部植被发育带及北部盐池区，占研究区总面积的29.70%；仅极少地区有轻微或显著的减少趋势，蒸散量减少地区约占研究区总面积的1%；其余轻微增长地区占研究区总面积的0.72%。同时为分析格尔木河流域山前平原区各季节蒸散量变化趋势，将一年12个月按照春(3—5月)、夏(6—8月)、秋(9—11月)、冬(12—2月)划分为四季。分别估计四季蒸散量稳定时的值，统计了研究区2001-2016年四季占全年蒸散量的比重，并求其平均值，得到16年间四个季节在年蒸散量中所占比例值(表1)，结合实际，最终确定四季的Slopesta值：春季Slopesta=1.5，夏季Slopesta=6，秋季Slopesta=3，冬季Slopesta=0.5。

表1 研究区四季蒸散量占年蒸散量比例表Table 1 Proportion of four seasons evapotranspiration to annual evapotranspiration in the study area

按照确定的slope值得到研究区16年间四季蒸散量空间变化趋势图及显著性变化趋势图(图5、6)，结合显著性变化趋势分类的面积统计结果进行分析。

结果表明：格尔木河流域山前平原区夏季蒸散量增减幅度远高于其他三季，秋季和春季次之，冬季增减幅度最小。夏季和秋季均有大面积显著增长地区，主要集中在研究区北部盐池区及中南部植被发育带，研究区中北部及南部保持基本稳定趋势，小部分区域有轻微降低趋势。冬季研究区蒸散量除盐池部分有显著和轻微的增长外，其余大部分地区蒸散量都有轻微降低趋势。春季蒸散量也主要集中在北部盐池区，中南部植被发育带较夏季、秋季增幅减少，但也有少量显著增涨趋势。蒸散量变化趋势面积统计(图7)看，呈现显著增长趋势：夏季(26.46%)>秋季(25.49%)>春季(20.99%)>冬季(16.66%)；轻微增长趋势：冬季(4.87%)>春季(4.38%)>夏季(0.49%)>秋季(0.36%)；基本稳定趋势：秋季(73.66%)>夏季(71.54%)>春季(67.79%)>冬季(52.25%)；轻微降低趋势：冬季(22.49%)>春季(5.19%)>夏季(1.06%)>秋季(0.45%)；显著降低趋势：冬季(3.72%)>春季(1.64%)>夏季(0.45%)>秋季(0.04%)。总体上格尔木河流域山前平原区大部分地区蒸散量处于基本稳定状态，16年来北部盐池地区不断扩大导致盐池区域蒸散量显著增加。

图5 研究区蒸散量季节空间变化趋势图Fig.5 Seasonal spatial variation trend of evapotranspiration in the study area

图6 研究区蒸散量季节显著性变化趋势图Fig.6 Seasonal significance variation trend of evapotranspiration in the study area

图7 研究区全年及四季蒸散量变化趋势柱状图Fig.7 Histogram of annual and seasonal evapotranspiration variation trend in the study area

为进一步分析研究区蒸散量的时空分布，将研究区按山前戈壁砾质平原、绿洲平原、盐沼平原及盐湖区划分为四个地貌及水文地质单元，分别对不同地貌单元蒸散量进行统计分析(图8)。结果表明：研究区蒸散量空间分布与地貌类型空间格局具有较好的一致性。总体来看，研究区蒸散量盐湖区>盐沼平原>绿洲平原>山前戈壁砾质平原。四个地貌单元蒸散量从2001—2016年总体呈增长趋势，其中盐湖区增长速率最快。2010年后，由于洪水及耕地面积增加的影响，绿洲平原蒸散量增长速率加快，盐沼平原和山前戈壁砾质平原蒸散量16年间增长速率较慢。通过对不同地貌四季蒸散量的对比发现(图9)，任一地貌类型对应蒸散量随季节的变化与研究区16年间四季蒸散量的空间变化趋势(图5)一致，即除盐湖区外，其余三个地貌类型蒸散量均呈夏季>秋季>春季>冬季的变化规律。

图8 不同地貌单元对应蒸散量变化趋势Fig.8 Variation trend of evapotranspiration of different geomorphic units

图9 不同地貌单元对应蒸散量四季变化趋势Fig.9 Variation trend of evapotranspiration of different geomorphic units in four seasons

4 蒸散量影响因素分析

4.1 气象因素

本文利用研究区内格尔木气象站2001—2016年每月的气象数据与对应月的蒸散量进行相关性分析(图10)。结果表明：研究区蒸散量与气温相关性较好，确定系数达0.604。气温对研究区内盐湖区的影响较大。蒸散量与降水量相关性一般，确定系数为0.405。通过气温、降水与蒸散量的年内与年际变化趋势对比分析可以看出，蒸散量总体与气温、降水量变化趋势一致，呈正相关关系，即蒸散量随气温和降水量的增大而升高。

图10 气象因素对蒸散影响分析图Fig.10 Effect of meteorological factors on evapotranspiration

4.2 植被指数

植被蒸腾也是区域主要蒸发途径之一，研究区内植被主要发育带分布于绿洲平原中南部格尔木市东西两侧，因此植被对研究区蒸散量的影响主要表现在绿洲平原区。近年来，随着农田面积的不断扩大，绿洲平原区的蒸散量也在不断增大。为分析植被对研究区蒸发的影响，本文将连续序列的植被指数MODIS NDVI数据与计算得出的蒸散量进行相关性分析(图11)，结果表明：两者相关性较好，确定系数为0.612，呈正相关。

图11 植被指数与蒸散量相关性分析图Fig.11 Correlation analysis diagram of NDVI and evapotranspiration

图12 不同地貌蒸散量与NDVI变化对照图Fig.12 Comparison of variation of evapotranspiration and NDVI in different physiognomies

为进一步研究植被对研究区蒸散量的影响，本文利用图1中AB剖面同时提取2016年7月NDVI及蒸散量空间分布图中对应点的NDVI及蒸散量值，通过滑动平均法，对比分析不同地貌蒸散量与NDVI值的变化特征(图12)。结果表明：研究区由北向南AB剖面的不同地貌下，蒸散量和NDVI值总体变化趋势具有较好的一致性。除水体月蒸散量较高，最高可达90 mm外，NDVI值较低的裸土地区和荒漠地区蒸散量值也相对较低，在0～15 mm之间；NDVI值位于0.3～0.7的灌木地区，蒸散量也相应增大，波动于15～45 mm间；NDVI值位于0.2～0.4的农田地区，蒸散量值波动于20～40 mm间。因此，蒸散量与NDVI呈正相关，蒸散量值随NDVI的增大而增大。