基于NPP VIIRS数据和SEBS模型的河南省冬小麦蒸散量估算与时空特征*

李 颖,陈怀亮,梁 辰,苏 伟,贺 添

(1.中国气象局·河南省农业气象保障与应用技术重点开放实验室 郑州 450003; 2.河南省气象科学研究所 郑州 450003;3.哈尔滨市气象局 哈尔滨 150028; 4.郑州大学生态与环境学院 郑州 450001; 5.中国农业大学土地科学与技术学院北京 100083; 6.郑州大学地球科学与技术学院 郑州 450001)

农田蒸散量(ET)包括植株的蒸腾量和土壤的蒸发量[1],是农田水碳交换过程的关键变量[2]。准确估算农田ET,对于掌握农田水分的动态变化,监测预测旱情,制定科学的灌溉管理措施,提高农田水分利用效率,促进农业可持续发展等具有重要意义[3]。中国小麦(Triticum aestivum)种植面积占粮食作物种植总面积的22%左右,水分条件是制约小麦生长发育的关键因素之一[4],大尺度冬小麦ET 的准确估算对于保障我国粮食安全具有重要意义。利用大型蒸渗仪、涡动相关法、通量仪法等传统方法可获得地表点尺度ET,如胡程达等[5]利用大型蒸渗仪观测了充分灌水和自然降水两种状况下冬小麦的实际ET,但这类方法无法反映ET 的空间分布格局[6]。利用卫星遥感技术估算ET 具有快速、宏观、动态、空间连续监测等优点,近年来,基于遥感数据估算区域尺度地表ET 的研究取得了许多成果[7-8],其中一类代表性方法是利用遥感反演的地表温度和气象观测的大气温度进行感热通量计算,并进一步通过地表能量平衡方法计算潜热通量和ET[9-11]。

地表能量平衡系统(surface energy balance system,SEBS)是地表能量平衡方法中单层模型的一种,由荷兰Wageningen 大学的Su[12]在地表能量平衡算法(surface energy balance algorithm for land,SEBAL)基础上发展而来,近年来在国内外获得了广泛应用。Ma 等[13]使用SEBS 模型和ASTER 卫星数据估算了澳大利亚Coleamball 灌区的ET,并结合实测数据进行了验证。张亚丽等[14]提出一种基于面积比例加权聚合的空间尺度转换方法,分别利用HJ 卫星数据及中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer,MODIS)数据 ,结合SEBS 模型估算了伊洛河流域地表ET。Wu 等[15]利用MODIS 数据,通过优化空气热力学粗糙度参数化方案,提高了基于SEBS 模型的黑河中游ET 估算精度。张圆等[16]利用Landsat 8 影像,采用SEBS 模型估算了新疆呼图壁县的ET,并进行了时空格局分析。Mohammadian 等[17]和郑倩倩等[18]分别在伊朗中部半干旱地区和新疆精河流域开展了类似的研究。在上述研究中,中低空间分辨率的MODIS 数据和高空间分辨率的ASTER、HJ、Landsat 卫星数据均得到有效应用,其中MODIS 数据因高时间分辨率在大尺度ET 估算中具有应用优势。国家极地轨道卫星Suomi NPP (National Polar-orbiting Partnership)由美国国家航空航天局(NASA)等机构研发,其搭载的可见光红外成像辐射仪(visible infrared imaging radiometer suite,VIIRS)是对MODIS 传感器的继承与发展,具有优于400 m的星下点分辨率,每12 h 可提供一次全球影像[19],在陆地、海洋和大气参数的全球性观测等方面具有很大应用潜力。近年来,不少学者将NPP VIIRS 数据应用于气溶胶光学厚度、水体浊度、湖冰等的监测[20-23],苏城林等[24]将VIIRS 反演的气溶胶光学厚度(aerosol optical depth,AOD)对标MODIS AOD 产品进行检验,证明两者具有显著的相关性。将VIIRS数据用于估算农田ET 的研究尚鲜见报道[25]。

