基于MOD16的银川平原地表蒸散量时空特征及影响因素分析

2021-05-25王卓月孔金玲张在勇刘慧慧蒋镒竹钟炎伶张静雅

水文地质工程地质 2021年3期

王卓月，孔金玲，李英，张在勇，刘慧慧,5，蒋镒竹，钟炎伶，张静雅

（1.长安大学地球科学与资源学院，陕西西安 710054；2.长安大学地质工程与测绘学院，陕西西安710054；3.宁夏回族自治区水文环境地质勘察院，宁夏银川 750011；4.长安大学水利与环境学院，陕西西安 710054；5.宁夏回族自治区自然资源勘测调查院，宁夏银川 750002）

蒸散发是水文循环过程中极其重要的环节，同时也是自然界水资源管理中最基本、最不确定的要素。蒸散发将地表水分与地表能量平衡联系起来[1]，由同时发生的两个过程控制，即土壤蒸发和叶面蒸腾。潜在蒸散发是指一定气象条件下，充分供水地表（植被、裸土或开阔水域）的蒸散发，是表征区域干旱程度的重要指标之一。

由于地表的高度异质性，大多数传统方法通常只提供点尺度的估算，不能有效进行区域尺度的计算。相较传统的监测方法，遥感技术可以针对地表的异质性和传热过程的动态性进行区域尺度的蒸散量监测[2]，被认为是量化蒸散发时空变化与空间尺度扩展的重要工具。目前采用遥感技术进行区域尺度蒸散量反演的方法主要有经验统计公式、能量余项模型、特征空间法等[3]。基于区域尺度蒸散量计算方法，众多学者及团队研发了蒸散量遥感数据产品，如NASA基于陆面模型提供的全球陆面通量数据GLDAS[4]；英国布里斯托大学Miralles等[5]基于Priestley-Taylor方程发布的GLEAM产品；美国地质调查局（USGS）基于简单化的陆面能量平衡模型（SSEB），以Penman-Monteith公式为理论支撑从而获取实际蒸散量的SSEBop/ET产品[6]；美国航空航天局（NASA）基于MODIS卫星遥感数据，考虑多种阻力算法制作的全球陆地蒸散量产品MOD16等[7]。MOD16产品因具有高时空分辨率、获取方式简单、对蒸散量遥感监测有较好的反演精度等特点而广受国内外学者认可，Sebinasi Dzikiti等[8]以南非干旱区为研究区，比较了MOD16和PT-JPL两种蒸散量估算模型的精度，结果表明与PT-JPL模型相比，MOD16具有准确分析干旱环境下蒸散量的潜力；Hassan等[9]采用MOD16数据与SEBS模型反演的埃及月蒸散量进行比较，两者存在高度相关性；何慧娟等[10]探讨了陕西省地表蒸散量的时空变化特征，证实MOD16数据具有较高的精度，满足陕西省蒸散发研究要求。

由于降水量低、年平均气温高以及延长的干旱期等不利条件，银川平原地区长期存在水资源空间分布不均和对传统灌溉技术过于依赖等问题，因此需要对区域水资源进行精确地估算及管理[11]。以往蒸散发研究多针对场地尺度，然而受黄河水系分布、降水等因素影响，银川平原地表蒸散量空间差异性极大，场地尺度研究难以反映区域尺度上的时空变化规律。

本文以MOD16蒸散量数据、地面气象实测数据为基础，对2004—2019年银川平原地表蒸散量的时空变化特征进行分析、模拟，并预测了2024年地表蒸散量的发展趋势。研究结果可为掌握银川平原水循环规律、合理开发利用水资源以及保护生态环境提供科学依据。

1 研究区和数据来源

1.1 研究区概况

银川平原位于宁夏回族自治区中部黄河两岸，东依鄂尔多斯高原，西至贺兰山，南接黄土高原，北止石嘴山，南北长度约为165 km，东西长10～50 km（图1）。银川平原属于温带干旱区，年均日照时数约3 000 h，年平均降水量仅约200 mm，是我国典型的引黄灌溉区，黄河过境水量年均达300×108m3。与此同时，银川平原年蒸发量高达1 600 mm，强烈的蒸发浓缩作用加剧了当地盐渍化的发展，成为制约农业发展的一个重要因素[12]。

图1 研究区概况Fig.1 Overview of the study area

1.2 数据源及数据预处理

本文主要采用MOD16数据进行银川平原蒸散量时空特征分析，并结合气象实测数据、Landsat数据探讨影响因素。

MOD16数据是NASA基于Penman-Monteith模型[7]，综合考虑植被蒸腾、土壤表面蒸发等条件提供的全球蒸散量数据集。本文选用2004—2019年MOD16A3数据，空间分辨率500 m，时间分辨率为年（https://search.earthdata.nasa.gov）。利用MRT工具对该数据集进行投影转换、波段提取等预处理，并根据NASA提供的MOD16A3产品说明对水体、建筑用地等无植被覆盖区进行无值化处理，将其转换成WGS-84基准面、UTM投影的蒸散量、潜在蒸散量初始数据。

