基于涡动相关通量观测的农田蒸散发产品精度验证

刘萌，彭中，黄凌霄，李召良，段四波，唐荣林

1.中国农业科学院农业资源与农业区划研究所/北方干旱半干旱耕地高效利用全国重点实验室,北京 100081;

2.中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室,北京 100101;

3.中国科学院大学,北京 100049

1 引言

地表蒸散发ET（Evapotranspiration）是指从地表传输到大气中的水分，包括土壤表面水分蒸发、植被水分蒸腾和冠层截留蒸发。水分蒸发的过程需要吸收能量来冷却地表，这部分能量被称为潜热通量LE（Latent heat flux ），是蒸散发的能量表达形式，二者可以相互转化（Wang 和Dickinson，2012）。蒸散发和潜热连接了水循环和能量循环，是地表水热循环中的关键过程参量，也是其中最难估算的地表特征参量。遥感技术因其大范围观测能力，为获取丰富的区域地表特征参量提供了可能，遥感反演蒸散发是获取区域、国家和全球尺度蒸散发时空分布的有效方法（Li 等，2009）。目前，国内外学者已经先后生产发布了不同时空分辨率的区域及全球蒸散发遥感反演产品，如MOD16产品（Mu 等，2011）、LSA−SAF 产品（Ghilain 等，2011）、SSEBop产品（Senay等，2013）、GLASS产品（Yao等，2014）、ETMonitor产品（Hu和Jia，2015）、BESS 产 品（Jiang 和Ryu，2016）、GLEAM 产 品（Martens 等，2017）、PML−V2 产品（Zhang 等，2019）、FLUXCOM 产品（Jung 等，2019）和CR−ET 产品（Ma 等，2019）等。相关研究已经详细介绍了现有的过程驱动的蒸散发遥感产品（李佳 等，2021）和数据驱动的蒸散发遥感产品（刘萌 等，2021）。由于不同产品之间存在反演方法机理差异、产品特征差异等，在使用这些产品前，需要对产品的精度进行验证和评估（熊育久 等，2021）。然而，由于近地表大气的复杂多变性、地表下垫面的异质性等特点，蒸散发产品的验证目前仍面临诸多困难和挑战。张圆等（2020）详细综述了蒸散发真实性检验的发展现状和研究进展，将蒸散发真实性检验方法分为像元尺度和区域尺度的直接检验方法以及交叉检验、多尺度逐级检验和时空变化趋势分析等间接检验方法，并指出像元尺度地表蒸散发相对真值的获取是其核心问题。现有的蒸散发相对真值获取的通量观测方法包括蒸渗仪方法、涡动相关技术、能量平衡波文比法、大孔径闪烁仪方法和微波闪烁仪方法。其中，涡动相关技术（EC）利用热通量、水汽变化和风速间的协方差来观测显热通量和潜热通量，是目前全球使用最广泛的通量观测方法，被广泛应用于蒸散发精度验证研究中。EC 得到的通量其观测源区与风速风向、仪器观测高度和下垫面状况等因素有关，受大气湍流的影响而时刻变化，EC 观测源区大约在几十米到几百米之间（Jia 等，2012）。对于高分辨率蒸散发估算，可以考虑通量观测源区来进行源区内多个像元的精度验证（Jia 等，2012；Bai 等，2015；Liu 等，2016）；针对中低分辨率蒸散发估算，目前主要是站点观测作为像元相对真值来直接验证（Jung 等，2010；Mu 等，2011；Tang 和Li，2017；Zhang 等，2019），或者用高分辨率蒸散发估算结果作为区域真值进行验证（Xu等，2018；Li等，2018）。对于公里尺度的蒸散发验证，Mu 等（2011）、Yao 等（2014）、Jiang 和Ryu（2016）、Zhang 等（2019）使用全球FLUXNET 通量观测数据在站点尺度分别对MOD16算法估算ET（1 km/8 d）、GLASS 算法估算ET（1 km/8 d）、BESS产品ET（1 km/8 d）以及PML−V2算法估算ET（500 m/8 d）进行了直接验证，估算ET和观测ET的RMSE分别为：0.9 mm/d、35.3 W/m2、0.78 mm/d、0.69 mm/d。

