劉萌，彭中，黃凌霄，李召良，段四波，唐榮林

1.中國農(nóng)業(yè)科學院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所/北方干旱半干旱耕地高效利用全國重點實驗室,北京 100081;

2.中國科學院地理科學與資源研究所 資源與環(huán)境信息系統(tǒng)國家重點實驗室,北京 100101;

3.中國科學院大學,北京 100049

1 引言

地表蒸散發(fā)ET（Evapotranspiration）是指從地表傳輸?shù)酱髿庵械乃?，包括土壤表面水分蒸發(fā)、植被水分蒸騰和冠層截留蒸發(fā)。水分蒸發(fā)的過程需要吸收能量來冷卻地表，這部分能量被稱為潛熱通量LE（Latent heat flux ），是蒸散發(fā)的能量表達形式，二者可以相互轉(zhuǎn)化（Wang 和Dickinson，2012）。蒸散發(fā)和潛熱連接了水循環(huán)和能量循環(huán)，是地表水熱循環(huán)中的關(guān)鍵過程參量，也是其中最難估算的地表特征參量。遙感技術(shù)因其大范圍觀測能力，為獲取豐富的區(qū)域地表特征參量提供了可能，遙感反演蒸散發(fā)是獲取區(qū)域、國家和全球尺度蒸散發(fā)時空分布的有效方法（Li 等，2009）。目前，國內(nèi)外學者已經(jīng)先后生產(chǎn)發(fā)布了不同時空分辨率的區(qū)域及全球蒸散發(fā)遙感反演產(chǎn)品，如MOD16產(chǎn)品（Mu 等，2011）、LSA?SAF 產(chǎn)品（Ghilain 等，2011）、SSEBop產(chǎn)品（Senay等，2013）、GLASS產(chǎn)品（Yao等，2014）、ETMonitor產(chǎn)品（Hu和Jia，2015）、BESS 產(chǎn) 品（Jiang 和Ryu，2016）、GLEAM 產(chǎn) 品（Martens 等，2017）、PML?V2 產(chǎn)品（Zhang 等，2019）、FLUXCOM 產(chǎn)品（Jung 等，2019）和CR?ET 產(chǎn)品（Ma 等，2019）等。相關(guān)研究已經(jīng)詳細介紹了現(xiàn)有的過程驅(qū)動的蒸散發(fā)遙感產(chǎn)品（李佳 等，2021）和數(shù)據(jù)驅(qū)動的蒸散發(fā)遙感產(chǎn)品（劉萌 等，2021）。由于不同產(chǎn)品之間存在反演方法機理差異、產(chǎn)品特征差異等，在使用這些產(chǎn)品前，需要對產(chǎn)品的精度進行驗證和評估（熊育久 等，2021）。然而，由于近地表大氣的復雜多變性、地表下墊面的異質(zhì)性等特點，蒸散發(fā)產(chǎn)品的驗證目前仍面臨諸多困難和挑戰(zhàn)。張圓等（2020）詳細綜述了蒸散發(fā)真實性檢驗的發(fā)展現(xiàn)狀和研究進展，將蒸散發(fā)真實性檢驗方法分為像元尺度和區(qū)域尺度的直接檢驗方法以及交叉檢驗、多尺度逐級檢驗和時空變化趨勢分析等間接檢驗方法，并指出像元尺度地表蒸散發(fā)相對真值的獲取是其核心問題。現(xiàn)有的蒸散發(fā)相對真值獲取的通量觀測方法包括蒸滲儀方法、渦動相關(guān)技術(shù)、能量平衡波文比法、大孔徑閃爍儀方法和微波閃爍儀方法。其中，渦動相關(guān)技術(shù)（EC）利用熱通量、水汽變化和風速間的協(xié)方差來觀測顯熱通量和潛熱通量，是目前全球使用最廣泛的通量觀測方法，被廣泛應用于蒸散發(fā)精度驗證研究中。EC 得到的通量其觀測源區(qū)與風速風向、儀器觀測高度和下墊面狀況等因素有關(guān)，受大氣湍流的影響而時刻變化，EC 觀測源區(qū)大約在幾十米到幾百米之間（Jia 等，2012）。對于高分辨率蒸散發(fā)估算，可以考慮通量觀測源區(qū)來進行源區(qū)內(nèi)多個像元的精度驗證（Jia 等，2012；Bai 等，2015；Liu 等，2016）；針對中低分辨率蒸散發(fā)估算，目前主要是站點觀測作為像元相對真值來直接驗證（Jung 等，2010；Mu 等，2011；Tang 和Li，2017；Zhang 等，2019），或者用高分辨率蒸散發(fā)估算結(jié)果作為區(qū)域真值進行驗證（Xu等，2018；Li等，2018）。對于公里尺度的蒸散發(fā)驗證，Mu 等（2011）、Yao 等（2014）、Jiang 和Ryu（2016）、Zhang 等（2019）使用全球FLUXNET 通量觀測數(shù)據(jù)在站點尺度分別對MOD16算法估算ET（1 km/8 d）、GLASS 算法估算ET（1 km/8 d）、BESS產(chǎn)品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。

