魏 濤，王云權(quán)

（中國地質(zhì)大學(xué)（武漢）環(huán)境學(xué)院，湖北 武漢 430074）

蒸散發(fā)的研究對于水資源的合理配置、農(nóng)作物的增產(chǎn)增量以及區(qū)域氣候的動態(tài)分析意義重大［1］。目前對于蒸散發(fā)的直接測量一般是集中在小尺度范圍，測量方法有蒸滲儀、波文比、渦度相關(guān)等，對于大區(qū)域尺度的蒸散發(fā)一般無法通過實(shí)際測得［2］。隨著遙感技術(shù)和經(jīng)驗(yàn)理論模型的發(fā)展，使得大區(qū)域尺度的實(shí)際蒸散發(fā)估算成為了可能。目前全球尺度與區(qū)域擬合效果較好的產(chǎn)品主要有基于PM（Penman-Monteith）公式的MOD16［3］、PML-V2［4］和基于PT（Priestley-Taylor）公式的GLEAM［5］蒸散發(fā)產(chǎn)品。

蒸散發(fā)的物理機(jī)制受多種氣候環(huán)境等地表參量的影響，通過對蒸散發(fā)在全國時空動態(tài)特征的分析，能夠從宏觀角度對全國地域氣候蒸散情況和未來蒸散發(fā)的變化情況作直觀的了解，為干旱預(yù)警、水資源調(diào)度提供數(shù)據(jù)來源和分析依據(jù)。2011 年劉昌明等［6］利用優(yōu)化后的PM 公式計(jì)算了中國的潛在蒸散發(fā)，并分析了相關(guān)氣象環(huán)境因子對蒸散發(fā)變化趨勢的影響；張巧鳳等［7］基于MODIS 蒸散發(fā)數(shù)據(jù)對錫林郭勒草原近14 a的蒸散發(fā)時空動態(tài)作分析，發(fā)現(xiàn)錫林郭勒草原蒸散發(fā)總體表現(xiàn)為由東北部和東南延線向西南部遞減且空間變異比較穩(wěn)定；鄧興耀等［8］利用MODIS 蒸散發(fā)數(shù)據(jù)集對天山山區(qū)2000 年到2014年蒸散發(fā)作了時空動態(tài)特征分析，發(fā)現(xiàn)蒸散發(fā)在空間上有西部大東部小、北部大南部小，且近15 a 全區(qū)蒸散發(fā)變異程度不明顯；鄭徳鳳等［9］利用河北省12 個氣象站的實(shí)際氣象觀測資料，采用PM公式計(jì)算了河北省1967—2019年的潛在蒸散量，分析發(fā)現(xiàn)研究時段內(nèi)蒸散發(fā)存在明顯的地帶性特征，其總體呈上升趨勢。對于蒸散發(fā)的時空動態(tài)分析，大多數(shù)學(xué)者主要集中在區(qū)域性尺度的蒸散發(fā)分析，或通過改進(jìn)PM 公式計(jì)算的潛在蒸散發(fā)進(jìn)行分析，或選取蒸散發(fā)產(chǎn)品分析時沒有考慮到蒸散發(fā)產(chǎn)品在區(qū)域的適宜性。本文通過全球通量站網(wǎng)發(fā)布的FLUXNET2015數(shù)據(jù)集中中國的9個通量站數(shù)據(jù)［10］，對上述3種蒸散發(fā)產(chǎn)品在中國區(qū)域的適用性進(jìn)行了評價。在此基礎(chǔ)上，利用效果最優(yōu)的蒸散發(fā)產(chǎn)品對中國蒸散發(fā)的時空動態(tài)特征進(jìn)行了分析，為中國的水資源規(guī)劃及生態(tài)水文研究提供參考及借鑒意義。

1 研究區(qū)概況

中國地處亞洲東南部，總體地勢西高東低，氣候類型多種多樣。根據(jù)中國科學(xué)院地理科學(xué)與資源研究所資源環(huán)境科學(xué)與數(shù)據(jù)中心氣候區(qū)劃數(shù)據(jù)，以日平均氣溫不低于10 ℃穩(wěn)定期的積溫和最冷月氣溫或極端最低氣溫多年平均值為熱量指標(biāo)，以干燥度（見氣候指數(shù)）為水分指標(biāo)將研究區(qū)劃分為10個一級氣候區(qū)進(jìn)行分區(qū)討論［11］（圖1），其中瓊南西沙區(qū)和南沙區(qū)各島嶼分別屬于中熱帶（VIII）和南熱帶（IX），島嶼面積較小且分散，對于總體研究意義不大，故本次研究不對2個區(qū)域的蒸散發(fā)作探討。