农田蒸散长期以来一直是农业、气象和水文等学科的研究焦点（姚云军 等，2012）。在农业旱情监测和灌溉管理中，蒸散发可以作为衡量作物缺水的指标。高精度的遥感反演农田蒸散发对于田间尺度的更精准的水分平衡量化和水分亏缺研究具有重要意义，对于农业精准灌溉和提高农田水分利用效率具有实用价值。在已发布的蒸散发产品中，MOD16 产品和PML−V2 产品空间分辨率（500 m）最接近于农田站点EC 观测源区大小，且500 m 像元上农田地表相对均一，空间代表性较好。对于均匀地表，通量站点观测值可被当作遥感像元上的相对真值来检验观测仪器所在蒸散产品像元值或者周围多个像元上的算术/加权平均值（张圆 等，2020）。故本研究选取这两种蒸散发产品来进行农田蒸散发的对比验证。

2 数据介绍

2.1 ET产品

2.1.1 MOD16产品

NASA 发布的MOD16 产品为用户分别提供了ET、LE、潜在ET 和潜在LE。冠层截留的水分蒸发是高叶面积指数（LAI）生态系统重要的水分通量，MOD16 算法 基于Penman−Monteith 方程，将ET 划分为湿冠层表面蒸发、干冠层表面植被蒸腾和土壤表面蒸发3部分，并在白天和夜晚分开估算求和。该算法将土壤表面分为饱和湿润表面和潮湿表面，并通过LAI将气孔导度转化为冠层导度来计算植被蒸腾，同时针对不同的地表覆被类型的阻抗分别参数化，以提高特定地表下垫面ET 的遥感反演精度（Mu等，2011）。相比于Mu等（2011），MOD16 V6产品生产进行了以下改进：如果白天净辐射为负值时，使其为0；对土壤热通量进行了限制保证土壤蒸发不为负值；针对持续严重的有云天导致的非洲西部雨林区没有有效的反照率的问题，设定其反照率为定值0.4；更新了不同地表覆被类型的生物群落属性查找表参数值；用3年平滑地表覆被产品MCDLCHKM 替换年地表覆被产品MOD12Q1。MOD16 产品生产使用的输入数据为GMAO MERRA 气象再分析数据（1/2°×2/3°）和MODIS 遥感产品（500 m）。MOD16 产品包括8 天、月和年3 种时间尺度的产品，本研究所采用的MOD16 产品为8 天合成的500 m 空间分辨率的MOD16A2 V6.1产品ET（Running等，2021）。

2.1.2 PML-V2产品

PML−V2产品是利用Penman−Monteith−Leuning模型生产的水热耦合的产品（张永强 等，2021），将ET 分为土壤蒸发、植被蒸腾和冠层截留蒸发3部分。其中，冠层截留蒸发通过适用性分析模型来估算，植被蒸腾和土壤蒸发通过LAI 来划分，PML−V1 模型结合Penman−Monteith 公式和Leuning导度模型估算植被蒸腾，在此基础上，PML−V2模型利用气孔导度耦合植被蒸腾和光合同化速率，采用碳水耦合冠层导度模型估算植被蒸腾和总初级生产力（GPP），使得ET 和GPP 得以相互制约，通过在全球分布的95 个EC 通量观测站分10 种植被类型分别进行参数率定，提高ET 和GPP 估算精度（Zhang等，2019）。PML−V2产品生产所采用的驱动数据包括全球陆面数据同化系统GLDAS（Global Land Data Assimilation System）气象数据（0.25°、3 h）和MODIS 遥感产品（500 m），将在通量观测站点上率定的模型参数由站点尺度扩展到全球逐网格尺度，从而实现全球ET 和GPP 产品生产。不同于MOD16 仅提供ET，PML−V2 产品分别提供土壤蒸发、植被蒸腾和冠层截留蒸发3部分。

2.2 站点观测数据

本论文所采用的用于农田类型蒸散发验证的通量观测站点包括20个全球FLUXNET 农田通量观测站点（表1）和8 个中国农田通量观测站点（表2），所采用的通量观测数据为日通量观测数据。FLUXNET 数据提供了日观测通量质量控制数据（QC），QC 的取值范围为0—1，表示每天半小时观测数据或者高质量插补数据所占的比例（Pastorello 等，2020）。参考Mu 等（2011），Yao 等（2014）、Hu 和Jia（2015）、Tang 等（2015）、Zhang等（2019）的研究，本研究在将8天内的日尺度ET观测累加获取8天尺度ET观测。在求8天累积ET观测时，对8天内的每天的观测质量控制数据求平均，当8 天内观测数据或者高质量插补数据大于80%（即平均QC大于0.8）才认为8天观测数据有效。中国通量观测数据未提供QC，在求8 天累积ET 观测时，8天内的每天都有观测值才认为8天观测数据有效。表1中将MOD16算法（Mu等，2011）和PML−V2算法（Zhang等，2019）曾使用过的算法验证站点分别用△和◇标识。虽然有11个FLUXNET站点在Zhang 等（2019）中使用，但是该研究是采用留一法（即每次使用10个农田站点来率定模型系数，将率定后的系数用在剩余一个站点上，对该站点上的算法精度进行验证）进行算法精度的验证，而不是对实际的PML−V2产品精度的验证。