農(nóng)田蒸散長期以來一直是農(nóng)業(yè)、氣象和水文等學科的研究焦點（姚云軍 等，2012）。在農(nóng)業(yè)旱情監(jiān)測和灌溉管理中，蒸散發(fā)可以作為衡量作物缺水的指標。高精度的遙感反演農(nóng)田蒸散發(fā)對于田間尺度的更精準的水分平衡量化和水分虧缺研究具有重要意義，對于農(nóng)業(yè)精準灌溉和提高農(nóng)田水分利用效率具有實用價值。在已發(fā)布的蒸散發(fā)產(chǎn)品中，MOD16 產(chǎn)品和PML?V2 產(chǎn)品空間分辨率（500 m）最接近于農(nóng)田站點EC 觀測源區(qū)大小，且500 m 像元上農(nóng)田地表相對均一，空間代表性較好。對于均勻地表，通量站點觀測值可被當作遙感像元上的相對真值來檢驗觀測儀器所在蒸散產(chǎn)品像元值或者周圍多個像元上的算術(shù)/加權(quán)平均值（張圓 等，2020）。故本研究選取這兩種蒸散發(fā)產(chǎn)品來進行農(nóng)田蒸散發(fā)的對比驗證。

2 數(shù)據(jù)介紹

2.1 ET產(chǎn)品

2.1.1 MOD16產(chǎn)品

NASA 發(fā)布的MOD16 產(chǎn)品為用戶分別提供了ET、LE、潛在ET 和潛在LE。冠層截留的水分蒸發(fā)是高葉面積指數(shù)（LAI）生態(tài)系統(tǒng)重要的水分通量，MOD16 算法 基于Penman?Monteith 方程，將ET 劃分為濕冠層表面蒸發(fā)、干冠層表面植被蒸騰和土壤表面蒸發(fā)3部分，并在白天和夜晚分開估算求和。該算法將土壤表面分為飽和濕潤表面和潮濕表面，并通過LAI將氣孔導度轉(zhuǎn)化為冠層導度來計算植被蒸騰，同時針對不同的地表覆被類型的阻抗分別參數(shù)化，以提高特定地表下墊面ET 的遙感反演精度（Mu等，2011）。相比于Mu等（2011），MOD16 V6產(chǎn)品生產(chǎn)進行了以下改進：如果白天凈輻射為負值時，使其為0；對土壤熱通量進行了限制保證土壤蒸發(fā)不為負值；針對持續(xù)嚴重的有云天導致的非洲西部雨林區(qū)沒有有效的反照率的問題，設定其反照率為定值0.4；更新了不同地表覆被類型的生物群落屬性查找表參數(shù)值；用3年平滑地表覆被產(chǎn)品MCDLCHKM 替換年地表覆被產(chǎn)品MOD12Q1。MOD16 產(chǎn)品生產(chǎn)使用的輸入數(shù)據(jù)為GMAO MERRA 氣象再分析數(shù)據(jù)（1/2°×2/3°）和MODIS 遙感產(chǎn)品（500 m）。MOD16 產(chǎn)品包括8 天、月和年3 種時間尺度的產(chǎn)品，本研究所采用的MOD16 產(chǎn)品為8 天合成的500 m 空間分辨率的MOD16A2 V6.1產(chǎn)品ET（Running等，2021）。

2.1.2 PML-V2產(chǎn)品

PML?V2產(chǎn)品是利用Penman?Monteith?Leuning模型生產(chǎn)的水熱耦合的產(chǎn)品（張永強 等，2021），將ET 分為土壤蒸發(fā)、植被蒸騰和冠層截留蒸發(fā)3部分。其中，冠層截留蒸發(fā)通過適用性分析模型來估算，植被蒸騰和土壤蒸發(fā)通過LAI 來劃分，PML?V1 模型結(jié)合Penman?Monteith 公式和Leuning導度模型估算植被蒸騰，在此基礎上，PML?V2模型利用氣孔導度耦合植被蒸騰和光合同化速率，采用碳水耦合冠層導度模型估算植被蒸騰和總初級生產(chǎn)力（GPP），使得ET 和GPP 得以相互制約，通過在全球分布的95 個EC 通量觀測站分10 種植被類型分別進行參數(shù)率定，提高ET 和GPP 估算精度（Zhang等，2019）。PML?V2產(chǎn)品生產(chǎn)所采用的驅(qū)動數(shù)據(jù)包括全球陸面數(shù)據(jù)同化系統(tǒng)GLDAS（Global Land Data Assimilation System）氣象數(shù)據(jù)（0.25°、3 h）和MODIS 遙感產(chǎn)品（500 m），將在通量觀測站點上率定的模型參數(shù)由站點尺度擴展到全球逐網(wǎng)格尺度，從而實現(xiàn)全球ET 和GPP 產(chǎn)品生產(chǎn)。不同于MOD16 僅提供ET，PML?V2 產(chǎn)品分別提供土壤蒸發(fā)、植被蒸騰和冠層截留蒸發(fā)3部分。

2.2 站點觀測數(shù)據(jù)