圖1 研究區(qū)示意圖Fig.1 Sketch map of the study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

2.1.1 通量站數(shù)據(jù)由于FLUXNET2015（https://fluxnet.org）數(shù)據(jù)集的完善性與缺失值或異常值處理方法的先進(jìn)性，本次研究選取數(shù)據(jù)集中位于中國區(qū)域的9 個共享數(shù)據(jù)站點(diǎn)（剔除了一個誤差較大的數(shù)據(jù)站點(diǎn)），對全球蒸散發(fā)產(chǎn)品在中國地區(qū)的適宜性進(jìn)行研究，站點(diǎn)分別是CN-Cha、CN-Cng、CN-Dan、CN-Din、CN-Du2、CN-Du3、CN-Ha2、CN-HaM、CNQia，站點(diǎn)分布如圖1，站點(diǎn)詳細(xì)信息如表1。

表1 通量站點(diǎn)基本信息Tab.1 Basic information of flux site

2.1.2 GLEAM蒸散發(fā)產(chǎn)品數(shù)據(jù)GLEAM（Global land evaporation amsterdam model）蒸散發(fā)產(chǎn)品數(shù)據(jù)集產(chǎn)品通過衛(wèi)星遙感數(shù)據(jù)估算0.25°分辨率每天的實(shí)際地表蒸散發(fā)和植物根系土壤水含量。該數(shù)據(jù)集算法由Miralles 于2011 年首次發(fā)布［12］，2014 年改進(jìn) 為GLEAM（v2）［13］，2017 修 訂 發(fā) 展 為 最 新 的GLEAM（v3）系列產(chǎn)品。GLEAM（v3）在保持水分通量與之前版本大致一致的情況下改進(jìn)了與91 個通量站點(diǎn)和2325 個土壤水分傳感器各層土壤水分含量的相關(guān)系數(shù)。數(shù)據(jù)下載網(wǎng)址www.GLEAM.eu。為了對產(chǎn)品進(jìn)行中國區(qū)域適宜性分析，本次研究根據(jù)通量站數(shù)據(jù)及產(chǎn)品數(shù)據(jù)時間范圍下載了2002 年1月—2010年12月GLEAM（v3）每月的蒸散發(fā)數(shù)據(jù)。

2.1.3 MOD16數(shù)據(jù)MOD16 數(shù)據(jù)集是由美國國家航空局通過Aqua 和Terra 衛(wèi)星平臺數(shù)據(jù)通過MOD16算法生產(chǎn)的全球地表蒸散發(fā)產(chǎn)品數(shù)據(jù)集，該數(shù)據(jù)集產(chǎn)品包含空間尺度500 m 和時間尺度8 d 的地表蒸散發(fā)，覆蓋面積為全球1.09×108km2的陸地，涵蓋了各種植被類型，該數(shù)據(jù)集的有效數(shù)據(jù)日期為2000年1月至今［14］。MOD16算法基于PM公式計(jì)算的蒸散發(fā)主要分3部分：土壤蒸散發(fā)、濕潤的冠層蒸散發(fā)和植物干冠層表面的蒸騰，通過計(jì)算各部分蒸散量的基礎(chǔ)上分別計(jì)算白天和夜間蒸散發(fā)并整合到天尺度。該產(chǎn)品數(shù)據(jù)通過AmeriFlux的46個通量站點(diǎn)數(shù)據(jù)與產(chǎn)品數(shù)據(jù)集對比顯示，改進(jìn)的算法有效地將平均絕對誤差（MAE）由0.39 mm·d-1提高到0.33 mm·d-1。為了對產(chǎn)品進(jìn)行中國區(qū)域適宜性分析，本次研究根據(jù)通量站數(shù)據(jù)及產(chǎn)品數(shù)據(jù)時間范圍下載了2002年1月—2010年12月每8 d的蒸散發(fā)數(shù)據(jù)，并通過取月均值的方法，將時間尺度處理為每月。數(shù)據(jù)通過GEE（Google Earth Engine）平臺下載。