表1 20个FLUXNET农田通量观测站点信息Table 1 Information of 20 tower sites for flux observation covered by cropland in FLUXNET2015

表2 8个中国农田通量观测站点信息Table 2 Information of eight tower sites for flux observation covered by cropland in China

2.3 其他辅助数据

本论文中使用的其他辅助数据包括土地利用/覆被数据和地面高程信息。土地利用/覆被数据使用的是500 m 空间分辨率的MODIS 地表覆被数据（MCD12Q1），采用其中的国际地圈—生物圈计划（IGBP）分类方案，该方案中定义了17 种不同的土地覆被类型，包括11种自然植被类别、3种人为改变类别以及3种非植被类别。本论文中使用的地面高程信息来源于ASTER GDEM数据，该数据是目前唯一覆盖全球陆地表面的高分辨率高程影像数据，其空间分辨率为30 m，垂直分辨率为1 m，数据可从NASA EARTH DATA下载获取，下载地址为：https://search.earthdata.nasa.gov/search[2022−01−10]。本论文所采用的遥感卫星影像数据来源于Google Earth，影像重采样分辨率为0.5 m，所用数据影像时间已在图1中标注。

图1 28个农田通量观测站点所在MODIS像元及周围3 × 3像元（500 m）Fig.1 The satellite images from Google Earth of the surrounding 3 × 3 MODIS pixels（500 m）for the 28 flux tower sites

2.4 评价指标

本论文中使用到的评价指标包括反映了产品8 天ET 与EC 观测8 天ET 互相关联程度的拟合度R2、反映产品8天ET 与EC 观测8天ET 偏离程度的平均偏差bias 和均方根误差RMSE。不同评价指标的计算方法如下：

式中，n代表观测8 天ET 或产品8 天ET 的样本数量，ETi是第i个8 天的累积ET 产品值，ETobs，i表示第i个8天的累积ET观测值。

3 精度对比验证结果

3.1 直接验证

在开展遥感产品真实性检验前，应先评价地表的空间异质性，对于相对均一的地表，可将站点观测值作为像元尺度相对真值，用于观测仪器所在像元的蒸散发产品值的验证（张圆 等，2020）。

3.1.1 地表均一性评价

为了开展ET 产品精度的综合评估，首先进行了通量观测站点的地表均一性评价。图1（a）、（b）、（c）分别展示了28个农田通量观测站点所对应的500 m 分辨率的MODIS单像元及周围3像元×3像元的地表类型分布、高程变化和从Google Earth上获取的卫星遥感影像（一个格子代表一个像元）。对每个站点，图1（a）中选用已获取通量观测数据时间范围内最新一年MCD12Q1的数据；图1（c）中尽量在Google Earth 上选取已获取通量观测数据年份的卫星遥感影像，如果观测年份没有可获取影像时，才选择其他时间的影像作为辅助检验；每个站点所用的地表分类数据及卫星遥感影像时间均标注在图上。28 个农田站点中，BE−Lon、DE−Geb、DE−Kli、FR−Gri、IT−Bci、US−ARM、US−CRT、US−Ne1、US−Ne2、US−Ne3、LC、CW、GT、YC 和SX 共15 个观测站点在站点附近3×3 个像元内MODIS地表覆被类型完全一致，均为农田。DK−Fou站点所在像元地表覆被类型为农田/自然植被混合，从Google Earth 上可以看出该像元地表覆被类型多样化，地表异质性较明显。虽然CH−Oe2、DE−Seh、IT−CA2、US−Lin、US−Tw3、US−Twt 和MY 等7 个站点的所在像元MODIS 地表覆被类型是农田，但是在站点附近3×3个像元内含有非农田类型的像元，农田像元所占比例依次为67%、78%、89%、67%、78%、67%、33%，除MY站点和CH−Oe2站点（左上角像元内含有建设用地）外的其余站点上，MODIS 地表覆被类型显示不是农田的像元在Google Earth 影像上显示像元大面积上基本均为农田。尽管MODIS 地表覆被数据显示DE−RuS、FI−Jok、US−Tw2、CN−Du1 和TW−Tar 等5 个站点（在表1 和表2 中用*标识）其站点所在像元地表覆被类型不是农田，依次为建设用地、稀树草原、草地、草地和稀树草原，但是在Google Earth 影像上显示这些像元均为农田地类。值得注意的是，MOD16 不计算城市区域，故而MOD16 产品在该像元上没有数据；US−Tw2 站点所在像元下垫面植被覆被类型14 年之前为农田，14 年之后为湿地。对不同下垫面，MOD16 产品和PML−V2 产品其算法参数化方案均不同，MODIS 地表覆被类型的错分，使得对原本的农田地类像元按照非农田地类的参数化方案开展该像元ET 的估算和产品的生产，使得这些像元上两种产品的ET 值存在较大的不确定性。