本論文所采用的用于農(nóng)田類型蒸散發(fā)驗證的通量觀測站點包括20個全球FLUXNET 農(nóng)田通量觀測站點（表1）和8 個中國農(nóng)田通量觀測站點（表2），所采用的通量觀測數(shù)據(jù)為日通量觀測數(shù)據(jù)。FLUXNET 數(shù)據(jù)提供了日觀測通量質(zhì)量控制數(shù)據(jù)（QC），QC 的取值范圍為0—1，表示每天半小時觀測數(shù)據(jù)或者高質(zhì)量插補數(shù)據(jù)所占的比例（Pastorello 等，2020）。參考Mu 等（2011），Yao 等（2014）、Hu 和Jia（2015）、Tang 等（2015）、Zhang等（2019）的研究，本研究在將8天內(nèi)的日尺度ET觀測累加獲取8天尺度ET觀測。在求8天累積ET觀測時，對8天內(nèi)的每天的觀測質(zhì)量控制數(shù)據(jù)求平均，當8 天內(nèi)觀測數(shù)據(jù)或者高質(zhì)量插補數(shù)據(jù)大于80%（即平均QC大于0.8）才認為8天觀測數(shù)據(jù)有效。中國通量觀測數(shù)據(jù)未提供QC，在求8 天累積ET 觀測時，8天內(nèi)的每天都有觀測值才認為8天觀測數(shù)據(jù)有效。表1中將MOD16算法（Mu等，2011）和PML?V2算法（Zhang等，2019）曾使用過的算法驗證站點分別用△和◇標識。雖然有11個FLUXNET站點在Zhang 等（2019）中使用，但是該研究是采用留一法（即每次使用10個農(nóng)田站點來率定模型系數(shù)，將率定后的系數(shù)用在剩余一個站點上，對該站點上的算法精度進行驗證）進行算法精度的驗證，而不是對實際的PML?V2產(chǎn)品精度的驗證。

表1 20個FLUXNET農(nóng)田通量觀測站點信息Table 1 Information of 20 tower sites for flux observation covered by cropland in FLUXNET2015

表2 8個中國農(nóng)田通量觀測站點信息Table 2 Information of eight tower sites for flux observation covered by cropland in China

2.3 其他輔助數(shù)據(jù)

本論文中使用的其他輔助數(shù)據(jù)包括土地利用/覆被數(shù)據(jù)和地面高程信息。土地利用/覆被數(shù)據(jù)使用的是500 m 空間分辨率的MODIS 地表覆被數(shù)據(jù)（MCD12Q1），采用其中的國際地圈—生物圈計劃（IGBP）分類方案，該方案中定義了17 種不同的土地覆被類型，包括11種自然植被類別、3種人為改變類別以及3種非植被類別。本論文中使用的地面高程信息來源于ASTER GDEM數(shù)據(jù)，該數(shù)據(jù)是目前唯一覆蓋全球陸地表面的高分辨率高程影像數(shù)據(jù)，其空間分辨率為30 m，垂直分辨率為1 m，數(shù)據(jù)可從NASA EARTH DATA下載獲取，下載地址為：https://search.earthdata.nasa.gov/search[2022?01?10]。本論文所采用的遙感衛(wèi)星影像數(shù)據(jù)來源于Google Earth，影像重采樣分辨率為0.5 m，所用數(shù)據(jù)影像時間已在圖1中標注。

圖1 28個農(nóng)田通量觀測站點所在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 評價指標

本論文中使用到的評價指標包括反映了產(chǎn)品8 天ET 與EC 觀測8 天ET 互相關(guān)聯(lián)程度的擬合度R2、反映產(chǎn)品8天ET 與EC 觀測8天ET 偏離程度的平均偏差bias 和均方根誤差RMSE。不同評價指標的計算方法如下：

式中，n代表觀測8 天ET 或產(chǎn)品8 天ET 的樣本數(shù)量，ETi是第i個8 天的累積ET 產(chǎn)品值，ETobs，i表示第i個8天的累積ET觀測值。

3 精度對比驗證結(jié)果

3.1 直接驗證

在開展遙感產(chǎn)品真實性檢驗前，應先評價地表的空間異質(zhì)性，對于相對均一的地表，可將站點觀測值作為像元尺度相對真值，用于觀測儀器所在像元的蒸散發(fā)產(chǎn)品值的驗證（張圓 等，2020）。

3.1.1 地表均一性評價