为探索NPP VIIRS 数据在ET 反演中的适用性,我们将VIIRS 数据反演的地表温度等参数和优化计算的归一化植被指数(normalized differential vegetation index,NDVI)数据代入SEBS 模型,通过运行SEBS 模型估算了河南省冬小麦农田ET,并进行了精度验证,进而开展了河南冬麦区ET 时空特征分析。本文旨在通过系统性的研究,一方面探究NPP VIIRS 卫星数据应用于ET 估算的适用性,为利用NPP VIIRS 数据和SEBS 模型估算ET 提供方法、技术上的支持; 另一方面通过研究ET 的时空分布格局,为河南省冬小麦干旱监测评估、农田灌溉管理以及农业可持续发展等提供科学依据。

1 研究区域与数据源

1.1 研究区域

本研究选取我国冬小麦种植第一大省河南省为研究区,位于110°21′E～116°39′E,31°23′N～36°22′N,属暖温带至亚热带、湿润至半湿润季风气候,冬季寒冷、雨雪少,夏季炎热、雨量丰沛。全省年平均气温普遍为12～16 ℃,1月份平均气温普遍为-3～3℃,7月份平均气温普遍为24～29 ℃,大体呈东高西低、南高北低分布。河南麦区耕作模式主要为冬小麦和夏玉米(Zea mays)轮作,冬小麦播种期一般为10月上中旬,收获期一般为次年5月下旬至6月上旬。研究区内冬小麦生育时期监测数据来源于河南省内农业气象观测站,研究区位置及农业气象站分布情况如图1 所示。河南省作为全国粮食生产大省,水资源严重不足,在河南麦区掌握农田水分的动态变化规律,进行科学的灌溉管理尤为重要。

1.2 气象数据与预处理

本研究所需气象数据包括气温、气压、相对湿度、辐射量、2 m 高度风速等,来源于2016年至2018年冬小麦关键生育时期河南省119 个气象观测站的实测数据,其中辐射量根据文献[26]由日照时数计算得到,2 m 高度风速由气象观测站实测10 m 风速转换计算得到。采用反距离插值法将站点尺度的气象观测数据拓展为面尺度气象数据,在ArcGIS 10.2软件中实现。

1.3 遥感数据与预处理

1.3.1 NPP VIIRS 数据

NPP VIIRS 传感器设置有22 个波段,包括9 个可见光至近红外波段,8 个短波、中波红外波段,4个热红外波段和1 个日夜波段(Day-Night Band,DNB),每天提供两次全球影像。VIIRS 传感器相较MODIS 传感器具有更高的灵敏度,且大多数波段分辨率高于MODIS 传感器,其星下点空间分辨率为371 m,扫描幅宽为3000 km。从美国宇航局戈达尔德航天中心(LAADS DAAC)网站下载覆盖2016年至2018年冬小麦关键生育时期河南省范围的VIIRS 产品级数据,利用HEG (HDF-EOS to GeoTIFF)软件从VIIRS L2级数据VNP09GA、VNP21A1D 和VNP43A1 中分别提取地表温度(land surface temperature,LST)、比辐射率(emissivity)、地表反照率(albedo)、太阳天顶角(sun zenith angle)等波段,将这些地表参数数据统一重采样为0.01°,重投影为WGS84 坐标系,使用ENVI5.6 软件进行辐射定标和研究区裁剪。

1.3.2 GF-1 WFV 数据

GF-1 卫星搭载了4 台WFV 传感器,可获取空间分辨率为16 m 的多光谱影像,幅宽为800 km[27-28]。从中国资源卫星应用中心网站下载2018年冬小麦生长季覆盖河南省境内的晴空影像,产品等级为Level 2 几何校正影像产品。使用ENVI5.6 软件对WFV 影像进行辐射定标、大气校正和以省辖市行政区划为处理单元的晴空影像拼接。

利用支持向量机(support vector machine,SVM)方法对河南省各省辖市WFV 影像进行监督分类,将土地利用类型划分为冬小麦田、水田、林草地、水体、建筑用地、其他用地共6 类。在各省辖市的分类影像上提取冬小麦田并拼接为河南省冬小麦种植区掩模,用于在预处理后的VIIRS 数据上提取研究区,该方法精度李颖等[29]在前期研究中已验证。

