刘曼晴,胡德勇,*,于 琛,王莎莎

1 首都师范大学资源环境与旅游学院,北京 100048 2 资源环境与地理信息系统北京市重点实验室,北京 100048

湿地是水陆相互作用的独特自然综合体,在抵御洪水、调节气候、控制污染等方面具有不可替代的作用,被誉为“地球之肾”。水是维持湿地稳定和健康的关键生境因子[1]：湿地的水文过程通过改变湿地环境的生物、物理和化学特性来影响湿地的类型、结构和功能,进而控制着湿地的形成与演化[2]。由于受自然环境变化和人类活动的双重影响,湿地水循环过程发生了深刻变化,其生态水文过程成为了国内外学者研究的热点。

蒸发和散发(Evapotranspiration,简称“ET”)是湿地水分损失的主要途径[3],在很大程度上影响着湿地生态系统的水热平衡。因此,开展湿地蒸散量及其时空特征研究,对湿地生境保护、生物多样性保护、水资源管理、湿地退化修复、灾害防御等方面具有重要意义。

传统的ET测量主要采用水文、气象学原理和方法在站点或田间测定。该类方法主要的不足有2点：(1)大范围蒸散观测耗费大量人力、物力和财力；(2)生境因子的空间变异性使站点观测数据很难在区域尺度扩展。由于遥感数据具有快速、宏观、重复性观测等特点,所用光谱波段涵盖可见光、近红外和热红外等波段,使得遥感数据不仅能解决尺度扩展问题,而且能提供丰富且与蒸散密切相关的参数,从而使遥感技术在区域ET反演中得到广泛应用[4- 6]。目前,利用遥感反演ET的方法主要有指数法和余项法[7]。指数法将遥感反演的指数(植被指数、微波土壤湿度指数、能量平衡指数)与 Penman-Monteith公式结合计算蒸散量[8- 9]；余项法利用遥感技术获得的辐射和通量数据结合能量平衡方程估算蒸散量,主要分为单层、双层和多层模型[10]。单层模型将土壤和植被视为均匀的“大叶”模型[11],在众多单层模型中,“陆面能量平衡算法”(Surface Energy Balance Algorithm for Land,简称“SEBAL”)发展于美国植被密集的平原地区,以地表能量平衡方程为核心,结合遥感数据和气象资料便可估算区域尺度的蒸散量[12],适用于所有可见光、近红外、热红外传感器和不同的气候条件,特别是在我国东部地区已有很多成功的应用实例[10,13- 16]。由于研究区植被覆盖度较高且长势均匀,可视为单一层次,故选择SEBAL模型作为本文ET遥感反演模型。

辽河三角洲湿地在调节河流水位和保持水量平衡中起着重要的作用,在气候变化的大背景下,农业用水量逐年增加以及中下游水坝的修建,导致该区域生态用水不足。目前,已有多个学者在辽河三角洲开展了湿地ET研究,如王炳亮等[17]探究了1961—2010年潜在蒸散量的变化特征,于文颖等[18]估算了盘锦湿地2013—2015年生长季日蒸散量,并分析其空间和频率特征。然而,对于辽河三角洲湿地实际蒸散量特别是长时间序列的蒸散量变化及其时空格局研究甚少,加上以往研究中多采用MOD16产品数据[19]开展长时间序列区域ET研究,其数据空间分辨率对于小区域ET反演存在不足。

本文以辽河三角洲湿地为研究区,利用Landsat数据、气象数据等,基于SEBAL模型模拟湿地生长季蒸散量,并进一步分析蒸散量的时空特征及其影响因素,为增进对湿地水量平衡的理解和保护湿地生态平衡提供支持。

1 研究区和数据

1.1 研究区