為了開展ET 產(chǎn)品精度的綜合評估，首先進行了通量觀測站點的地表均一性評價。圖1（a）、（b）、（c）分別展示了28個農(nóng)田通量觀測站點所對應的500 m 分辨率的MODIS單像元及周圍3像元×3像元的地表類型分布、高程變化和從Google Earth上獲取的衛(wèi)星遙感影像（一個格子代表一個像元）。對每個站點，圖1（a）中選用已獲取通量觀測數(shù)據(jù)時間范圍內(nèi)最新一年MCD12Q1的數(shù)據(jù)；圖1（c）中盡量在Google Earth 上選取已獲取通量觀測數(shù)據(jù)年份的衛(wèi)星遙感影像，如果觀測年份沒有可獲取影像時，才選擇其他時間的影像作為輔助檢驗；每個站點所用的地表分類數(shù)據(jù)及衛(wèi)星遙感影像時間均標注在圖上。28 個農(nóng)田站點中，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 個像元內(nèi)MODIS地表覆被類型完全一致，均為農(nóng)田。DK?Fou站點所在像元地表覆被類型為農(nóng)田/自然植被混合，從Google Earth 上可以看出該像元地表覆被類型多樣化，地表異質(zhì)性較明顯。雖然CH?Oe2、DE?Seh、IT?CA2、US?Lin、US?Tw3、US?Twt 和MY 等7 個站點的所在像元MODIS 地表覆被類型是農(nóng)田，但是在站點附近3×3個像元內(nèi)含有非農(nóng)田類型的像元，農(nóng)田像元所占比例依次為67%、78%、89%、67%、78%、67%、33%，除MY站點和CH?Oe2站點（左上角像元內(nèi)含有建設用地）外的其余站點上，MODIS 地表覆被類型顯示不是農(nóng)田的像元在Google Earth 影像上顯示像元大面積上基本均為農(nóng)田。盡管MODIS 地表覆被數(shù)據(jù)顯示DE?RuS、FI?Jok、US?Tw2、CN?Du1 和TW?Tar 等5 個站點（在表1 和表2 中用*標識）其站點所在像元地表覆被類型不是農(nóng)田，依次為建設用地、稀樹草原、草地、草地和稀樹草原，但是在Google Earth 影像上顯示這些像元均為農(nóng)田地類。值得注意的是，MOD16 不計算城市區(qū)域，故而MOD16 產(chǎn)品在該像元上沒有數(shù)據(jù)；US?Tw2 站點所在像元下墊面植被覆被類型14 年之前為農(nóng)田，14 年之后為濕地。對不同下墊面，MOD16 產(chǎn)品和PML?V2 產(chǎn)品其算法參數(shù)化方案均不同，MODIS 地表覆被類型的錯分，使得對原本的農(nóng)田地類像元按照非農(nóng)田地類的參數(shù)化方案開展該像元ET 的估算和產(chǎn)品的生產(chǎn)，使得這些像元上兩種產(chǎn)品的ET 值存在較大的不確定性。

3.1.2 驗證結(jié)果

DE?Rus 站點所對應的地表覆被類型為建筑用地，MOD16 產(chǎn)品不提供該地類的ET 值，故圖2 展示了19 個FLUXNET 農(nóng)田站點上MOD16 產(chǎn)品和PML?V2 產(chǎn)品8 天累積ET 的EC 觀測直接驗證結(jié)果。本研究采用融合了核密度估計曲線、箱線圖及抖動散點圖的云雨圖來表示兩種產(chǎn)品8 天ET 的偏差分布。圖3 展示了19 個FLUXNET 農(nóng)田站點上MOD16產(chǎn)品和PML?V2產(chǎn)品8天累積ET與EC 觀測8天ET偏差的云雨圖。從整體上來看，PML?V2產(chǎn)品略高估8 天ET（bias：0.66 mm/8 d），MOD16 產(chǎn)品略低估8 天ET（bias：?1.79 mm/8 d），PML?V2產(chǎn)品精度（RMSE：8.85 mm/8 d）與MOD16 產(chǎn)品精度（RMSE：8.54 mm/8 d）相差不多；12個站點上PML?V2產(chǎn)品精度優(yōu)于MOD16產(chǎn)品精度，另外7個站點上MOD16 產(chǎn)品精度優(yōu)于PML?V2 產(chǎn)品精度；在BE?Lon、DE?Geb、DE?Kli、DK?Fou、FI?Jok和IT?Ca2 6 個站點，MOD16 產(chǎn)品（bias：0.67—15.05 mm/8 d）和PML?V2 產(chǎn) 品（bias：1.53—13.27 mm/8 d）均高估8天ET；在CH?Oe2、DE?Seh、US?CRT、US?Tw2、US?Tw3 和US?Twt 等6 個站點上，MOD16 產(chǎn)品（bias：?16.42—?1.04 mm/8 d）和PML?V2 產(chǎn)品（bias：?15.98—?2.28 mm/8 d）均低估8天ET。