2 研究方法

2.1 基于SEBS 模型的ET 估算

地表能量平衡理论是SEBS 模型的理论基础,地表能量平衡方程为[30]:

式中:Rn表示地表净辐射通量,H表示感热通量,G0表示土壤热通量,λE表示潜热通量。

地表净辐射通量Rn指长、短波辐射经过相互抵消后地表所得到的净能量,是驱动水分交换和大气运动的主要能量,其计算公式为:

式中: α为地表反照率,来源于遥感影像;Rswd为 下行的太阳短波辐射; εa为大气比辐射率; ε为地表比辐射率,来源于遥感影像;Rlwd为下行的大气长波辐射;σ为斯蒂芬-波尔兹曼(Stefan-Boltzmann)常数,取值为2.68×10-8W·m-2·k-4;T0为地表温度;Rnl为地表长波净辐射;Rns为地表短波净辐射。

土壤热通量(G0)计算公式为[31]:

式中: Γc为植被覆盖区参数,Γs为裸土区参数,根据荷兰国际地球观测研究所开发的陆地及水体集成系统ILWIS 软件中的SEBS 模块,这两个参数分别取值0.050 和0.315;fc为植被覆盖率,由VIIRS NDVI影像计算得到。

感热通量(H)的计算根据大气边界层相似理论,由以下公式推导得到:

式中:u为风速;u*为摩擦速度; k为卡尔曼常数,本研究取0.4;z为参考高度;d0为零平面位移高度;zom为动量传输粗糙长度; Ψm为动力学稳定度修正函数;L为由感热通量(H)计算得到的奥布霍夫长度。

SEBS 模型的输入参数来源于气象数据、DEM数据和遥感数据,通过公式(2)至(4)可分别计算得到地表净辐射通量、土壤热通量和感热通量,进而可通过公式(1)得到潜热通量 λE,进而估算ET。

2.2 VIIRS NDVI 优化

NDVI 数据是SEBS 模型的一项重要的输入参数,用于反演叶面积指数(leaf area index,LAI),并进一步计算得到冠层内风速廓线消光系数[12]。由遥感影像计算得到的单日NDVI 易受天气情况的影响,本研究在利用VNP09GA 产品的I1 和I2 波段计算NDVI 影像时做出了优化。采用最大值合成法(maximum value composite,MVC)对单日NDVI 数据进行校正[32],选择应用较广泛的周最大值合成策略,即利用研究日当天及其前后3 d 的影像分别计算NDVI,取1 周内的最大值合成VIIRS NDVI 取代单日NDVI。

2.3 彭曼公式

本研究利用彭曼(Penman-Monteith,P-M)公式结合气象观测数据计算河南省冬小麦农田ET,用于基于VIIRS 数据的ET 估算结果的验证。彭曼公式如下[33]:

式中: ET0为参考作物的潜在蒸散量,Rn为地表净辐射,u2为2 m 高处的风速,Δ为饱和水气压曲线斜率,γ为干湿系数,T为2 m 高处日均气温,es为饱和水汽压,ea为实际水汽压。

由彭曼公式计算得到参考作物的潜在蒸散量ET0,实际ET 可由ET0乘以作物系数得到,本研究根据文献[34]得到冬小麦不同生育时期相应的作物系数,进而计算得到实际ET,记为P-M ET。

3 结果与讨论

3.1 研究区ET 估算结果

利用荷兰国际地球观测研究所开发的陆地及水体集成系统ILWIS 软件中的SEBS 模块,输入气象观测数据、90 m 空间分辨率的SRTM DEM 数据和NPP VIIRS 数据,对河南省冬小麦ET 进行估算,空间分辨率为0.01°,研究时段为2016年至2018年冬小麦关键生育时期(返青期至灌浆期)。以2018年为例,以8 d 为时间间隔且尽量选取晴空遥感影像,对3月8日至5月19日河南省冬小麦日蒸散量进行估算(图2)。图中无ET 数据的区域除山区等非冬小麦种植区外,还有少量冬小麦种植区内的像元因局部云影响等遥感数据质量问题赋为空值。采用Arc-GIS 10.2 软件中的自然断点分级法将每幅影像的ET值分为4 级,可直观显示冬小麦农田ET 在空间上的分布情况。其在空间格局上大体呈现中部和东南部较高,向西北部和西南部逐渐降低的趋势,在时间格局上呈现由返青期(3月8日)开始上升,至抽穗期(3月24日)达最大值,灌浆期(5月19日)开始下降的趋势。