气象数据源自中国气象数据共享网（https://data.cma.cn/），筛选银川平原内9个站点数据，主要包括气温、降水、日照时数、相对湿度及蒸发量气象数据。统计分析2004—2019年实测蒸发皿数据并处理成年值数据，用于MOD16蒸散量在银川平原的适用性验证。

Landsat数据源自地理空间数据云网站http://www.gscloud.cn/，选取2004年、2014年7月份遥感影像，根据《土地利用现状分类》国家标准并结合银川平原实际土地利用情况，采用支持向量机及面向对象的分类方法，将银川平原土地划分为水田、旱田、林地、草地、水体、建设用地及荒漠化土地7类，并通过实测调查取点验证，分类精度满足研究需求。

2 研究方法

2.1 Theil-Sen Median趋势度

Theil-Sen Median趋势度方法主要用于蒸散量时序数据在像元尺度上的趋势分析[13-14]。该方法与常用的一元线性回归方法不同，其采用时序数据点的中值进行趋势分析，减小了异常值造成的误差。非参数Sen斜率度结果可用来评估研究时段内该时序数据单位时间的变化：

式中：Si、Sj——研究时段中第i、j年对应的数据。

趋势度结果β通过正负反映该时序数据增加、下降的变化趋势。

2.2 Mann-Kendall检验

Mann-Kendall检验是一种以非参数检验为基础，进行时间序列数据趋势分析的检验方法。世界气象组织（WMO）推荐使用这种方法来评估水文气象时间序列是否存在单调趋势[15]。该方法将非参数检验与非线性趋势检验相结合，使结果合理化，主要用于蒸散量时序数据趋势变化显著性分析及突变检验[16]。

Mann-Kandall方法通过检验统计量Z值评估蒸散量时序数据显著趋势的存在：

式中：S——Kendall检验统计量；

xi、xj——第i、j年对应的蒸散量；

n——时序数据的长度，本文为16。

在特定α显著性水平下，采用双边趋势检验，根据正态分布表中的值判断蒸散量趋势变化是否显著。若该时序数据趋势变化通过了α显著性水平。本文选取0.01，0.05，0.10的α水平进行显著性检验，Z1−2/α的值分别为2.576，1.960，1.645，根据3个断点将蒸散量变化趋势分为极显著、显著、弱显著以及不显著4种状态[17]。

Mann-Kandall可通过对蒸散量时序数据构造一组秩序列，用于总体趋势分析中的突变点检验。秩序列计算公式为：

此时设定蒸散量时序数据是独立、互不干扰的，UF检验统计量为：

通过计算，得到蒸散量UFi曲线。通过UFi曲线与α显著性水平下对应的正态分布表临界值Uα对比：若|UFi|＞Uα，意味着蒸散量时序数据存在显著的趋势变化。将逆序的蒸散量时序数据X(xn,···,x2,x1)重复以上计算步骤，并将计算结果乘以−1，得到统计量曲线UBi。通过两条统计量曲线在α置信区间内的交点个数及位置，可以判断该时序数据突变发生的次数及相应的时间。

2.3 CA-Markov模型

CA-Markov模型可用于蒸散量时序数据的模拟预测。元胞自动机（CA）设定元胞在时间、空间均处于独立、离散的状态，是一种能够模拟复杂系统空间演化过程的模型。该模型基本思想为：复杂系统内的任一元胞互不干扰，执行各自对应的局部转换规则。CA模型计算公式为：

式中：St、St+1——复杂系统内t、t+1时间的元胞状态；

N——设定的元胞邻域范围；

f——从t到t+1时刻元胞执行的转换规则。

CA-Markov模型在元胞自动机对空间格局精确模拟的基础上，结合Markov模型对时序数据预测的显著优势，实现了模拟结果在时间、空间上的合理化[18]。本文根据2004年、2009年、2014年、2019年4期年蒸散量数据统计特征进行分级，采用Markov模型计算相应时段的转移概率矩阵，以5 a为反演时段，对2024年银川平原的蒸散量进行模拟预测。

3 结果与分析

3.1 MOD16数据适用性分析

在采用MOD16数据对银川平原蒸散量时空特征分析之前，需对数据适用性进行检验。本文采用实测蒸发量与潜在蒸散量的相关分析进行适用性检验。由于MOD16数据是年尺度，因此，将气象站点蒸发皿日值数据处理成年值，解决了气象数据中异常值及不同口径大小对结果的影响，实现两种数据在时间尺度上的统一。提取该站点所在像元点的潜在蒸散量，在点尺度上进行各个气象站点蒸发量实测值与潜在蒸散量的相关分析。结果（图2）显示气象站实测蒸发量与潜在蒸散量在0.01水平上显著相关，相关系数为0.88，表明MOD16数据可用于银川平原地区蒸散量的相关研究。