3.1.2 验证结果

DE−Rus 站点所对应的地表覆被类型为建筑用地，MOD16 产品不提供该地类的ET 值，故图2 展示了19 个FLUXNET 农田站点上MOD16 产品和PML−V2 产品8 天累积ET 的EC 观测直接验证结果。本研究采用融合了核密度估计曲线、箱线图及抖动散点图的云雨图来表示两种产品8 天ET 的偏差分布。图3 展示了19 个FLUXNET 农田站点上MOD16产品和PML−V2产品8天累积ET与EC 观测8天ET偏差的云雨图。从整体上来看，PML−V2产品略高估8 天ET（bias：0.66 mm/8 d），MOD16 产品略低估8 天ET（bias：−1.79 mm/8 d），PML−V2产品精度（RMSE：8.85 mm/8 d）与MOD16 产品精度（RMSE：8.54 mm/8 d）相差不多；12个站点上PML−V2产品精度优于MOD16产品精度，另外7个站点上MOD16 产品精度优于PML−V2 产品精度；在BE−Lon、DE−Geb、DE−Kli、DK−Fou、FI−Jok和IT−Ca2 6 个站点，MOD16 产品（bias：0.67—15.05 mm/8 d）和PML−V2 产 品（bias：1.53—13.27 mm/8 d）均高估8天ET；在CH−Oe2、DE−Seh、US−CRT、US−Tw2、US−Tw3 和US−Twt 等6 个站点上，MOD16 产品（bias：−16.42—−1.04 mm/8 d）和PML−V2 产品（bias：−15.98—−2.28 mm/8 d）均低估8天ET。

图2 19个FLUXNET农田站点MOD16产品和PML−V2产品8天ET验证散点图（95%置信区间）Fig.2 Observed ET versus MOD16 ET and PML−V2 ET for 8−day estimates at the 19 cropland sites in FLUXNET（95% confidence interval）

图3 19个FLUXNET农田站点MOD16产品和PML−V2产品与EC观测8天ET偏差的云雨图Fig.3 Raincloud plots of the biases between the 8−day ET（of the MOD16 or PML−V2 product）and the observed 8−day ET at the 19 cropland sites in FLUXNET