2.1.4 PML-V2數(shù)據(jù)PML-V2數(shù)據(jù)集是中國科學(xué)院地理科學(xué)與資源研究所Zhang等［4］2019年利用谷歌地球引擎GEE 平臺，以遙感數(shù)據(jù)為驅(qū)動，基于Penman-Monteith-Leuning 蒸散發(fā)模型計(jì)算的全球地表蒸散發(fā)數(shù)據(jù)，該數(shù)據(jù)集的空間分辨率為500 m，時間分辨率為8 d［4］。該數(shù)據(jù)集通過全球91 個廣泛分布的通量站點(diǎn)（涵蓋了10種植被類型）數(shù)據(jù)驗(yàn)證，產(chǎn)品數(shù)據(jù)集的蒸散發(fā)均方根差和偏差有效的減小到0.69 mm·d-1和-1.8%。PML-V2 數(shù)據(jù)集時間范圍為2002年7月—2020年12月，本次選取該時間段的全部數(shù)據(jù)，并通過取月均值的方法，獲取2002年7月—2020年12月的月尺度數(shù)據(jù)，用以適應(yīng)性分析對比。

2.1.5 其他參考數(shù)據(jù)中國區(qū)域的DEM高程數(shù)據(jù)為美國國家航空航天局噴氣推進(jìn)實(shí)驗(yàn)室SRTM v3 產(chǎn)品［15］，分辨率為30 m，通過GEE平臺下載。

2.2 研究方法

2.2.1 蒸散發(fā)產(chǎn)品區(qū)域適宜性評價為了分析GLEAM、MOD16 和PML-V2 3 種蒸散發(fā)產(chǎn)品在中國區(qū)域的適宜性，本研究通過FLUXNET2015 數(shù)據(jù)集中中國區(qū)域通量站點(diǎn)的通量站蒸散發(fā)數(shù)據(jù)與3種蒸散發(fā)產(chǎn)品的站點(diǎn)蒸散發(fā)數(shù)據(jù)作對比，利用數(shù)據(jù)擬合的皮爾遜相關(guān)系數(shù)、均方根差分析3 種蒸散發(fā)產(chǎn)品在中國的適宜性［16-17］。

（1）皮爾遜相關(guān)系數(shù)。皮爾遜相關(guān)系數(shù)廣泛用于度量2 個變量之間的相關(guān)程度，本研究通過對比站點(diǎn)實(shí)測數(shù)據(jù)與蒸散發(fā)產(chǎn)品數(shù)據(jù)的相關(guān)性強(qiáng)度作為蒸散發(fā)產(chǎn)品在中國區(qū)域適宜性的評價指標(biāo)，計(jì)算公式如下：

式中：r為皮爾遜相關(guān)系數(shù)；xi、yi分別為樣本值；、分別為兩組數(shù)據(jù)的均值；n為樣本總數(shù)。

（2）均方根差。均方根差（RMSE）作為一種常用測量數(shù)值之間差異的量度，本研究通過站點(diǎn)觀測值與蒸散發(fā)產(chǎn)品模型值的RMSE來評價蒸散發(fā)產(chǎn)品模擬中國區(qū)域蒸散發(fā)實(shí)際值，計(jì)算公式如下：

式中：xobs,i為第i個觀測數(shù)據(jù)；xmod,i為第i個模擬數(shù)據(jù)。

2.2.2 時空動態(tài)特征分析蒸散發(fā)時空動態(tài)特征分析主要是利用各種分析方法，對蒸散發(fā)在時間和空間上的動態(tài)變化進(jìn)行分析及預(yù)測，發(fā)現(xiàn)其中的規(guī)律，為更好規(guī)劃和進(jìn)行水資源配置提供參考依據(jù)［18-19］，分析方法如下表2［20-26］。

表2 分析方法Tab.2 Analytical method

3 結(jié)果與分析

3.1 蒸散發(fā)產(chǎn)品適宜性評價

本文9個通量站點(diǎn)蒸散發(fā)擬合情況如圖2、圖3所示：