辽河三角洲湿地位于辽宁省西南部辽河平原,是我国重要的河口湿地之一,本文研究区如图1所示。它包括盘锦、营口、锦州部分区域,由来自辽河、大凌河、小凌河的泥沙沉积而成,介于约121°10′—122°30′E,40°30′—41°30′N之间,面积约6000km2。该区域海拔小于7m,地势平坦。辽河三角洲属于温带半湿润半干旱季风气候,年平均温度约为8.6℃；年平均降雨量约为631mm,雨热同期；年蒸发量约为1390—1705mm。辽河三角洲是我国重要的芦苇、粮食和石油产地,在地区经济中起着支柱作用。1985年以来,由于社会经济快速发展,居民和工矿用地面积大幅增加,湿地景观破碎化程度加剧,人类活动正导致着辽河三角洲湿地生境特征和生态结构转变[20]。

1.2 数据

1.2.1遥感数据

本文所采用的遥感数据为Landsat系列卫星数据,包括Landsat- 5 TM和Landsat- 8 OLI/TIRS数据。在数据选择方面,由于生长季内研究区植被覆盖度较高,基于SEBAL模型反演的蒸散量更为准确,故本文对区域生长季蒸散量加以研究。但由于研究区地处海边,遥感影像大多有云层覆盖,加之Landsat卫星重访周期为16天,选择同日或同月数据较为困难。综上,本文选择1985—2017年7—9月(植被生长期)晴空少云资料共16景,对数据进行严格预处理,使之尽可能真实地反映地面属性。具体影像资料介绍见表1。

表1 遥感数据

1.2.2气象数据