对于在站点附近3×3 个像元内MODIS 地表覆被类型完全一致均为农田的10 个站点：在BE−Lon、DE−Geb、DE−Kli、US−CRT、US−Ne1、US−Ne2 和US−Ne3 等7 个站点上，PML−V2 产品精度（RMSE：4.42—7.23 mm/8 d）优于或略优于MOD16产品精度（RMSE：7.22—8.33 mm/8 d）；在FR−Gri、IT−Bci 和US−ARM 3 个站点上，MOD16 产品精度（RMSE：4.24—6.07 mm/8 d）优于或略优于PML−V2产品精度（RMSE：5.22—11 mm/8 d）。对于其站点所在像元MODIS 地表覆被类型不是农田的3 个站点：站点所在像元为建筑用地的DE−RuS站点，MOD16 产品没有产品估算结果，故在图4和图5 中未予以评价；所在像元为稀树草原的FI−Jok 站点上，MOD16 产品精度较差RMSE 大于10 mm/8 d，PML−V2 产品RMSE 仍在5 mm/8 d 内，两种产品根据MODIS 地类分类计算，计算ET 时按稀树草原的参数计算，可能是造成估算误差大的原因；所在像元为草地的US−Tw2 站点，尽管所处位置非常靠近下方农田像元，但产品值应该是按草地的参数计算的草地ET 而非农田ET，两种产品 均低估8 天ET，PML−V2 产品精度（RMSE：4.45 mm/8 d，bias：−2.72 mm/8 d）优于MOD16 产品精度（RMSE：8.06 mm/8 d，bias：−5.35 mm/8 d）。对于虽然站点所在像元MODIS 地表覆被类型是农田，但是在站点附近3×3个像元内含有非农田类型像元的6 个站点：在CH−Oe2、DE−Seh 和IT−CA2等3 个站点上，两种产品RMSE 均在10 mm/8 d 以内，MOD16产品精度（RMSE：3.81—6.38 mm/8 d，bias：−2.78—0.67 mm/8 d）优于或略优于PML−V2产品精度（RMSE：5.75—9.58 mm/8 d，bias：−4.34—4.25 mm/8 d），在IT−CA2 站点上，MOD16 产品大 部分8 天累积ET 与观测8 天ET 偏差在−10—10 mm/8 d，高估和低估的分布相对均匀，而PML−V2产品部分8天累积ET高估可达将近30 mm/8 d；在US−Tw3 和US−Lin两个站点上，PML−V2 产品（RMSE：5.57—8.03 mm/8 d，bias：−5.39—0.66 mm/8 d）精度略优于MOD16 产品精度（RMSE：6.81—11.03 mm/8 d，bias：−3.44—−8.55 mm/8 d）；地表异质性较明显的US−Twt 站点上，MOD16 产品（RMSE：21.47 mm/8 d，bias：−16.42 mm/8 d）和PML−V2产品（RMSE：22.4 mm/8 d，bias：−15.98 mm/8 d）精度均较差，严重低估8 天ET，该站点是唯一的水稻站点，证明了现有产品对水稻ET的严重低估。

图4 中国8个农田站点MOD16产品和PML−V2产品8天ET验证（95%置信区间）Fig.4 Observed ET versus MOD16 ET and PML−V2 ET for 8−day estimates at the 8 cropland sites in China（95% confidence interval）

图5 中国8个农田站点MOD16产品和PML−V2产品与EC观测8天ET偏差云雨图Fig.5 Raincloud plots of the biases between the 8−day ET（of the MOD16 or PML−V2 product）and the observed 8−day ET at the 8 cropland sites in China

图4和图5分别展示了中国8个农田站点的8天ET 直接验证结果，以及MOD16 产品和PML−V2 产品8天累积ET与EC观测8天ET偏差的云雨图。在LC、YC、GT和MY这4个站点上，PML−V2产品8天ET精度（R2：0.271—0.811，RMSE：3.3—9.2 mm/8 d）明显优 于MOD16 产 品8 天ET 精度（R2：0.018—0.600，RMSE：9.74—13.45 mm/8 d），PML−V2 产品在GT 站点高估8 天ET（bias：3.52 mm/8 d），在其余3个站点低估8天ET（bias：−1.72—−0.57 mm/8 d），MOD16产品在4个站点均低估8天ET（bias：−8.4—−6.21 mm/8 d）；在CN−Du1站点，两种产品8天ET验证精度相差不大且整体精度较其他站点相对较高，MOD16产品略低估8天ET（bias：−2.49 mm/8 d，RMSE：5.89 mm/8 d），而PML−V2 产品略高估8 天ET（bias：1.08 mm/8 d，RMSE：6.76 mm/8 d）；在CW 和SX两个站点，MOD16产品8天ET精度（R2：0.143—0.439，RMSE：5.73—10.12 mm/8 d，bias：−2.29—2.37 mm/8 d）优于PML−V2产品8天ET精度（R2：0.014—0.282，RMSE：10.08—11.25 mm/8 d，bias：7.42—8.04 mm/8 d）；TW−Tar 站点有效数据太少，统计意义不足，在仅有的3 个8 天，PML−V2产品略优于MOD16产品。

图6展示了中国8个农田站点和全球28个农田站点上MOD16 产品和PML−V2 产品8 天ET 总体验证结果。对于中国8 个农田站点，PML−V2 产品精度（R2：0.480，RMSE：8.46 mm/8 d）优于MOD16产品精 度（R2：0.407，RMSE：10.77 mm/8 d）；MOD16 产品在7 个站点上均低估8 天ET，仅在SX站点上高估8天ET，整体低估8 天ET（bias：−6.29 mm/8 d）；而PML−V2 产品在4 个站点上高估8天ET，在另外4个站点上低估8天ET，整体略低估8 天ET（bias：−0.66 mm/8 d）。对于本研究中评估的全球共28个农田站点，MOD16产品精度（R2：0.452，RMSE：8.82 mm/8 d）和PML−V2 产品精度（R2：0.455，RMSE：8.81 mm/8 d）整体精度相当，MOD16 产品整体略低估8 天ET（bias：−2.31 mm/8 d），而PML−V2 产品整体略高估8 天ET（bias：0.51 mm/8 d）。