3.2 蒸散量动态特征

2004—2019年银川平原蒸散量年际变化见图3。蒸散量整体上波动较大，增加趋势明显。年际蒸散量为270.3～396.1 mm，多年平均蒸散量为330.7 mm。研究时段内蒸散量距平相对变化率波动强烈，为−18.3%～19.8%。2004—2009年蒸散量低于研究时段内平均水平，2010—2019年（2013年除外）蒸散量均高于多年平均水平。其中2005年蒸散量最小，为270.3 mm，2016年蒸散量最大，为396.1 mm，差值达到了125.8 mm。

图2 MOD16潜在蒸散量与蒸发皿实测年蒸发量拟合关系Fig.2 Correlation between potential evapotranspiration of MOD16 and measured annual evapotranspiration of pan

图3 2004—2019年银川平原蒸散量年际变化Fig.3 Interannual variation of evapotranspiration in the Yinchuan Plain from 2004 to 2019

在分析时序数据时，除分析蒸散量总体变化趋势外，还需要判断并检验该时序数据在研究时段内突变发生的时间以及频率。在Matlab中采用Mann-Kendall法对2004—2019年银川平原年蒸散量序列的突变情况进行检验（图4）。结果表明：（1）近16年来，年平均蒸散量整体呈现波动增加的趋势；（2）2010—2019年UF曲线超过α=0.05显著性水平临界值（1.96），说明该时段年蒸散量增加趋势显著；（3）UF与UB曲线在置信区间（α=0.05）内有1个交点，位于2010年，说明银川平原年平均蒸散量在2010年左右发生了明显突变，突变后年蒸散量增加了45.2 mm。突变前2004—2009年时间段年蒸散量均低于多年平均水平，2010年突变后高于多年平均水平。

图4 银川平原年蒸散量Mann-Kendall统计量曲线Fig.4 Mann-Kendall statistical curve of annual evapotranspiration in the Yinchuan Plain

3.3 蒸散量空间变化特征

分析银川平原多年平均蒸散量、潜在蒸散量基本空间特征（图5）可知，2004—2019年银川平原多年平均蒸散量、潜在蒸散量空间分布具有明显差异。多年平均蒸散量变化范围为175.7～577.8 mm，整体呈现出南高北低、中间高两侧低的空间格局。从各市行政单元看，吴忠市蒸散量最高，银川市次之，石嘴山市蒸散量最低。多年平均潜在蒸散量为1 429.3～1 929.1 mm，整体呈现南低北高、中间低两侧高与蒸散量分布特征相反的差异特性。Bouchet等[19]认为受区域能量的影响，实际蒸散量减少，剩余能量将用于增加潜在蒸散量。这一区域蒸散互补相关理论可以解释蒸散量、潜在蒸散量空间分布的差异。多年平均潜在蒸散量以银川平原西侧荒漠带及石嘴山市最高，低值位于吴忠市、银川市内。蒸散干旱指数（EDI）[20]，即整数1与蒸散量/潜在蒸散量之差，可初步描述银川平原的干旱程度，该指数理论上量值范围为0～1，EDI值越大，说明该地区地表水分缺失越严重。图5（c）显示，银川平原地区整体干旱程度严重，EDI波动范围在0.66～0.90，石嘴山市最为干旱，银川市次之。吴忠市内黄河支系发达，整体地表水分较为充足，受水分胁迫程度较轻。

从银川平原2004—2019年蒸散量和潜在蒸散量趋势显著性结果分析可知（图6），银川平原蒸散量整体呈增加趋势，银川、吴忠、青铜峡城市周边区域蒸散量呈下降趋势。潜在蒸散量变化趋势与蒸散量趋势相反，石嘴山市零星区域潜在蒸散量呈上升趋势。整体上近16年来银川平原潜在蒸散量下降的同时，蒸散量在相对增加，符合区域蒸散互补相关理论，并从蒸散干旱指数方向分析，该地区干旱情况较之前有明显好转。

图5 银川平原2004—2019年年平均蒸散量、潜在蒸散量空间分布Fig.5 Spatial distribution of annual average evapotranspiration and potential evapotranspiration in the Yinchuan Plain from 2004 to 2019

图6 2004—2019年银川平原蒸散量、潜在蒸散量变化趋势空间分布Fig.6 Spatial distribution of evapotranspiration and potential evapotranspiration in the Yinchuan Plain from 2004 to 2019