研究区内有两个气象站：大洼和营口(图1)。本文所采用的气象数据来自中国气象科学数据共享服务网(http://cdc.cma.gov.cn/),包括太阳辐射、风速、气温、相对湿度和蒸散发数据,采用反距离权重法对输入模型的气象数据进行了空间插值[21]。

1.2.3土地利用数据

研究区土地利用/覆被专题信息提取于1985—2017年的Landsat影像,共分为20类(图2)。为方便后续统计,按照景观类型差异,将研究区分为四类：水体区(水库、盐池、河流、水池、养殖水体)、湿地植被区(水田沼泽、草甸、苇田)、非湿地植被区(旱地、灌丛、林地)、非植被区(滩涂、浅水沙洲、裸地、城镇居民点),滩地、工矿、道路、防护堤由于面积小,统计误差较大,因此不做分析。

2 方法和技术流程

2.1 研究方法

2.1.1SEBAL模型及其参数确定

陆面能量平衡方程如下：

Rn-G-H-λET=0

(1)

式(1)中,λET为潜热通量,W/m2；ET为蒸散量,mm；Rn为地表净辐射,W/m2；G为土壤热通量,W/m2；H为显热通量,W/m2。

图1 研究区示意图Fig.1 Schematic map of research area

图2 土地利用/覆被类型图Fig.2 Map of land use/cover types

(1)净辐射通量计算方法[22]

Rn=(1-α)RS↓+(RL↓-RL↑)-(1-ε)RL↓

(2)

式中,ε为地表发射率,无量纲；α为地面反照率,无量纲；RS↓为入射的短波辐射,W/m2；RL↓为入射的长波辐射,W/m2；RL↑为反射的长波辐射,W/m2。其具体参数获取详见参考文献[23]。

(2)土壤热通量估算方法[24]

(3)

式(3)中,I为归一化植被指数,无量纲；c11的取值与卫星过境时刻地方时有关,12:00前,c11=0.9,12:00—14:00,c11=1.0,14:00—16:00,c11=1.1。

(3)显热通量估算方法

(4)

式(4)中,ρair为空气密度,kg/m3；Cp为空气定压比热容(1004 J kg-1K-1)；rah为空气动力学阻抗,s/m；dT是地气温差,K。

(4)瞬时蒸散量估算方法

基于得到的Rn、H和G,由式(1)得卫星过境时刻的瞬时蒸散量(ET)：

(5)

式(5)中,λ是水的汽化潜热,计算公式为[23]：

λ=2.501-0.002361×(Ts-273.15)

(6)

2.1.2日蒸散量的估算方法

以上各通量均是基于遥感影像计算所得的瞬时通量,需转换为日蒸散量。假设24h内,蒸发比EF(潜热通量λET与可利用能量Rn-G之比)相对稳定[25]。

(7)

式(7)中,EF为瞬时蒸发比,无量纲；EFd为日蒸发比,无量纲；λETd为日潜热通量,W/m2；Rnd为日总净辐射量,W/m2,通过Allen等推荐的气象学方法计算[26]；Gd为日土壤热通量,W/m2,白天热量从地表向土壤传输,取正值,晚上则相反,且二者数量相当,因此Gd=0[26]。故日蒸散量ETd为：

(8)

2.1.3蒸散量反演结果的精度评价方法

本文拟选用气象站实测的小型蒸散发数据对蒸散量的反演结果进行验证。该数据采用小型蒸发皿(平底圆筒形、器口直径20cm、面积314cm2、深度10cm)测得,因蒸发皿的蒸发条件与自然水体不同,所以该数据不能准确代表自然水体的蒸散量。首先,设置不同的折算系数对蒸发皿数据进行校正,其转换方法如下[27- 28]：

E0=E1×0.98

(9)

E1=E2×k

(10)

式(9)—(10)中,E0为自然水体蒸散量；E1为大型蒸发器的蒸散量；E2为小型蒸发器的蒸散量,k是折算系数。1、2、3、11、12月由于冬季结冰故不进行观测,所以只提供生长季k值,如表2所示。

表2 辽宁E2对E1蒸散量的平均折算系数

其次,在蒸散反演结果中选择离气象站点最近的水体作为样点,并选择以水体样点为中心的3像元×3像元窗口内的9个像元蒸散量的平均值作为SEBAL模型的水面蒸发估测值；最后,进行蒸散量的SEBAL模型估测值与实测值之间的对比分析。

2.1.4气象因子的加权综合分析方法

本文拟通过对气象因子的加权综合分析来探究气象因子对蒸散量变化的影响。该方法主要利用蒸散量与几个关键气象因子(如平均气温、太阳辐射、风速、相对湿度)的相关性所占比例作为权重,以多个气象因子的归一化值作为统计参数,进行加权综合,得到气象因子归一化加权值,用于分析气象因子对蒸散量的影响。

蒸散量与气象因子的相关性用相关系数表达：

(11)

2.2 技术流程

本文技术流程如图3所示。首先,利用遥感数据反演地表参数(地表温度、地面反照率、植被指数、地表发射率),结合气象数据(气温)和DEM数据,计算研究区净辐射通量、土壤热通量、显热通量,基于辐射平衡估算瞬时蒸散量；再根据蒸发比恒定原理,估算研究区日蒸散量；最后,对所得估算结果进行精度验证,分析日蒸散量的时空格局及其影响因素。

3 结果与分析

3.1 蒸散量反演结果及精度评价

由以上方法反演得到研究区1985—2017年8期日蒸散量结果。于文颖等[18]利用涡度相关数据验证了蒸散量的反演结果,但在许多情况下,涡度相关数据不易获得。为此,本文选用郭生练和任芝花等的转换公式[27- 28],将从中国气象科学数据共享服务网获得的小型蒸散发数据转换为自然水体蒸散量,与蒸散反演结果中水面蒸发量进行对比分析。选取9组符合统计要求的数据进行比较,得到图4和表3。

图3 技术流程图Fig.3 Technique flow chart

图4 实测值与估测值相关性分析Fig.4 Analysis of correlation between measured and estimated values

由图4可知,估测值和实测值具有一定相关性,相关系数为0.61(P<0.1)。表3得出SEBAL模型估测值的相对误差为9.01%,与于文颖等[18]模拟精度接近,说明SEBAL模型反演蒸散量在本研究区可行。

表3 精度评价

3.2 蒸散量时空格局分析

3.2.11985—2017年日蒸散量均值变化特点

在分析长时间序列日蒸散量均值变化趋势前,本文首先对遥感数据的选择做验证,由于蒸散量遥感反演结果与遥感数据所选年份内生长季(7—9月)气象站实测数据平均值有较好的相关性,表明本文所选8期遥感数据基本能代表相应年份生长季蒸散量的平均水平,可以用于进行蒸散量的时空变化分析。

图5 蒸散量的时间变化特征Fig.5 Temporal change characteristics of ET

日蒸散量不同年份的变化过程及其相对变化率的年际波动情况如图5所示。1985—2017年,日蒸散量均值在3.09—3.70mm/d之间波动,整体呈双峰态势。研究时段,日蒸散量最高值出现在1989年(3.70mm/d)；最低值出现在2005年(3.09mm/d)；多年日均蒸散量为3.40mm/d(实线)。相对变化率波动情况也呈双峰态势,其中最低为-17.21%(2005年),最高达到8.31%(1989年)。

比较分析各年日蒸散量均值发现,日蒸散量超出多年日均蒸散量的年份出现在1985、1989、1995、2014年,其他年份均低于多年日蒸散量均值,以2005年最为突出,低于0.50mm。据此将研究时段分为：1985—2005年,日蒸散量呈显著降低趋势；2005—2017年,日蒸散量呈显著上升趋势。从整体上看,1985—2017年蒸散量微弱下降(虚线)。

3.2.2日蒸散量均值的空间格局

图6 平均蒸散量空间分布图Fig.6 Spatial distribution map of average ET

1985—2017年,日蒸散量空间分布如图6所示。总体而言,研究区日蒸散量具有较高的空间分异性,呈现水陆交界处最低、西部较低、中东部和南部高的趋势。由图6和图2可知,蒸散量的空间分异性表现为：研究区水陆交界处主要是潮间带光滩、裸地和浅水沙洲,植被覆盖度低,在不受海水影响的情况下蒸散量也相对较低；南部主要是水体,包括水库、河流等,蒸散过程以水面蒸发为主,蒸散量较高；东部主要是水稻田,水面蒸发和植被蒸腾都很重要；中部主要是芦苇沼泽,由于芦苇高度和覆盖度较高,以植被蒸腾为主,蒸散量较大。该结果表明在一定程度上蒸散与区域内土地利用/覆被类型有较大联系。

3.3 蒸散量影响因素分析

影响区域蒸散量时空格局的因素极有可能源于下垫面特性的土地利用/覆被和长时间的区域气象条件变动[10],它们驱使蒸散量形成独特的时空特征,下面分别对此进行分析。

3.3.1土地利用/覆被类型对蒸散量的影响

以2000年为例,统计该年不同土地利用/覆被类型的日蒸散量(图7)。总体来看,水体及其周围湿地的蒸散量较大,植被区蒸散量大于非植被区,蒸散主要与植被覆盖度及水分供给状况有关[29],蒸散量从大到小排序依次为：水体区>湿地植被区>非湿地植被区>非植被区(除水体外)。相关研究指出,水体蒸发是湿地蒸散的主体[30]；作为研究区面积最大的两种植被的水田沼泽和苇田,水田沼泽日蒸散量高于苇田的结论与于文颖等[18]结论相似。

图8为1985—2017年不同土地利用/覆被类型的蒸散量变化趋势,水体蒸散量最高,非植被区蒸散量在大多数年份最低,按蒸散量从大到小排序为：水体区>湿地植被区>非湿地植被区>非植被区。由于非植被区中的滩涂、浅水沙洲等多位于浅海水域附近,该区域受潮汐的影响较大,使得蒸散量大小不可避免受水体影响。多年土地利用/覆被类型下蒸散量平均值(虚线)变化趋势与图5所示的蒸散量年际变化趋势一致,均呈双峰态势。

图7 不同土地利用/覆被类型平均蒸散量Fig.7 Average ET of different land use/cover types1: 草甸；2: 居民点；3: 灌丛；4: 旱地；5: 河流；6: 林地；7: 裸地；8: 浅水沙洲；9: 水池；10: 水库；11: 水田沼泽；12: 滩地；13: 滩涂；14: 苇田；15: 盐池；16：养殖水体

为进一步探究土地利用/覆被类型变化对蒸散量的影响,本文分析了不同土地利用/覆被类型面积的变化情况及由于面积的改变导致日总蒸散量的变化情况(图9)。可以看出,水体、非湿地植被区、非植被区面积均呈增加趋势,增加量分别为12.75km2、103.89km2、25.36km2；湿地植被区面积减小,且变化量最大(167.53km2)。研究时段内,水体、湿地植被区、非湿地植被区、非植被区单位面积平均蒸散量相差不大,但由于较大的面积变化量也会给日总蒸散量带来较大影响,需先剔除不同土地/覆被类型单位面积平均蒸散量的影响。从图9中曲线可知,水体、非湿地植被区、非植被区日总蒸散量增加,增加量分别为4.34×104m3、31.25×104m3、7.52×104m3；湿地植被区日总蒸散量减小,且变化量最大(54.45×104m3)。不同土地利用/覆被类型日总蒸散量的变化特征与其面积变化特征一致。结果表明,研究区土地利用/覆被变化主要是湿地植被的转型,日总蒸散量随土地利用/覆被面积的变化而变化。

图8 不同土地利用/覆被类型蒸散量变化趋势Fig.8 Trends of ET in different land use/cover types

图9 不同土地利用/覆被面积与日总蒸散量变化情况Fig.9 Changes of land use/cover area and total daily ET

3.3.2气象因子对蒸散量的影响

蒸散主要受下垫面性质、太阳辐射、气温、风速、相对湿度的影响[31]。图10为气象因子与日蒸散量的相关关系,结果表明：太阳辐射、气温、风速、相对湿度与日蒸散量均呈正相关关系,相关系数分别为0.37、0.73、0.24、0.30,说明气温可能是影响辽河三角洲湿地蒸散量变化的主要因子。

日蒸散量的变化原因复杂,是多种因素综合作用的结果。图11为蒸散量与气象因子加权值的关系。图中曲线均是多项式拟合曲线,两者吻合度较好,相关系数为0.69,说明四个气象因子加权能解释蒸散量变化的原因。由图5和图8可知,2005年蒸散量为曲线拐点。而2005年平均气温、风速、太阳辐射和相对湿度分别为16.70℃、1.70m/s、17.65M J m-2和53%,除相对湿度外,其他三个气象因子值都是序列中的最低值。由此,我们可以推断：太阳辐射低,对土壤、植被和空气的增温作用较弱；气温较低,使得下垫面的水分运动较弱；而近地面风速较小,减弱了近地面空气流动过程伴随的水汽扩散和热能传递能力,这可能最终导致蒸散速率较低,植物蒸腾和地表蒸发量较少。

图10 蒸散量与气象因子的关系Fig.10 The relationship between ET and meteorological factors

图11 蒸散量与气象因子加权值的关系 Fig.11 Relationship between ET and weighted values of meteorological factors

4 结论

本文基于SEBAL模型反演辽河三角洲湿地1985—2017年生长季蒸散量,分析其时空格局,探讨不同土地利用/覆被类型、气象因子对蒸散量的影响,主要得到以下结论：

(1) 蒸散量估测值与实测值相比,相对误差为9.01%,相关系数为0.61。SEBAL模型能被用来估算辽河三角洲湿地蒸散量,估算结果可应用于时空变化研究。

(2) 时间变化特征表现为：1985—2017年平均蒸散量及相对变化率呈双峰态势,其中极小值出现在2005年,极大值出现在1989年和2014年,区域蒸散量呈微弱下降趋势；空间变化特征表现为：辽河三角洲湿地蒸散量具有显著的空间分异特征。

(3) 土地利用/覆被类型和气象因子是湿地蒸散量时空变化的主要影响因素。多年土地利用/覆被类型蒸散量平均值的变化趋势与蒸散量的年际波动趋势一致,均呈双峰态势；不同土地利用/覆被类型按蒸散量从大到小排序依次为：水体区>湿地植被区>非湿地植被区>非植被区(除水体外)；土地利用/覆被变化能直接影响区域日总蒸散量。在区域尺度上,选取与蒸散量变化关系显著的气象因子(太阳辐射、气温、风速、相对湿度),采用因子加权方法能在一定程度上解释蒸散量的时间变化特征。

本文选择气象站点的蒸散发实测数据来进行精度评价,这种验证可能存在一定的不确定性；湿地蒸散同时存在蒸发和蒸腾两种过程,两种方式对总蒸散的影响仍需进一步研究；同时,气象因子的加权方法解释蒸散量的时间变化特征仅被初步应用在辽河三角洲,需要在其他区域进行验证。