圖3 站點(diǎn)擬合Taylor圖Fig.3 Taylor diagram of site fitting

由圖2 可知，除了CN-Du3 站點(diǎn)外，各通量站點(diǎn)數(shù)據(jù)與GLEAM、MOD16和PML-V2蒸散發(fā)產(chǎn)品數(shù)據(jù)都有較好相關(guān)性，蒸散發(fā)隨時間變化的趨勢大體一致。由圖3 的Taylor 圖可以觀察相關(guān)系數(shù)、標(biāo)準(zhǔn)差和均方根差3 個指標(biāo)的綜合效果，GLEAM 與PMLV2 與觀測值距離較近，PML-V2 除CN-Du3 站點(diǎn)以外，3 個指標(biāo)均不如CLEAM 和MOD16 產(chǎn)品。綜上，PML-V2各項(xiàng)評價指標(biāo)相對其他兩種產(chǎn)品均表現(xiàn)出一定優(yōu)勢，在時間和空間尺度上與站點(diǎn)蒸散發(fā)數(shù)據(jù)的擬合效果較好，趨勢一致，在中國區(qū)域有較好的適宜性。因此，對中國區(qū)域的蒸散發(fā)時空動態(tài)特征分析選取了PML-V2產(chǎn)品數(shù)據(jù)集。

3.2 2003—2020年中國蒸散發(fā)空間分布特征

將PML-V2 蒸散發(fā)產(chǎn)品8 d×500 m 分辨率數(shù)據(jù)匯總計(jì)算全國多年平均蒸散發(fā)。將多年平均蒸散發(fā)分為3級（圖4、圖5）：低蒸散發(fā)區(qū)（ET＜200 mm·a-1）、中蒸散發(fā)區(qū)（200 mm·a-1≤ET≤400 mm·a-1）、高蒸散發(fā)區(qū)（ET＞400 mm·a-1）［8］。蒸散發(fā)全國分布由西北向南呈逐漸增加趨勢，區(qū)域內(nèi)最大值一般出現(xiàn)在水域附近。由變異系數(shù)可知，北溫帶區(qū)（I）、中亞熱帶區(qū)（V）、南亞熱帶區(qū)（VI）、北熱帶區(qū)（VII）、中熱帶區(qū)（VIII）空間蒸散發(fā)變異程度較低，蒸散發(fā)區(qū)域變化不大。中溫帶區(qū)（II）、南溫帶區(qū)（III）、北亞熱帶區(qū)（IV）、高原氣候區(qū)（H）變異程度較大，蒸散發(fā)在空間分布相差較大（表3）。貧瘠地區(qū)[至少60%的面積為植被不足10%的無植被貧瘠（沙、巖、土）地區(qū)]主要分布在中國西北地區(qū)，受氣候和覆蓋類型的共同作用蒸散發(fā)較低［27］。在同一緯度地區(qū)，不同土地覆蓋類型和所處的地理氣候類型區(qū)蒸散發(fā)也有不同程度的差異［28-29］。氣候分區(qū)圖與蒸散發(fā)多年均值有較好的一致性趨勢。各氣候分區(qū)內(nèi)蒸散發(fā)的最大值隨分區(qū)類型由西北像東南逐漸增加。

表3 2003—2020年中國蒸散發(fā)統(tǒng)計(jì)分析Tab.3 Statistical analysis of evapotranspiration in China from 2003 to 2020

圖4 2003—2020年中國平均蒸散發(fā)分級Fig.4 Average evapotranspiration classification of China from 2003 to 2020

圖5 2003—2020年中國多年平均蒸散發(fā)空間分布Fig.5 Spatial distribution of annual mean evapotranspiration in China from 2003 to 2020

3.3 2003—2020年中國蒸散發(fā)時序變異特征

由2003—2020 年中國蒸散發(fā)時序變異系數(shù)空間分布（圖6）和變異系數(shù)數(shù)據(jù)統(tǒng)計(jì)（表4）可以看出，北溫帶區(qū)（I）、中溫帶區(qū)（II）、南溫帶區(qū)（III）、高原氣候區(qū)（H）3 個區(qū)變異系數(shù)分別為0.37、0.35、0.38、0.34，在時間序列上的變異相對其他區(qū)較大。與上部分空間分析的變異性類似，蒸散發(fā)時空變異較大區(qū)域主要為中國西北部，在沿海地區(qū)城市與較高的時間變異系數(shù)相關(guān)性較強(qiáng)。變異系數(shù)較大區(qū)域大多集中在盆地邊緣山脈附近，濱海城市亦有較大的變異系數(shù)區(qū)。