3.2 模型验证与分析

3.2.1 结合彭曼公式进行精度验证

在返青期、拔节期和抽穗期分别采用P-M ET对本研究利用SEBS 模型和VIIRS 数据估算的VIIRS ET 进行验证。VIIRS ET 为面尺度数据,P-M ET 则为点尺度数据,为使两种数据表征相同位置的冬小麦农田ET,根据周边环境选取研究区内有代表性的9 个国家气象站,确保站点位于冬小麦种植区掩膜内。利用彭曼公式计算其站点ET0,根据文献[33]得到冬小麦返青期、拔节期和抽穗期的平均作物系数分别为0.8、0.9 和1.1,进而计算得到P-M ET; 再根据站点经纬度提取其对应像元的VIIRS ET,对比两种方法估算的日蒸散发如表1 所示。

表1 显示,2018年3月8日、3月24日和4月17日3 个时相的日VIIRS ET 与P-M ET 对比,平均绝对偏差分别为0.124 mm·d-1、0.063 mm·d-1和0.721 mm·d-1,总的平均绝对偏差为0.303 mm·d-1; 平均相对偏差分别为12.3%、2.7%及15.8%,总的平均相对偏差为10.3%。本研究使用SEBS 模型和VIIRS 数据估算的VIIRS ET 与采用联合国粮食及农业组织(FAO)推荐的彭曼公式计算得到的P-M ET 的相对偏差在作物不同生育时期在2%～16%的范围内波动。相对偏差的波动较大,主要因为难以在时间上和空间上准确确定作物系数。作物系数的测定过程复杂,适用于反映特定区域多年平均值。总体来看,VIIRS ET 与P-M ET 的总体偏差在10%左右,且两者的变化趋势具有较高的一致性。

3.2.2 结合大型土壤蒸渗仪实测蒸散量进行精度验证

研究区内郑州农业气象试验站(113°39′0′′,34°43′12′′)安装有大型称重式蒸渗仪,可获得冬小麦试验田ET 的实测数据,记为Real ET,用于本研究VIIRS ET 精度的验证。于2018年冬小麦返青期至灌浆期,以8 d 为间隔,在郑州农业气象试验站对比VIIRS ET 与Real ET 如表2 所示。

表2 显示,2018年3月8日至5月19日期间,VIIRS ET 对比郑州农业气象试验站实测ET 平均绝对偏差为0.20 mm·d-1,平均相对偏差为6.74%。以Real ET 为基准,进一步计算VIIRS ET 的均方根误差RMSE 为0.203 mm·d-1。证明本研究方法得到的VIIRS ET 具有较高精度。

表2 2018年冬小麦返青期(3月8日)至灌浆期(5月19日)可见光红外成像辐射仪数据估算蒸散(VIIRS ET)与大型称重式蒸渗仪实测蒸散(Real ET)对比Table 2 Comparison of between daily evapotranspiration estimated by visible infrared imaging radiometer suite data (VIIRS ET)and macro-weighting lysmeter measured ET (Real ET) from winter wheat regeneration stage (March 8th ) to filling stage(May 19th) in 2018

3.2.3 与MODIS ET 进行对比