圖2 19個FLUXNET農(nóng)田站點MOD16產(chǎn)品和PML?V2產(chǎn)品8天ET驗證散點圖（95%置信區(qū)間）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農(nóng)田站點MOD16產(chǎn)品和PML?V2產(chǎn)品與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 個像元內(nèi)MODIS 地表覆被類型完全一致均為農(nóng)田的10 個站點：在BE?Lon、DE?Geb、DE?Kli、US?CRT、US?Ne1、US?Ne2 和US?Ne3 等7 個站點上，PML?V2 產(chǎn)品精度（RMSE：4.42—7.23 mm/8 d）優(yōu)于或略優(yōu)于MOD16產(chǎn)品精度（RMSE：7.22—8.33 mm/8 d）；在FR?Gri、IT?Bci 和US?ARM 3 個站點上，MOD16 產(chǎn)品精度（RMSE：4.24—6.07 mm/8 d）優(yōu)于或略優(yōu)于PML?V2產(chǎn)品精度（RMSE：5.22—11 mm/8 d）。對于其站點所在像元MODIS 地表覆被類型不是農(nóng)田的3 個站點：站點所在像元為建筑用地的DE?RuS站點，MOD16 產(chǎn)品沒有產(chǎn)品估算結(jié)果，故在圖4和圖5 中未予以評價；所在像元為稀樹草原的FI?Jok 站點上，MOD16 產(chǎn)品精度較差RMSE 大于10 mm/8 d，PML?V2 產(chǎn)品RMSE 仍在5 mm/8 d 內(nèi)，兩種產(chǎn)品根據(jù)MODIS 地類分類計算，計算ET 時按稀樹草原的參數(shù)計算，可能是造成估算誤差大的原因；所在像元為草地的US?Tw2 站點，盡管所處位置非?？拷路睫r(nóng)田像元，但產(chǎn)品值應該是按草地的參數(shù)計算的草地ET 而非農(nóng)田ET，兩種產(chǎn)品 均低估8 天ET，PML?V2 產(chǎn)品精度（RMSE：4.45 mm/8 d，bias：?2.72 mm/8 d）優(yōu)于MOD16 產(chǎn)品精度（RMSE：8.06 mm/8 d，bias：?5.35 mm/8 d）。對于雖然站點所在像元MODIS 地表覆被類型是農(nóng)田，但是在站點附近3×3個像元內(nèi)含有非農(nóng)田類型像元的6 個站點：在CH?Oe2、DE?Seh 和IT?CA2等3 個站點上，兩種產(chǎn)品RMSE 均在10 mm/8 d 以內(nèi)，MOD16產(chǎn)品精度（RMSE：3.81—6.38 mm/8 d，bias：?2.78—0.67 mm/8 d）優(yōu)于或略優(yōu)于PML?V2產(chǎn)品精度（RMSE：5.75—9.58 mm/8 d，bias：?4.34—4.25 mm/8 d），在IT?CA2 站點上，MOD16 產(chǎn)品大 部分8 天累積ET 與觀測8 天ET 偏差在?10—10 mm/8 d，高估和低估的分布相對均勻，而PML?V2產(chǎn)品部分8天累積ET高估可達將近30 mm/8 d；在US?Tw3 和US?Lin兩個站點上，PML?V2 產(chǎn)品（RMSE：5.57—8.03 mm/8 d，bias：?5.39—0.66 mm/8 d）精度略優(yōu)于MOD16 產(chǎn)品精度（RMSE：6.81—11.03 mm/8 d，bias：?3.44—?8.55 mm/8 d）；地表異質(zhì)性較明顯的US?Twt 站點上，MOD16 產(chǎn)品（RMSE：21.47 mm/8 d，bias：?16.42 mm/8 d）和PML?V2產(chǎn)品（RMSE：22.4 mm/8 d，bias：?15.98 mm/8 d）精度均較差，嚴重低估8 天ET，該站點是唯一的水稻站點，證明了現(xiàn)有產(chǎn)品對水稻ET的嚴重低估。

圖4 中國8個農(nóng)田站點MOD16產(chǎn)品和PML?V2產(chǎn)品8天ET驗證（95%置信區(qū)間）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個農(nóng)田站點MOD16產(chǎn)品和PML?V2產(chǎn)品與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個農(nóng)田站點的8天ET 直接驗證結(jié)果，以及MOD16 產(chǎn)品和PML?V2 產(chǎn)品8天累積ET與EC觀測8天ET偏差的云雨圖。在LC、YC、GT和MY這4個站點上，PML?V2產(chǎn)品8天ET精度（R2：0.271—0.811，RMSE：3.3—9.2 mm/8 d）明顯優(yōu) 于MOD16 產(chǎn) 品8 天ET 精度（R2：0.018—0.600，RMSE：9.74—13.45 mm/8 d），PML?V2 產(chǎn)品在GT 站點高估8 天ET（bias：3.52 mm/8 d），在其余3個站點低估8天ET（bias：?1.72—?0.57 mm/8 d），MOD16產(chǎn)品在4個站點均低估8天ET（bias：?8.4—?6.21 mm/8 d）；在CN?Du1站點，兩種產(chǎn)品8天ET驗證精度相差不大且整體精度較其他站點相對較高，MOD16產(chǎn)品略低估8天ET（bias：?2.49 mm/8 d，RMSE：5.89 mm/8 d），而PML?V2 產(chǎn)品略高估8 天ET（bias：1.08 mm/8 d，RMSE：6.76 mm/8 d）；在CW 和SX兩個站點，MOD16產(chǎn)品8天ET精度（R2：0.143—0.439，RMSE：5.73—10.12 mm/8 d，bias：?2.29—2.37 mm/8 d）優(yōu)于PML?V2產(chǎn)品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 站點有效數(shù)據(jù)太少，統(tǒng)計意義不足，在僅有的3 個8 天，PML?V2產(chǎn)品略優(yōu)于MOD16產(chǎn)品。

圖6展示了中國8個農(nóng)田站點和全球28個農(nóng)田站點上MOD16 產(chǎn)品和PML?V2 產(chǎn)品8 天ET 總體驗證結(jié)果。對于中國8 個農(nóng)田站點，PML?V2 產(chǎn)品精度（R2：0.480，RMSE：8.46 mm/8 d）優(yōu)于MOD16產(chǎn)品精 度（R2：0.407，RMSE：10.77 mm/8 d）；MOD16 產(chǎn)品在7 個站點上均低估8 天ET，僅在SX站點上高估8天ET，整體低估8 天ET（bias：?6.29 mm/8 d）；而PML?V2 產(chǎn)品在4 個站點上高估8天ET，在另外4個站點上低估8天ET，整體略低估8 天ET（bias：?0.66 mm/8 d）。對于本研究中評估的全球共28個農(nóng)田站點，MOD16產(chǎn)品精度（R2：0.452，RMSE：8.82 mm/8 d）和PML?V2 產(chǎn)品精度（R2：0.455，RMSE：8.81 mm/8 d）整體精度相當，MOD16 產(chǎn)品整體略低估8 天ET（bias：?2.31 mm/8 d），而PML?V2 產(chǎn)品整體略高估8 天ET（bias：0.51 mm/8 d）。