表4 2003—2020年中國蒸散發(fā)時序變異系數(shù)統(tǒng)計(jì)分析Tab.4 Statistical analysis of temporal variation coefficients of evapotranspiration in China from 2003 to 2020

圖6 2003—2020年中國蒸散發(fā)時序變異系數(shù)空間分布Fig.6 Spatial distribution of temporal variation coefficients of evapotranspiration in China from 2003 to 2020

2003—2020 年中國全區(qū)蒸散發(fā)時間上的變化特征如圖7所示。并對各區(qū)每年平均蒸散發(fā)進(jìn)行一元線性回歸，擬合公式為：ET=B0+B1×X，其中X為年份，B0為擬合回歸常數(shù)項(xiàng)，B1為擬合回歸一次項(xiàng)系數(shù)，各擬合回歸系數(shù)如表5。

表5 擬合回歸系數(shù)Tab.5 Fit regression coefficient

圖7 年均蒸散發(fā)時間變化圖Fig.7 Diagram of annual evapotranspiration changing

從圖7 可以看出，2003—2020 年各區(qū)蒸散發(fā)有一定的波動。各區(qū)蒸散發(fā)波動范圍與上文空間分析趨勢一致，由西北向東南區(qū)逐漸增大。

全區(qū)蒸散發(fā)2003—2020年波動微弱，為增大趨勢，變化率為0.89 mm·a-1。為了研究蒸散發(fā)在時間上的變化趨勢及趨勢程度的顯著性和可信度，通過Theil-Sen median（T-S）和Mann-Kendall（M-K）分析方法，進(jìn)一步計(jì)算了2003—2020年的趨勢斜率（T-S分析系數(shù)如圖8）和歸一化檢驗(yàn)統(tǒng)計(jì)量（M-K 分析Z值如圖9）。

圖8 2003—2020年中國T-S分析系數(shù)空間分布Fig.8 Spatial distribution of T-S analytical index in China from 2003 to 2020

圖9 2003—2020年中國M-K分析Z值空間分布Fig.9 Spatial distribution of M-K analytical Z value in China from 2003 to 2020

對中國各區(qū)進(jìn)行區(qū)域統(tǒng)計(jì)分析如表6：

表6 中國不同分區(qū)T-S分析系數(shù)和M-K分析Z值統(tǒng)計(jì)Tab.6 Statistics of T-S analytical index and M-K analytical Z value of each zone in China

Theil-Sen median 趨勢分析與一元線性回歸分析基本一致，全區(qū)蒸散發(fā)程增長趨勢。由M-K分析Z值圖和相應(yīng)的區(qū)域統(tǒng)計(jì)分析表可知，除了北溫帶區(qū)（I）的置信度水平較高外，其他區(qū)域的置信度水平較低，即蒸散發(fā)增加或減少在全區(qū)的趨勢不顯著。M-K 趨勢顯著性分析時一般采取雙邊趨勢檢驗(yàn)，通常將置信度90%、95%和99%對應(yīng)Z絕對值1.65、1.96 和2.58 作為臨界值，將顯著性劃分為：極顯著（|Z|＞2.58）、顯著（1.96＜|Z|≤2.58）、微顯著（1.65＜|Z|≤1.96）、不顯著（0＜|Z|≤1.65）和無變化（Z=0）。在中國區(qū)域內(nèi)，2003—2020年蒸散發(fā)的增長或下降一般趨勢不顯著，但是在局部區(qū)域內(nèi)，蒸散發(fā)存在明顯的增長或下降趨勢。

3.4 蒸散發(fā)的預(yù)測

