南京地区蒸散发降尺度研究<br/>——基于增强型时空自适应反射融合模型

南京地区蒸散发降尺度研究
——基于增强型时空自适应反射融合模型

2022-08-31尉毓姣曹鑫宇王文科龚建师余慧琳

生态学报 2022年15期

尉毓姣,朱琳,*,曹鑫宇,王文科,龚建师,余慧琳,孟丹

1 首都师范大学资源环境与旅游学院,北京 100048 2 城市环境过程与数字模拟国家重点实验室培育基地,北京 100048 3 首都师范大学水资源安全北京实验室,北京 100048 4 长安大学水利与环境学院,西安 710054 5 中国地质调查局南京地质调查中心,南京 210016 6 自然资源部流域生态地质过程重点实验室,南京 210016

蒸散发是水文循环中自降水到达地面后由液态或固态转化为水汽返回大气的过程,包括水面蒸发、土壤蒸发和植被蒸腾,是陆地水分循环和能量交换的重要组成部分,能够调节全球能量和水量平衡。精确估算蒸散发以及获取高时空分辨率的蒸散发数据对区域乃至全球气候变化、水资源评价和陆地碳水循环具有重要意义[1—2]。

针对大范围、快速变化的地表信息监测,单一传感器无法满足需求,许多学者针对这一问题,提出了切实可行的遥感数据时空融合技术[3]。2006年,Gao等人[4]首次提出STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)融合算法,开创了遥感数据融合的先河。Yi等人[5]使用STARFM模型对MODIS和ASTER蒸散发结果进行数据融合,用于获取农田尺度上连续的日蒸散发,发现在非均质农业区融合精度低于均质区域。Yang等人[6]使用STARFM模型融合基于Landsat和MODIS数据估算的蒸散发,发现精度与仅使用Landsat估算的蒸散发相似或略低。STARFM算法进行遥感数据的融合可以较好的捕捉植被的物候变化,适用于地表覆盖同质的像元,而对于复杂、异质性的地表覆盖,预测效果不佳。Zhu等人[7]提出增强型时空自适应反射融合模型(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model,ESTARFM),该模型考虑了像元之间的光谱和空间相似性,保留了空间细节,适合复杂地物区域的数据预测。程筱茜等人[8]利用ESTARFM模型模拟了内蒙古红碱淖1987—2018年间缺失的Landsat影像数据,用于反映水体的变化情况,该模型模拟精度达到0.95,弥补了Landsat数据在时间分辨率上的不足。Zhou等人[9]利用ESTARFM模型对归一化植被指数(Normalized Difference Vegetation Index,NDVI)和地表温度(Land Surface Temperature,LST)分别进行融合,用于提高干旱监测的精度,发现预测结果与实际数据的相关系数分别为0.9和0.74。Heimhuber等人[10]分别利用STARFM和ESTARFM模型对Landsat影像和MODIS影像进行了8 d间隔、30 m分辨率的数据预测,发现ESTARFM结果优于STARFM,在漫滩地区融合精度最佳。邬明权等人[3]将STARFM、ESTARFM以及统计回归等模型进行对比,发现ESTARFM模型的融合精度最高。

目前,对于蒸散发的研究以中国北部地区居多,且研究内容以蒸散发估算、时空特征分析为主,而对于南部地区蒸散发降尺度的研究较少。因此,本文选择南京市为案例区,结合Landsat- 8遥感影像数据和气象数据,采用基于能量平衡原理的SEBS(Surface Energy Balance System)模型估算南京市的地表日蒸散量。在此基础上,选取南京市部分区域(400像元×400像元,面积144 km2),采用ESTARFM模型将高空间分辨率(30 m)的Landsat蒸散发估算结果与高时间分辨率(8 d)的MOD16A2产品数据进行时空融合,并评价模型的融合精度。研究结果可为南京市水资源的合理规划提供基础数据。

1 研究区概况