圖6 MOD16產(chǎn)品和PML?V2產(chǎn)品8天ET整體驗證結(jié)果Fig.6 The overall validation of MOD16 ET and PML?V2 ET for 8?day estimates

3.2 趨勢分析檢驗

圖7展示了20個FLUXNET 農(nóng)田站點上MOD16產(chǎn)品和PML?V2 產(chǎn)品8 天累積ET 與EC 觀測8 天累積ET 的時序變化。在8 天中EC 觀測ET 其QC 顯示觀測和高質(zhì)量插補數(shù)據(jù)不足80%的在圖中用紅色*表示。該圖7 中R2、RMSE 與圖2 中R2和圖3 中RMSE的差別在于，該圖7中對MOD16產(chǎn)品有效數(shù)據(jù)和PML?V2有效數(shù)據(jù)分別計算R2和RMSE，而圖2和圖3 中只展示了MOD16 產(chǎn)品和PML?V2 產(chǎn)品數(shù)據(jù)值均有效的8 天，其計算R2和RMSE 相應的8 天的天數(shù)兩種產(chǎn)品是一致的。整體而言，PML?V2產(chǎn)品時間連續(xù)性優(yōu)于MOD16 產(chǎn)品，相比于考慮兩種產(chǎn)品有效數(shù)據(jù)時間一致的情況（圖2和圖3），大多數(shù)站點上PML?V2產(chǎn)品精度有所提升，在15個站點上，PML?V2產(chǎn)品精度優(yōu)于MOD16產(chǎn)品精度，但是MOD16產(chǎn)品可以相對更好的捕捉到年中因為作物輪種而引起的ET下降和上升趨勢，而PML?V2產(chǎn)品普遍在年中達到當年最高峰值。對于BE?Lon、DE?Geb和DE?Kli等3個站點，PML?V2產(chǎn)品與EC觀測ET二者曲線的吻合度相對更好（R2：0.78—0.8），而MOD16產(chǎn)品（R2：0.54—0.62）幾乎在每年的年中波峰處存在嚴重高估ET 的現(xiàn)象；對于CE?Oe2 站點，MOD16產(chǎn)品和PML?V2產(chǎn)品幾乎在每年ET高值處存在低估ET 的現(xiàn)象；在DK?Fou 站點，MOD16 和PML?V2 兩種產(chǎn)品時間變化特點相似（二者R2為0.85，趨勢均為先上升，在6月下旬至7月上旬達到波峰，而后下降），相比EC觀測ET（二者與觀測R2分別為0.28和0.31，趨勢先上升，在8月下旬至9月上旬達到波峰，而后下降）而言均高估ET（bias>10 W/m2）；在US?Lin站點，兩種產(chǎn)品均未捕捉到觀測ET 峰值，相比觀測ET 提前達到各自峰值，且PML?V2產(chǎn)品ET整體高于MOD16產(chǎn)品ET；US?Ne1、US?Ne2和US?Ne3等3個以玉米作物為主的站點，雖然相比于觀測8 天ET，PML?V2 產(chǎn)品整體精度高于MOD16產(chǎn)品，但MOD16產(chǎn)品8天ET對于細節(jié)的捕捉能力略優(yōu)于PML?V2產(chǎn)品，在上半年MOD16產(chǎn)品與EC觀測ET二者曲線的吻合度更好，而PML?V2產(chǎn)品明顯高估EC觀測8天ET。以US?Ne1站點為例，在3—6月，從與EC觀測ET時序曲線的吻合度上來看，MOD16產(chǎn)品（R2：0.70，bias：?4.02 mm/8 d，RMSE：49.02 mm/8 d）略優(yōu)于PML?V2產(chǎn)品（R2：0.64，bias：?4.53 mm/8 d，RMSE：58.82 mm/8 d）；在IT?Bci、IT?CA2 和US?ARM 等3 個站點，MOD16 產(chǎn)品與EC觀測8天ET變化趨勢更為相似，PML?V2產(chǎn)品高估8天ET，而MOD16產(chǎn)品可以相對較好的捕捉到年中的下降或上升趨勢，如在IT?CA2站點上，MOD16產(chǎn)品在年中更好地反映了ET 先下降后上升的趨勢；US?Tw2、US?Tw3 和US?Twt 等3 個站點上，兩種產(chǎn)品均低估8 天ET 表現(xiàn)不佳，但是MOD16 產(chǎn)品仍能相對更好捕捉到年中ET 的波動（下降或上升趨勢），US?Twt 站點上MOD16 產(chǎn)品對最高峰值所處時間的刻畫明顯優(yōu)于PML?V2產(chǎn)品，MOD16產(chǎn)品和觀測值均在下半年達到年度最高峰值，而PML?V2產(chǎn)品則提前在上半年達到年度最高峰值。

圖7 20個FLUXNET農(nóng)田站點MOD16產(chǎn)品、PML?V2產(chǎn)品和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個農(nóng)田站點上MOD16產(chǎn)品和PML?V2 產(chǎn)品8 天累積ET 與EC 觀測8 天累積ET 的時序變化。整體來看，PML?V2產(chǎn)品時間連續(xù)性優(yōu)于MOD16產(chǎn)品，PML?V2產(chǎn)品8天ET與EC觀測8天ET 整體吻合性更好（R2更高、RMSE 更低），而MOD16產(chǎn)品對部分站點時序細節(jié)的刻畫優(yōu)于PML?V2 產(chǎn)品。例如MOD16 產(chǎn)品可以捕捉到Y(jié)C 站點和LC 站點每年度年中的低谷，刻畫年中8 天ET 的先下降后上升趨勢，而在年中EC 觀測8天ET 處于低谷時PML?V2產(chǎn)品仍處于其年度峰值狀態(tài)。