图6 MOD16产品和PML−V2产品8天ET整体验证结果Fig.6 The overall validation of MOD16 ET and PML−V2 ET for 8−day estimates

3.2 趋势分析检验

图7展示了20个FLUXNET 农田站点上MOD16产品和PML−V2 产品8 天累积ET 与EC 观测8 天累积ET 的时序变化。在8 天中EC 观测ET 其QC 显示观测和高质量插补数据不足80%的在图中用红色*表示。该图7 中R2、RMSE 与图2 中R2和图3 中RMSE的差别在于，该图7中对MOD16产品有效数据和PML−V2有效数据分别计算R2和RMSE，而图2和图3 中只展示了MOD16 产品和PML−V2 产品数据值均有效的8 天，其计算R2和RMSE 相应的8 天的天数两种产品是一致的。整体而言，PML−V2产品时间连续性优于MOD16 产品，相比于考虑两种产品有效数据时间一致的情况（图2和图3），大多数站点上PML−V2产品精度有所提升，在15个站点上，PML−V2产品精度优于MOD16产品精度，但是MOD16产品可以相对更好的捕捉到年中因为作物轮种而引起的ET下降和上升趋势，而PML−V2产品普遍在年中达到当年最高峰值。对于BE−Lon、DE−Geb和DE−Kli等3个站点，PML−V2产品与EC观测ET二者曲线的吻合度相对更好（R2：0.78—0.8），而MOD16产品（R2：0.54—0.62）几乎在每年的年中波峰处存在严重高估ET 的现象；对于CE−Oe2 站点，MOD16产品和PML−V2产品几乎在每年ET高值处存在低估ET 的现象；在DK−Fou 站点，MOD16 和PML−V2 两种产品时间变化特点相似（二者R2为0.85，趋势均为先上升，在6月下旬至7月上旬达到波峰，而后下降），相比EC观测ET（二者与观测R2分别为0.28和0.31，趋势先上升，在8月下旬至9月上旬达到波峰，而后下降）而言均高估ET（bias>10 W/m2）；在US−Lin站点，两种产品均未捕捉到观测ET 峰值，相比观测ET 提前达到各自峰值，且PML−V2产品ET整体高于MOD16产品ET；US−Ne1、US−Ne2和US−Ne3等3个以玉米作物为主的站点，虽然相比于观测8 天ET，PML−V2 产品整体精度高于MOD16产品，但MOD16产品8天ET对于细节的捕捉能力略优于PML−V2产品，在上半年MOD16产品与EC观测ET二者曲线的吻合度更好，而PML−V2产品明显高估EC观测8天ET。以US−Ne1站点为例，在3—6月，从与EC观测ET时序曲线的吻合度上来看，MOD16产品（R2：0.70，bias：−4.02 mm/8 d，RMSE：49.02 mm/8 d）略优于PML−V2产品（R2：0.64，bias：−4.53 mm/8 d，RMSE：58.82 mm/8 d）；在IT−Bci、IT−CA2 和US−ARM 等3 个站点，MOD16 产品与EC观测8天ET变化趋势更为相似，PML−V2产品高估8天ET，而MOD16产品可以相对较好的捕捉到年中的下降或上升趋势，如在IT−CA2站点上，MOD16产品在年中更好地反映了ET 先下降后上升的趋势；US−Tw2、US−Tw3 和US−Twt 等3 个站点上，两种产品均低估8 天ET 表现不佳，但是MOD16 产品仍能相对更好捕捉到年中ET 的波动（下降或上升趋势），US−Twt 站点上MOD16 产品对最高峰值所处时间的刻画明显优于PML−V2产品，MOD16产品和观测值均在下半年达到年度最高峰值，而PML−V2产品则提前在上半年达到年度最高峰值。

图7 20个FLUXNET农田站点MOD16产品、PML−V2产品和EC观测8天ET时序图Fig.7 Temporal evolution of the 8−day ET from MOD16，PML−V2 and observation at the 20 cropland sites in FLUXNET

图8展示了中国8个农田站点上MOD16产品和PML−V2 产品8 天累积ET 与EC 观测8 天累积ET 的时序变化。整体来看，PML−V2产品时间连续性优于MOD16产品，PML−V2产品8天ET与EC观测8天ET 整体吻合性更好（R2更高、RMSE 更低），而MOD16产品对部分站点时序细节的刻画优于PML−V2 产品。例如MOD16 产品可以捕捉到YC 站点和LC 站点每年度年中的低谷，刻画年中8 天ET 的先下降后上升趋势，而在年中EC 观测8天ET 处于低谷时PML−V2产品仍处于其年度峰值状态。