南京市(31°14′—32°37′N,118°22′—119°14′E)位于中国东部、长江下游中部地区,是长三角辐射带动中西部地区发展的国家重要门户城市,总面积6587 km2。区内四季分明,雨水充沛,多年平均温度为16.6 ℃(2000—2019年),年平均降雨量为1150.8 mm,属亚热带季风湿润气候。长江穿城而过,沿江岸线总长近200 km,境内共有大小河道120条,水域面积占总面积的11%。南京市是中国重要的农业和商品粮基地之一,植物资源丰富、种类繁多。全市森林覆盖率27%,建成区绿化覆盖率45%,人均公共绿地面积13.7 m2。研究区示意图如图1所示。

图1 研究区示意图Fig.1 Location map of study area

2 数据与研究方法

2.1 数据来源与预处理

本文采用的数据包括遥感影像数据和专题数据,其中专题数据包括：气象数据、土地利用类型数据、数字高程模型以及MOD16A2蒸散发产品数据。

2.1.1遥感影像数据

由于时空数据融合需要尽量选择日期相同或相近的两种数据源,且融合需要三个时期的数据,本次经过数据筛选后,发现Landsat- 8遥感影像在2017年7月21日、10月9日和12月12日数据质量较好(云量小于2%),且能够与MOD16A2蒸散发产品保持日期一致。故选择覆盖南京市的2017年7月21日、10月9日和12月12日三期Landsat- 8遥感影像(地理空间数据云：http://www.gscloud.cn/),空间分辨率为30 m,时间分辨率为16 d。在ENVI 5.3中对影像数据进行辐射定标、大气校正、图像裁剪等预处理,通过波段运算工具计算归一化植被指数NDVI、地表比辐射率和反照率参数,采用大气校正法计算地表温度[11]。

2.1.2专题数据

(1)气象数据

气象数据来源于中国气象数据网(http://data.cma.cn/),选取南京市及周边城市共计11个站点的日风速、气温、气压、相对湿度等气象要素,作为SEBS模型的输入数据。在ArcGIS 10.5中采用反距离加权插值方法将站点数据进行空间插值,得到30 m分辨率的栅格数据。

(2)土地利用类型数据

土地利用类型数据为2018年1 km空间分辨率的栅格数据,来源于中科院资源环境科学与数据中心网站(http://www.resdc.cn/Default.aspx)。原数据包括耕地、林地、草地、水域、居民地和未利用土地6个一级类型和25个二级类型。在ArcGIS 10.5中对二级类型进行合并处理,并将合并结果重分类为耕地、林地、水域、草地和其他五类(表1),通过裁剪得到南京市土地利用类型图(图2),其中其他类别主要包括建设用地和未利用土地。

表1 南京市土地利用类型分类表

图2 研究区土地利用类型图Fig.2 Land use type map of the study area

(3)数字高程模型(DEM)

数字高程模型(Digital Elevation Model,DEM)为30 m空间分辨率的SRTM产品,下载于Google Earth Engine平台(https://earthengine.google.com/),并通过裁剪处理得到南京市DEM图。

(4)MOD16A2蒸散发产品数据

MOD16A2蒸散发产品数据来源于美国国家航空航天局(NASA)网站(https://ladsweb.modaps.eosdis.nasa.gov),空间分辨率为500 m,时间分辨率为8 d。本次通过MRT(MODIS Reprojection Tool)工具对MOD16A2产品数据进行格式转换、投影转换等预处理,并重采样至30 m格网大小,便于后续与基于Landsat-遥感影像反演的蒸散发结果进行融合计算。此外,该产品数据对建设用地、水域、裸地等区域赋值无意义数值,为无效值(https://lpdaac.usgs.gov/products/mod16a2v006/)。为了减少对融合结果的影响,本次采用邻域均值代替无效值,并选取无效值较少的区域(400像元×400像元,面积144 km2)作为后续的融合计算区。

2.2 研究方法

2.2.1SEBS模型

2002年,荷兰学者苏中波[12]提出基于能量平衡原理的单层模型—SEBS模型,该模型具有较高的估算精度和更好的实用价值,目前广泛应用于干旱监测、水资源评价等方面[13—16]。该模型假设每个像元的显热通量都介于该像元的干限和湿限显热通量之间,一定程度上避免了气象数据空间插值不确定性带来的误差[17—18]。SEBS模型的地表能量平衡方程[19]为：

Rn=G0+H+λE

(1)

式中,Rn为地表净辐射通量(W/m2)；G0为土壤热通量(W/m2)；H为感热通量(W/m2)；λE为潜热通量(W/m2),其中E为蒸散量(mm),λ为水的汽化潜热(Wm-2mm-1),与地表温度Ts(K)之间的关系式为[20]：

λ=[2.501-0.002361×(Ts-273.15)]×106

(2)

(1)地表净辐射通量Rn

地表净辐射(辐射平衡)指地表向上和向下的总辐射之差,是地面所获取的净辐射能量[21]。计算公式为：

(3)

ε=0.004×Pv+0.986

(4)

式中,Pv为植被覆盖度,公式计算如下：

(5)

式中,NDVI为归一化植被指数；NDVISoil为裸土或无植被覆盖区的NDVI值；NDVIVeg为植被完全覆盖区的NDVI值。本次取累计概率为5%和95%的NDVI值作为NDVISoil和NDVIVeg,对植被覆盖度进行计算[23]。

(2)土壤热通量G0

土壤热通量是单位面积土壤在单位时间内的热交换量[24],是由于传导导致的储存到土壤和植被中的热量的比率[20],经验统计公式为：

(6)

(3)感热通量H

感热通量也叫做显热通量,是指由于温度变化而引起的大气与下垫面之间发生的湍流形式的热交换[25],是由地表温度和参考高度的气象参数获取的,与温度差值成正比[26],计算公式为：

H=ρCp(Ts-Ta)/γa

(7)

式中,ρ=1.29(kg/m3)为空气密度；Cp=1004(Jkg-1K-1)为空气定压比热容；Ta为参考高度的温度(K)；γa为空气动力学阻抗(s/m)。空气动力学阻抗γa与风速、粗糙度和大气层结构有关,根据大气边界层相似理论,存在以下关系：

(8)

(9)

(10)

(4)蒸发比

基于能量平衡公式计算的蒸散发是卫星过境时刻的瞬时值,本文采用恒定蒸发比法将瞬时蒸散发扩展到日蒸散发[29]。该时间尺度扩展方法简单易用,适用于晴朗或云量保持恒定不变的白天[30]。计算公式为：

(11)

(12)

式中,ET24为日实际蒸散发量(mm)；EF为蒸发比；Rn24为日净辐射通量(W/m2)；G024为日土壤热通量(W/m2)；ρw为水的密度。

2.2.2ESTARFM时空融合模型

ESTARFM时空融合模型是对STARFM模型的改进,其原理是考虑同一时间同一区域内像元空间和光谱的相似性,依据预测日期前后至少两对高低分辨率影像,以及预测日期当天的一幅低分辨率影像,计算对应像元之间的权重和转换系数,从而模拟出预测日期的高时空分辨率数据[31]。本次利用待预测时期(2017年10月9日)前后两个时期(2017年7月21日、12月12日)基于Landsat-8遥感影像反演的蒸散发数据和MOD16A2产品数据以及待预测时期的MOD16A2产品数据,模拟出待预测时期30 m的蒸散发数据[32]。计算公式为[5]：

ET(xw/2,yw/2,tp)=TmETm(xw/2,yw/2,tp)+TnETn(xw/2,yw/2,tp)

(13)

(14)

式中,ET为最终预测时期的蒸散发；tp为预测时期；ETm为Tm时期预测的蒸散发；ETn为Tn时期预测的蒸散发；ETL为基于Landsat-8遥感影像反演的蒸散发；ETM为MODI6A2蒸散发；ETk为Tk时期预测的蒸散发；(xw/2,yw/2)为中心像元的位置；w为相似像元搜索窗口；(xi,yi)为第i个相似像元的位置；N′为相似像元的个数；Vi为转换系数；Wi为综合权重因子,计算公式如下：

(15)

Di=(1-Ri)di

(16)

(17)

式中,di为距离权重；Ri为光谱相似权重。

Tk为Tm、Tn时期的时间权重因子,计算公式为：

(18)

3 结果与分析

3.1 SEBS模型结果验证

以往研究利用通量站数据、蒸发皿观测数据或P-M公式计算的潜在蒸散发量与蒸散发模型估算的结果进行对比验证,以证明结果的合理性[33—35]。因未收集到通量站实测数据,本文采用南京市气象站点(31°55′N,118°54′E)在3个研究时期的蒸发皿折算数据以及MOD16A2产品数据进行SEBS模型30 m Landsat蒸散发估算结果的验证。由于蒸发皿实测数据为水面蒸发,本次参考相关文献确定该研究区水面蒸发的折算系数[36—37],将水面蒸发转化为实际蒸发后与SEBS模型结果进行对比。统计结果如表2所示,模型估算结果与蒸发皿折算数据的平均相对误差为0.14 mm/d,与MOD16A2产品数据的平均相对误差为0.22 mm/d。三者变化趋势相同,并与已有研究[38—39]的蒸散发验证结果的表现规律一致,可认为SEBS模型的估算结果合理。

表2 SEBS模型计算结果与蒸发皿折算数据、MOD16A2产品数据结果对比

3.2 SEBS模型估算的蒸散发时空特征

研究区3个不同时期的日蒸散量空间分布如图3所示,7月21日平均蒸散量为7 mm,最大值为10 mm；10月9日平均蒸散量为3.9 mm,最大值为5.6 mm；12月12日平均蒸散量为1.1 mm,最大值仅为2.8 mm。由此得出,蒸散量在季节上差异较大,呈现夏季>秋季>冬季的变化规律。7月21日为植被的生长季中期,气温高,降水充沛,日照时间长,植被生长比较旺盛,尤其是农田水分灌溉充足,植被覆盖度高,蒸腾量增加等都是导致该时间蒸散量较高的原因。10月9日为植被的生长季末期,气温降低,降水减少,气候已发生明显变化,并且农田水分灌溉减少,植被覆盖度降低,因而蒸散量有所降低。12月12日正值冬季,气温较低从而导致蒸散量很低。

图3 研究区日蒸散量空间分布图Fig.3 Spatial distribution of daily evapotranspiration in the study area

在空间上,南京市各区在7月21日的平均蒸散量差异较大,最大为六合区的7.7 mm,最小为秦淮区的5.4 mm。六合区耕地面积626km2,水域面积202km2,植被蒸腾量和水面蒸发量在夏季都比较大,所以该地区蒸散量最大。秦淮区属于城市中心区域,只有少部分自然河、人工河错落和极少的城市绿地分布,蒸散量最小。

为了进一步分析研究区不同土地利用类型的蒸散发差异,选取特征最为明显的7月21日的蒸散发结果,将土地利用与Landsat蒸散发数据进行叠加并分区统计,得到研究区不同土地利用类型的蒸散量。由表3可知,南京市不同土地利用类型的蒸散量大小为：水域>林地>耕地>草地>其他。水域是蒸散量最高的地物类型,因其水量充足,水面蒸发量最大；林地具有涵养水源的能力及较高的植被蒸腾,虽然面积占整个南京市面积的比例不大,但蒸散量较高；耕地面积最大,约占南京市总面积的51%,并且南部高淳区分布大量水田,夏季较高的气温、充足的灌溉条件等使得耕地蒸散量较高。

表3 研究区7月21日各土地利用类型平均蒸散量统计

3.3 ESTARFM模型结果及精度评价

本次输入ESTARFM模型的初始数据为2017年7月21日、12月12日的两组Landsat、MOD16A2数据以及预测时期10月9日的MOD16A2数据,以此融合预测10月9日的Landsat蒸散发数据,融合区域范围见图1。

图4为融合区10月9日基于ESTARFM模型的融合结果与基于Landsat-8遥感影像反演的蒸散发结果,图中亮白色代表蒸散发较高,为植被覆盖度较高的区域；黑灰色代表蒸散发较低,主要为其他类型的地物。经统计,融合结果最大值为5.8 mm,最小值为0 mm,平均值为4 mm；Landsat蒸散发结果最大值为5.5 mm,最小值为0 mm,平均值为3.6 mm。整体上看,融合结果数值偏高,但二者在空间分布上具有一定的相似性。为了进一步度量时空融合模型的预测结果精度,本次选择融合区域内均匀分布的10000个验证点,并计算二者的相关系数。图中(图5)散点分布在1：1线附近,相关系数为0.74,表明基于ESTARFM模型融合的高分辨率蒸散发结果可靠。

图4 融合结果与Landsat蒸散发对比Fig.4 Comparison of fusion result and Landsat evapotranspiration

图5 Landsat蒸散发与融合结果相关性分析 Fig.5 Correlation analysis between Landsat evapotranspiration and fusion result

4 讨论与结论

4.1 讨论

4.1.1SEBS模型的不确定性分析

本文利用Landsat- 8遥感影像数据进行地表参数的反演,采用SEBS模型估算地表日蒸散量。研究结果经过多种方法的验证,表明SEBS模型估算的蒸散发具有合理性。但是模型本身会受到输入气象要素与地表参数精度的影响,带来模型反演结果的不确定性。因此,保证模型输入参数的准确性对于提高模型的精度极其重要。同时,受地面通量站实测数据的限制,精度验证的结果受到了一定的制约,获取充足的蒸散发实测验证数据是需要进一步解决的问题。此外,基于地表能量平衡的蒸散发模型分为单层模型和双层模型,其中单层模型又称为大叶模型,把地表视为一张大叶与外界进行水分和能量的交换,在下垫面均匀的陆面反演蒸散发的精度较高。本文采用的SEBS模型是单层模型的典型代表,虽然应用较为广泛,但实际上研究区下垫面并非单一均匀组成,这在一定程度上也会对反演结果造成误差。而双层模型是单层模型的延伸,能够定量描述陆面非均匀性对地表通量的影响,分别获得土壤蒸发和植被蒸腾[40—41]。故在后续研究中将考虑下垫面异质性问题,采用特征空间的双层蒸散模型进行区域蒸散发反演方面的研究。

4.1.2ESTARFM模型的不确定性分析

本文采用ESTARFM模型进行蒸散发产品数据的时空融合,研究结果表明,该模型可以有效地进行数据时空降尺度。为了进一步讨论融合结果的质量,将融合区土地利用分别与ESTARFM模型融合结果和Landsat蒸散发估算结果进行叠加并分区统计,分析二者在土地利用类型上的蒸散发差异。融合区不同土地利用类型在10月9日的平均蒸散量统计结果如表4所示,由于融合区域范围小,没有涉及草地这一地物类型,故没有做相应的统计。由表可以看出,二者均呈现水域>林地>耕地>其他的特点。ESTARFM模型融合结果与基于Landsat-8遥感影像反演的蒸散发结果的土地利用蒸散量差异在0.2—0.4 mm之间,其中,差异最大的土地利用类型为其他类型和水域,这可能与MOD16A2产品本身在建设用地、水域、裸地等区域存在无效值有关。由此可见,输入ESTARFM模型的遥感数据质量对融合结果的质量会产生直接影响。

本次采用的Landsat- 8遥感影像为2017年7月21日、10月9日和12月12日三期,长时间间隔对融合结果会产生一定的影响。此外,融合结果会受到模型本身参数的影响,导致相似像元的选取会产生一定的误差,使模型存在一定的不确定性。本次只针对小部分区域进行了时空数据融合降尺度研究,在后续的工作中可以尝试将该模型应用到更大的尺度上,并且尽量选择同一季节、质量更好的产品数据进行融合,为获取长时间序列的蒸散发数据集提供数据保证,这也是未来进一步研究的方向。