利用MODIS 数据估算ET 的精度已得到广泛验证[35-37],可作为基准来评价其他遥感数据源估算ET的效果,如赵红等[38]将FY-3 VIRR 数据反演的日ET与同区域、同时段的EOS MODIS 数据反演的日ET进行对比分析,通过证实两种产品的一致性较好,指出利用国产FY-3 卫星资料估算ET 是可行的。因MODIS ET 产品MOD16 是空间分辨率1 km 的8 d、月、年尺度数据,与本研究0.01°空间分辨率的单日VIIRS ET 数据时空分辨率不一致,故对比分析时未直接使用MOD16 产品,而是采用本研究中VIIRS ET 的计算方法处理同一日期的MODIS 地表反射率产品MOD09GA 数据等,运行SEBS 模型计算MODIS ET,空间分辨率统一设为0.01°,以准确评估两种传感器数据反演ET 的一致性。本研究选择2016年至2018年冬小麦抽穗期代表性的1 天(4月17日),将相同技术方法反演的当日VIIRS ET 与MODIS ET 进行对比分析,每年在研究区随机且均匀选择40～50 个影像质量较好的像元作为对比样本点,2016年至2018年VIIRS ET 与MODIS ET 对比图如图3所示。

图3 中VIIRS 估算日蒸散量与MODIS 估算日蒸散量对比散点图显示,大部分样本点分布在1∶1线附近,对VIIRS ET 与MODIS ET 进行线性回归分析,2016年、2017年和2018年相关系数R分别为0.897、0.857 和0.896,决定系数R2分别为0.804、0.734 和0.802。以MODIS ET 为基准验证VIIRS ET的 RMSE 分别为 0.222 mm·d-1、 0.158 mm·d-1和0.211 mm·d-1。结果表明,VIIRS ET 与MODIS ET 具有较好的一致性。

图3 可见光红外成像辐射仪数据估算日蒸散量(VIIRS ET)与中分辨率成像光谱仪数据估算日蒸散量(MODIS ET)对比图Fig.3 Comparison between daily evapotranspiration estimated by visible infrared imaging radiometer suite data (VIIRS ET) and by MODIS data (MODIS ET)

3.3 河南省冬小麦农田蒸散量时空特征分析

本研究利用VIIRS 数据,结合SEBS 模型估算2016年至2018年河南省冬小麦关键生育时期农田ET,进而分析研究区日蒸散量的时空分布格局。首先计算2016年至2018年河南省冬小麦关键生育时期代表性日期的麦区平均日蒸散量(表3),分析返青期至灌浆期农田ET 的时间变化特征。河南省冬小麦农田ET 的时间变化趋势在2016年、2017年和2018年保持一致,在关键生育时期均以返青期的平均日ET 为最低,其后日ET 逐渐上升,抽穗期达到最大值,灌浆期开始下降。本研究得到的河南麦区平均日ET 时间变化特征与姚瑶等[39]研究得到的华北平原冬小麦不同生育阶段的ET 特征基本一致。结合郑珍等[1]的研究,分析冬小麦农田ET 时间变化特征与冬小麦生长发育规律的关系,在冬小麦生育初期,ET 主要来自于土壤蒸发,在快速生长期至生育中期,植株的蒸腾量变大,特别是抽穗期至灌浆期,冬小麦日均ET 达到最大,此时作物对水分亏缺敏感,应保证该生育阶段水分充分供应以满足作物生长需要,到了作物生育后期,冬小麦灌浆期至成熟期,随着籽粒开始成熟,叶子变黄,植株的蒸腾量随之降低,应减少该生育阶段的田间灌水量,有助提高作物产量。

表3 2016年至2018年河南省冬小麦关键生育时期麦区平均日蒸散量Table 3 Average daily evapotranspiration in key growth periods of winter wheat in Henan Province from 2016 to 2018mm·d-1

将河南省冬小麦农田ET 的时间特征与空间特征协同分析,冬小麦农田ET 在不同生育时期的空间分布均大体呈现中部和东南部较高,向西北部和西南部逐渐降低的趋势。分析河南省的地形地貌,西部、西南部和西北部多山地与丘陵,这些区域的冬小麦种植区多为雨养区,易发生春旱和秋旱,供给蒸发与散发的水分不足,各生育时期ET 均相对较低;中部和东南部则多平原,这些区域的冬小麦种植区多为灌溉区,水热条件较好,各生育时期ET 均相对较高。与张圆等[16]对新疆呼图壁县ET 时空格局受人类活动影响的分析一致,河南省冬小麦区ET 的时空变化特征与灌溉条件亦具有较好对应性,受农业管理活动影响明显。进一步对农业干旱进行评估时,可利用VIIRS ET 与潜在蒸散的关系计算作物缺水指数(crop water stress index,CWSI)[40],结合其他农业干旱监测指标,得到河南省干旱时空变化特征。