图8 中国8个农田站点MOD16产品、PML−V2产品和EC观测8天ET时序图Fig.8 Temporal evolution of the 8−day ET from MOD16，PML−V2 and observations at the 8 cropland sites in China

为了更具体展示两种产品对于时序曲线细节的刻画，本研究选取IT−CA2、US−ARM、US−NE1、US−Twt、LC 和YC 等6 个站点，对MOD16产品和PML−V2 产品8 天累积ET 与EC 观测8 天累积ET 开展了时序变化分解，分析其季节变化趋势（图9）。本研究采用基于贝叶斯集成的BEAST时序分解算法（Zhao 等，2019）来对ET 时序进行季节变化分解，相比于其他时序分解算法，该算法允许空缺值的存在，在MOD16 产品时序不完整的情况下，不需要进行插值处理就能获得较为理想的结果。此处主要是为了对ET 时序变化的波峰和波谷所处时刻进行探究，故只展示了分解后的季节性波动信号变化。从图9 中可以清晰看出，虽然MOD16 产品连续性较差，但是仍能较好的反映出年中ET的下降和上升趋势。在IT−CA2和US−ARM两个站点上，MOD16 产品与EC 观测ET 季节变化波动曲线贴合度更好（IT−CA2 和US−ARM 的R2分别为0.83 和0.76），均在上半年（4—5 月）达到波峰，而PML−V2 产品（IT−CA2 和US−ARM 的R2分别为0.43 和0.69）则在年中（6—7 月）达到波峰；在US−Ne1 站点上，由于MOD16 产品数据在年初和年末的缺失，使得其没有较好的捕捉到两年交替之际ET 时序曲线的低谷，从与EC 观测ET 时序曲线的吻合度上来看，MOD16产品（R2：0.92）整体略优于PML−V2产品（R2：0.88），在上半年3—6月，MOD16 产品（R2：0.91）优于PML−V2 产品（R2：0.82），MOD16 产品和EC 观测ET 约从5—6 月开始急剧上升，上升曲线的坡度更陡（MOD16产品在2—6 月、3—6 月、4—6 月、5—6 月、6 月的平均坡度分别为0.09、0.13，0.23、0.36、0.55，EC观测在2—6 月、3—6 月、4—6 月、5—6 月、6 月的平均坡度分别为0.18、0.23，0.32、0.45、0.63），而PML−V2产品ET约从3月开始稳步上升，上升曲线的坡度相对缓和（平均坡度一直维持在0.23—0.26范围内）；在US−Twt站点上，MOD16产品能够捕捉到EC观测ET的两处小高峰，并且在下半年达到其最高峰，尽管PML−V2产品在2012年—2014年也能够捕捉到两个波峰，但是在上半年达到最高峰，与EC观测ET和MOD16 产品相反；在LC 和YC 两个站点上，PML−V2产品均在年中达到年度峰值，未捕捉到冬小麦、夏玉米轮种而导致的年中ET先下降再上升的趋势，而MOD16产品成功捕捉到冬小麦和夏玉米生长期内的两个ET波峰，并在下半年达到最高峰（夏玉米生长季），与EC观测ET恰好相反，EC观测ET在上半年达到最高峰（冬小麦生长季），两个站点上MOD16产品均低估了冬小麦ET。

图9 选取的6个农田站点EC观测8天ET、MOD16产品和PML−V2产品8天ET时序分解季节波动信号Fig.9 The seasonal fluctuant signals decomposed from the 8−day time−series of MOD16 ET，PML−V2 ET and observed ET at the selected 6 cropland sites

4 结论