银川平原绝大部分区域蒸散量呈增加趋势，占比达92.96%，极显著增加的范围最大，达75.49%；其次分别是不显著增加、显著增加及弱显著增加（表1）。蒸散量呈减少趋势的区域中以不显著减少为主，面积占比为4.41%；与蒸散量变化趋势相反，银川平原99.87%的区域潜在蒸散量呈下降趋势，极显著下降变化的面积占比最大，达53.68%，潜在蒸散量呈上升趋势的区域面积占比仅0.13%。

表1 银川平原蒸散量、潜在蒸散量变化趋势及显著性统计Table 1 Change trend and significance of evapotranspiration and potential evapotranspiration in the Yinchuan Plain

3.4 基于CA-Markov模型的蒸散量预测

基于CA-Markov模型，以2004年、2009年、2014年、2019年4期年蒸散量为基础数据，按照银川平原实际蒸散量统计特征，以200，300，400，500，600 mm为断点分为6个等级。通过Markov模型计算2004—2009年、2009—2014年、2014—2019年转移概率矩阵，选取合适的滤波器，进行预测次数的迭代，完成银川平原2009年、2014年、2019年蒸散量的预测。对模型模拟结果与实际蒸散量数据进行Kappa系数精度验证，精度检验结果显示2009年、2014年、2019年总体模拟精度分别为0.84，0.78，0.82，模拟精度较好，CA-Markov模型预测结果可信。在此基础上，以2019年蒸散量为初始数据，对银川平原2024年蒸散量空间分布进行模拟。

从2024年银川平原蒸散量模拟结果看（图7），低值主要分布在银川市、石嘴山市；与2019年相比，蒸散量整体上由低等级向高等级转移，其中银川平原北部以200～400 mm低蒸散量等级向400～500 mm较高蒸散量等级转移为主，南部主要由较高蒸散量等级向500～600 mm高蒸散量等级转移。

图7 2024年银川平原蒸散量 CA-Markov预测图Fig.7 CA-Markov prediction of evapotranspiration in the Yinchuan Plain in 2024

银川平原蒸散量等级面积占比统计结果见表2，100～300 mm蒸散量较低值区面积逐渐减少，300～400 mm蒸散量中等区面积基本持平，400～700 mm 蒸散量高值区面积持续增加。2024年CA-Markov模型模拟结果表明蒸散量整体呈增加的趋势。

表2 银川平原蒸散量等级面积占比Table 2 Proportion of the evapotranspiration grade area in the Yinchuan Plain

3.5 蒸散量时空变化影响因素分析

蒸散量的时空变化受自然因素与人类活动的共同影响，气候变化影响区域水热传输过程，土地利用格局影响地表水分的供给能力。

2004—2019年年均蒸散量与各气象因子年均值进行相关分析（表3），在时间尺度上探究银川平原蒸散量与各气象因子的关系。研究结果表明：蒸散量与气温、降水、日照时数呈正相关关系，相关系数分别为0.910，0.905，0.768；与相对湿度呈负相关关系，相关系数为−0.746，相关系数均通过了P=0.01的显著性检验。相关系数表明银川平原蒸散量对气温、降水具有更高的敏感性。银川平原日照充足，年均气温呈现逐渐增加趋势，这为蒸散发过程提供了良好的热力环境，而相对湿度增加限制了该地区的蒸散需求，由于蒸散量对气温具有高敏感性，银川平原蒸散量呈增加的趋势。银川平原年降水量较少，分布不均，除南部及黄河沿岸地区外，降水是蒸散发过程主要的水分来源，受降水控制，蒸散量呈现与降水分布一致的南高北低的空间格局。

为进一步掌握土地利用结构对蒸散量空间分布格局的影响程度，利用2004年、2014年银川平原土地利用遥感解译图分析不同土地利用类型的蒸散量特征。

图8 2004年、2014年银川平原土地利用类型Fig.8 Land use types of the Yinchuan Plain in 2004 and 2014

从2004年、2014年银川平原土地利用类型空间分布图（图8）来看，研究区土地类型以水田、旱田为主，呈现出南部水田、北部旱田的空间格局，这与多年平均蒸散量南高北低的空间分布特征相同。耕地的蒸散量是银川平原区域蒸散量的主要组成部分[21]，其中水田受地表水分胁迫程度较小，蒸散量高于旱田。荒漠化土地主要为裸土地、沙地等，地表裸露，土壤水分含量较少，蒸散量相对较低。草地多分布在荒漠化土地与其他土地类型过渡地带，降水是该类型地表水分的主要来源，水分缺乏，蒸散作用受限制。采用空间分析法统计各土地类型的平均年蒸散量（图9）。由于MOD16数据在水体及城镇用地的蒸散量无数据，为避免统计误差，在对各土地类型进行统计时不考虑这两类。两期土地利用类型的年蒸散量均呈现出水田＞旱田＞林地＞草地＞荒漠的规律，且2014年各土地类型的蒸散量均高于2004年，这一现象与银川平原近16年蒸散量整体变化趋势相同。