圖8 中國8個農(nóng)田站點MOD16產(chǎn)品、PML?V2產(chǎn)品和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

為了更具體展示兩種產(chǎn)品對于時序曲線細節(jié)的刻畫，本研究選取IT?CA2、US?ARM、US?NE1、US?Twt、LC 和YC 等6 個站點，對MOD16產(chǎn)品和PML?V2 產(chǎn)品8 天累積ET 與EC 觀測8 天累積ET 開展了時序變化分解，分析其季節(jié)變化趨勢（圖9）。本研究采用基于貝葉斯集成的BEAST時序分解算法（Zhao 等，2019）來對ET 時序進行季節(jié)變化分解，相比于其他時序分解算法，該算法允許空缺值的存在，在MOD16 產(chǎn)品時序不完整的情況下，不需要進行插值處理就能獲得較為理想的結(jié)果。此處主要是為了對ET 時序變化的波峰和波谷所處時刻進行探究，故只展示了分解后的季節(jié)性波動信號變化。從圖9 中可以清晰看出，雖然MOD16 產(chǎn)品連續(xù)性較差，但是仍能較好的反映出年中ET的下降和上升趨勢。在IT?CA2和US?ARM兩個站點上，MOD16 產(chǎn)品與EC 觀測ET 季節(jié)變化波動曲線貼合度更好（IT?CA2 和US?ARM 的R2分別為0.83 和0.76），均在上半年（4—5 月）達到波峰，而PML?V2 產(chǎn)品（IT?CA2 和US?ARM 的R2分別為0.43 和0.69）則在年中（6—7 月）達到波峰；在US?Ne1 站點上，由于MOD16 產(chǎn)品數(shù)據(jù)在年初和年末的缺失，使得其沒有較好的捕捉到兩年交替之際ET 時序曲線的低谷，從與EC 觀測ET 時序曲線的吻合度上來看，MOD16產(chǎn)品（R2：0.92）整體略優(yōu)于PML?V2產(chǎn)品（R2：0.88），在上半年3—6月，MOD16 產(chǎn)品（R2：0.91）優(yōu)于PML?V2 產(chǎn)品（R2：0.82），MOD16 產(chǎn)品和EC 觀測ET 約從5—6 月開始急劇上升，上升曲線的坡度更陡（MOD16產(chǎn)品在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產(chǎn)品ET約從3月開始穩(wěn)步上升，上升曲線的坡度相對緩和（平均坡度一直維持在0.23—0.26范圍內(nèi)）；在US?Twt站點上，MOD16產(chǎn)品能夠捕捉到EC觀測ET的兩處小高峰，并且在下半年達到其最高峰，盡管PML?V2產(chǎn)品在2012年—2014年也能夠捕捉到兩個波峰，但是在上半年達到最高峰，與EC觀測ET和MOD16 產(chǎn)品相反；在LC 和YC 兩個站點上，PML?V2產(chǎn)品均在年中達到年度峰值，未捕捉到冬小麥、夏玉米輪種而導致的年中ET先下降再上升的趨勢，而MOD16產(chǎn)品成功捕捉到冬小麥和夏玉米生長期內(nèi)的兩個ET波峰，并在下半年達到最高峰（夏玉米生長季），與EC觀測ET恰好相反，EC觀測ET在上半年達到最高峰（冬小麥生長季），兩個站點上MOD16產(chǎn)品均低估了冬小麥ET。

圖9 選取的6個農(nóng)田站點EC觀測8天ET、MOD16產(chǎn)品和PML?V2產(chǎn)品8天ET時序分解季節(jié)波動信號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 結(jié)論