本文采用20 个FLUXNET 通量观测站点和8 个中国的通量观测站点共计28 个农田站点评估了当前国际上仅有的两种500 m 分辨率的8 天ET 产品：MOD16产品和PML−V2产品。整体而言，28个农田站点的验证表明MOD16产品（RMSE：8.82 mm/8 d）和PML−V2产品（RMSE：8.81 mm/8 d）总体精度相当；MOD16产品低估8天ET（bias：−2.31 mm/8 d），PML−V2 产品略高估8 天ET（bias：0.51 mm/8 d）；PML−V2 产品时间连续性优于MOD16 产品，但是部分站点MOD16 产品在时序变化细节（如到达峰值的季节、年中的下降上升趋势）的刻画上优于PML−V2产品。20个FLUXNET通量观测站点评估结果表明，两种产品精度相当，PML−V2产品（RMSE：4.42—22.4 mm/8 d，bias：−15.98—13.27 mm/8 d）略低估8天ET，MOD16产品（RMSE：3.81—21.48 mm/8 d，bias：−16.42—15.05 mm/8 d）略高估8 天ET；PML−V2产品时间连续性较好，在不考虑两种产品验证点数一致的情况下，PML−V2 产品精度（RMSE：4.06—20.36 mm/8 d）有所提升。PML−V2算法利用FLUXNET站点观测通量对农田地类进行了单独参数率定，采用留一法在11 个农田站点上进行了算法精度验证，结果表明对于农田地类该算法估算RMSE 变化范围为0.54—2.38 mm/d，R2变化范围为0.42—0.86，MOD16 产品RMSE 变化范围为0.48—2.4 mm/d，R2变化范围为0.49—0.76（Zhang 等，2019）。本研究与该研究中对农田地类的8 天ET 估算精度验证结果相当，均显示US−Twt站点上水稻ET 的估算结果最差。本文采用的8 个中国农田通量站点并未用于PML−V2 算法参数率定，但中国站点的验证结果表明两种产品均低估8天ET，PML−V2 产品（RMSE：2.88—11.25 mm/8 d，bias：−1.72—8.04 mm/8 d）仍具有相对可靠的精度，且精度优于MOD16 产品（RMSE：5.73—13.45 mm/8 d，bias：−8.4—2.37 mm/8 d），在不考虑两种产品验证点数一致的情况下，PML−V2产品精度（R2：0.28—0.87，RMSE：3.23—9.58 mm/8 d）有所提升，但MOD16 产品对年中8 天ET 的下降和上升趋势刻画优于PML−V2产品。

本文的研究仍存在不足，虽然结合MODIS 地表覆被产品和Google Earth 卫星遥感影像对这些站点附近的地类和地表均一性进行了判断，但是仍显粗糙，未来建议结合高空间分辨率的植被指数、叶面积指数等植被生态参数来进行同时期更精细的地表均一性评价。由于地表的空间异质性，地面观测和遥感估算的蒸散发之间存在着空间尺度不匹配的问题，如何从地面站点观测中获取遥感像元尺度地表蒸散发的相对真值，是目前蒸散发研究的热点和难点。对于站点像元空间异质性强的情况，需要考虑观测源区，且可能要进行空间尺度转换，但在实际操作中并不容易；站点所处500 m 像元的空间代表性或均一性较好，可以直接进行验证。已有研究报道，馆陶农田站点观测空间范围大约在450—460 m（Liu 等，2011），多伦农田站点的观测空间范围大约在300—400 m，寿县农田站点的观测空间范围大约在200—400 m（王有恒，2010）。出于对农田站点观测源区一般小于500 m 的假设，本研究只选取了500 m 空间分辨率的MOD16 产品和PML−V2 产品，未来应考虑结合其他时空分辨率产品进行交叉检验，如1 km空间分辨率的GLASS ET 产品、BESS ET 产品和ETMonitor ET 产品，并结合高时间分辨率的ET 产品如GLEAM 产品，进行更大尺度像元的时空趋势验证。另外，站点观测通量的30 min 到日尺度、8 天尺度、月尺度和年尺度的时间尺度转换方法亦是目前研究的难点。500 m 尺度像元上的地表空间相对均一，可以考虑直接简单平均，而对于异质地表，可以考虑通量观测源区的时间和空间加权来获取8 天地面观测ET，但是在实际操作时并非易事。这些难点需要在未来的研究中予以更好的考虑，探索更合理的空间尺度和时间尺度转化方法，以获取更可靠的蒸散发像元尺度验证结果。

志 谢感谢FLUXNET 全球通量观测研究网络、ChinaFLUX 中国通量观测研究网络、中国生态系统研究网络、国家生态科学数据中心、国家青藏高原数据中心、台湾农业研究所等研究网络和机构提供宝贵的农田站点通量观测数据，感谢中国科学院地理科学与资源研究所于贵瑞研究员、北京师范大学刘绍民教授、中国科学院青藏高原研究所李新研究员等科研学者共享数据及在推动观测网络建设和数据共享上作出的贡献，同时感谢各个农田站点工作人员在数据观测整理过程中付出的辛勤劳动，在此表示衷心的感谢！