通過中國區(qū)域2003—2020年P(guān)ML-V2的蒸散發(fā)數(shù)據(jù)計(jì)算中國區(qū)域的赫斯特指數(shù)（圖10），對中國區(qū)域赫斯特指數(shù)進(jìn)行統(tǒng)計(jì)分析有：北溫帶區(qū)（I）赫斯特指數(shù)略微大于0.5，即區(qū)域未來蒸散發(fā)與過去趨勢一致，呈略微增長趨勢。南溫帶區(qū)（III）赫斯特指數(shù)為0.44，未來蒸散發(fā)與過去呈相反趨勢，即略微增長。北亞熱帶區(qū)（IV）、中亞熱帶區(qū)（V）、北熱帶區(qū)（VII）、中熱帶區(qū)（VIII）、高原氣候區(qū)（H）赫斯特指數(shù)均小于0.5，未來蒸散發(fā)與過去呈相反趨勢，預(yù)測未來蒸散發(fā)可能會減小。中溫帶區(qū)（II）和中亞熱帶區(qū)（V）由于兩種方法趨勢性分析結(jié)果不一致，趨勢變化率較小，且置信度水平較低，未來趨勢可能在穩(wěn)定水平波動。中國總體赫斯特指數(shù)小于0.5，即未來趨勢與過去相反。由上文，蒸散發(fā)的時間變化特征可知中國蒸散發(fā)為增長趨勢，未來蒸散發(fā)即可能與過去增長相反，逐漸減小。

圖10 2003—2020年中國赫斯特指數(shù)空間分布Fig.10 Spatial distribution of Hurst index in China from 2003 to 2020

4 結(jié)論

（1）中國多年平均蒸散發(fā)由西北向東南逐漸增大，主要受地理氣象要素影響。各分區(qū)蒸散發(fā)中中溫帶區(qū)（II）、南溫帶區(qū)（III）、北亞熱帶區(qū)（IV）和高原氣候區(qū)（H）空間變異系數(shù)較大，主要原因是分區(qū)空間上的異質(zhì)性較高，其他分區(qū)變異系數(shù)較低且空間地理環(huán)境比較穩(wěn)定。

（2）對中國蒸散發(fā)的時間特征分析顯示，在相對缺水干燥的北方地區(qū)，時間變異系數(shù)一般較大，而南方濕潤區(qū)相對系數(shù)則相對偏低。濱海城市和盆地邊緣山脈附近，時間變異系數(shù)較大，受不規(guī)則的局部氣候條件影響較大。同時，各分區(qū)內(nèi)也存在時間變異系數(shù)差別較大區(qū)域，主要與土地覆蓋類型及地理環(huán)境有關(guān)。各區(qū)蒸散發(fā)在長時序上呈小范圍波動變化，各區(qū)年蒸散發(fā)最大值與最小值出現(xiàn)時間并不完全一致，但是在相關(guān)鄰域一般最大值與最小值出現(xiàn)的時間較為相關(guān)。

（3）通過兩種方法得出北溫帶區(qū)（I）、北亞熱帶區(qū)（IV）、中亞熱帶區(qū)（V）、北熱帶區(qū)（VII）、中熱帶區(qū)（VIII）、高原氣候區(qū)（H）蒸散發(fā)均為增長趨勢，南溫帶區(qū)（III）為減小趨勢。中溫帶區(qū)（II）和南亞熱帶區(qū)（VI）兩種分析方法所得的變化趨勢不一致，且置信度水平較低，主要原因可能是兩區(qū)蒸散發(fā)相對穩(wěn)定，波動較小，增長或減小趨勢不明顯。結(jié)合時間變化趨勢和赫斯特指數(shù)對未來中國蒸散發(fā)進(jìn)行預(yù)測顯示，赫斯特指數(shù)除了在北溫帶區(qū)（I）大于0.5，即未來的變化狀況與過去趨勢一致，但赫斯特指數(shù)偏小，持續(xù)性較小。其他區(qū)域赫斯特均值均小于0.5，即未來趨勢與過去相反，但值與0.5接近，持續(xù)性也不強(qiáng)。本次研究發(fā)現(xiàn)，地理氣候因素可能是影響大尺度區(qū)域蒸散發(fā)的主要因素。

（4）研究尺度的選擇對于問題的研究息息相關(guān)，為了探究影響蒸散發(fā)在時間空間上的不均勻性及影響因子，需要通過不同尺度，不同分析數(shù)據(jù)進(jìn)行更深層次的研究。