本文采用20 個FLUXNET 通量觀測站點和8 個中國的通量觀測站點共計28 個農(nóng)田站點評估了當前國際上僅有的兩種500 m 分辨率的8 天ET 產(chǎn)品：MOD16產(chǎn)品和PML?V2產(chǎn)品。整體而言，28個農(nóng)田站點的驗證表明MOD16產(chǎn)品（RMSE：8.82 mm/8 d）和PML?V2產(chǎn)品（RMSE：8.81 mm/8 d）總體精度相當；MOD16產(chǎn)品低估8天ET（bias：?2.31 mm/8 d），PML?V2 產(chǎn)品略高估8 天ET（bias：0.51 mm/8 d）；PML?V2 產(chǎn)品時間連續(xù)性優(yōu)于MOD16 產(chǎn)品，但是部分站點MOD16 產(chǎn)品在時序變化細節(jié)（如到達峰值的季節(jié)、年中的下降上升趨勢）的刻畫上優(yōu)于PML?V2產(chǎn)品。20個FLUXNET通量觀測站點評估結(jié)果表明，兩種產(chǎn)品精度相當，PML?V2產(chǎn)品（RMSE：4.42—22.4 mm/8 d，bias：?15.98—13.27 mm/8 d）略低估8天ET，MOD16產(chǎn)品（RMSE：3.81—21.48 mm/8 d，bias：?16.42—15.05 mm/8 d）略高估8 天ET；PML?V2產(chǎn)品時間連續(xù)性較好，在不考慮兩種產(chǎn)品驗證點數(shù)一致的情況下，PML?V2 產(chǎn)品精度（RMSE：4.06—20.36 mm/8 d）有所提升。PML?V2算法利用FLUXNET站點觀測通量對農(nóng)田地類進行了單獨參數(shù)率定，采用留一法在11 個農(nóng)田站點上進行了算法精度驗證，結(jié)果表明對于農(nóng)田地類該算法估算RMSE 變化范圍為0.54—2.38 mm/d，R2變化范圍為0.42—0.86，MOD16 產(chǎn)品RMSE 變化范圍為0.48—2.4 mm/d，R2變化范圍為0.49—0.76（Zhang 等，2019）。本研究與該研究中對農(nóng)田地類的8 天ET 估算精度驗證結(jié)果相當，均顯示US?Twt站點上水稻ET 的估算結(jié)果最差。本文采用的8 個中國農(nóng)田通量站點并未用于PML?V2 算法參數(shù)率定，但中國站點的驗證結(jié)果表明兩種產(chǎn)品均低估8天ET，PML?V2 產(chǎn)品（RMSE：2.88—11.25 mm/8 d，bias：?1.72—8.04 mm/8 d）仍具有相對可靠的精度，且精度優(yōu)于MOD16 產(chǎn)品（RMSE：5.73—13.45 mm/8 d，bias：?8.4—2.37 mm/8 d），在不考慮兩種產(chǎn)品驗證點數(shù)一致的情況下，PML?V2產(chǎn)品精度（R2：0.28—0.87，RMSE：3.23—9.58 mm/8 d）有所提升，但MOD16 產(chǎn)品對年中8 天ET 的下降和上升趨勢刻畫優(yōu)于PML?V2產(chǎn)品。

本文的研究仍存在不足，雖然結(jié)合MODIS 地表覆被產(chǎn)品和Google Earth 衛(wèi)星遙感影像對這些站點附近的地類和地表均一性進行了判斷，但是仍顯粗糙，未來建議結(jié)合高空間分辨率的植被指數(shù)、葉面積指數(shù)等植被生態(tài)參數(shù)來進行同時期更精細的地表均一性評價。由于地表的空間異質(zhì)性，地面觀測和遙感估算的蒸散發(fā)之間存在著空間尺度不匹配的問題，如何從地面站點觀測中獲取遙感像元尺度地表蒸散發(fā)的相對真值，是目前蒸散發(fā)研究的熱點和難點。對于站點像元空間異質(zhì)性強的情況，需要考慮觀測源區(qū)，且可能要進行空間尺度轉(zhuǎn)換，但在實際操作中并不容易；站點所處500 m 像元的空間代表性或均一性較好，可以直接進行驗證。已有研究報道，館陶農(nóng)田站點觀測空間范圍大約在450—460 m（Liu 等，2011），多倫農(nóng)田站點的觀測空間范圍大約在300—400 m，壽縣農(nóng)田站點的觀測空間范圍大約在200—400 m（王有恒，2010）。出于對農(nóng)田站點觀測源區(qū)一般小于500 m 的假設，本研究只選取了500 m 空間分辨率的MOD16 產(chǎn)品和PML?V2 產(chǎn)品，未來應考慮結(jié)合其他時空分辨率產(chǎn)品進行交叉檢驗，如1 km空間分辨率的GLASS ET 產(chǎn)品、BESS ET 產(chǎn)品和ETMonitor ET 產(chǎn)品，并結(jié)合高時間分辨率的ET 產(chǎn)品如GLEAM 產(chǎn)品，進行更大尺度像元的時空趨勢驗證。另外，站點觀測通量的30 min 到日尺度、8 天尺度、月尺度和年尺度的時間尺度轉(zhuǎn)換方法亦是目前研究的難點。500 m 尺度像元上的地表空間相對均一，可以考慮直接簡單平均，而對于異質(zhì)地表，可以考慮通量觀測源區(qū)的時間和空間加權(quán)來獲取8 天地面觀測ET，但是在實際操作時并非易事。這些難點需要在未來的研究中予以更好的考慮，探索更合理的空間尺度和時間尺度轉(zhuǎn)化方法，以獲取更可靠的蒸散發(fā)像元尺度驗證結(jié)果。

志 謝感謝FLUXNET 全球通量觀測研究網(wǎng)絡、ChinaFLUX 中國通量觀測研究網(wǎng)絡、中國生態(tài)系統(tǒng)研究網(wǎng)絡、國家生態(tài)科學數(shù)據(jù)中心、國家青藏高原數(shù)據(jù)中心、臺灣農(nóng)業(yè)研究所等研究網(wǎng)絡和機構(gòu)提供寶貴的農(nóng)田站點通量觀測數(shù)據(jù)，感謝中國科學院地理科學與資源研究所于貴瑞研究員、北京師范大學劉紹民教授、中國科學院青藏高原研究所李新研究員等科研學者共享數(shù)據(jù)及在推動觀測網(wǎng)絡建設和數(shù)據(jù)共享上作出的貢獻，同時感謝各個農(nóng)田站點工作人員在數(shù)據(jù)觀測整理過程中付出的辛勤勞動，在此表示衷心的感謝！