4 讨论与结论

本文通过MVC 方法优化计算VIIRS NDVI,将VIIRS NDVI 和VIIRS 数据反演的地表温度等参数与SRTM DEM、气象观测数据等输入SEBS 模型,估算了2016年至2018年河南省冬小麦返青期至灌浆期的VIIRS ET,验证了VIIRS ET 与P-M ET、Real ET 和MODIS ET 结果的一致性。前人研究中主要利用MODIS 数据和Landsat 卫星数据估算ET[41-42],本文利用NPP VIIRS 数据实现了河南省冬小麦ET的高精度估算,展示了这种新型遥感数据源用于农田ET 反演的适用性,为今后研究中利用NPP VIIRS数据和SEBS 模型估算ET 提供了方法、技术上的支持。张志高等[41]利用MODIS 蒸散产品MOD16研究了河南省2001-2019年蒸散量时空变化特征,得出河南省ET 多年均值西南部山区最高,东南部较高,西南部和中北部较低的空间分布特点,本文与张志高等结论里东南部ET 较高、西南部ET 较低的趋势相同,不同点在于本文得到中部ET 较高、西北部ET 较低的特点,原因在于本文是去除西南部和西北部等山地、丘陵后的冬小麦ET,且本文仅研究冬小麦关键生育时期的田间ET。河南中部地区灌溉条件良好,在冬小麦关键生育时期由于作物生长旺盛且灌溉供水充足,田间ET 较高,西北部则为冬小麦雨养区,田间ET 相对较低。张志高等[41]研究了多年平均ET,中部耕地冬季ET 较低,森林覆盖类型的年均ET 高于其他土地覆盖类型。通过上述讨论,在排除研究对象和研究时段不同的影响后,本文得到的河南省冬小麦ET 分布特点与前人研究具有良好的一致性。本文开展了针对农作物关键生育时期的田间ET 估算,可为指导河南省农业水资源管理、分配和高效利用等提供重要依据。本文的主要结论如下:

1)利用NPP VIIRS 数据,结合SEBS 模型估算了2016年至2018年河南省冬小麦关键生育时期ET。将本研究方法估算的VIIRS ET 与彭曼公式计算的PM ET 进行比较,两者的总体偏差在10%左右,且两者的变化趋势具有较高的一致性。将大型土壤蒸渗仪实测Real ET 用于精度验证,计算VIIRS ET 的均方根误差RMSE 为0.203 mm·d-1,验证结果表明,新型遥感数据NPP VIIRS 数据适用于ET 反演,本研究方法估算的VIIRS ET 具有较高精度。

2)将利用相同技术方法反演的VIIRS ET 与MODIS ET 进行线性回归分析,2016年至2018年R2分别为0.804、0.734 和0.802; 以MODIS ET 为基准验证VIIRS ET,3年的RMSE 分别为0.222 mm·d-1、0.158 mm·d-1和0.211 mm·d-1。 表明VIIRS ET 与MODIS ET 具有较好的一致性,VIIRS 数据可以作为MODIS 数据的有效补充与替代遥感数据源,在农田ET 估算中发挥积极作用。

3)河南冬小麦区ET 空间分布大体呈现中部和东南部较高,向西北部和西南部逐渐降低的趋势,ET的空间变化特征与灌溉条件具有较好对应性。研究区冬小麦关键生育时期农田ET 的时间变化特征以返青期的平均日ET 为最低,其后日ET 逐渐上升,抽穗期达到最大值,灌浆期开始下降。河南省冬小麦农田ET 的时空特征与冬小麦生长发育规律和田间管理密切相关,准确估算河南省冬小麦农田ET,可为设计灌溉管理制度提供科学依据,对保障区域粮食安全具有重要意义。