韓再惠，陳 燕，吳志俊，任 黎，苗 平，達古拉，李瓊芳，5

(1.內蒙古自治區(qū)鄂爾多斯市杭錦旗水利事業(yè)發(fā)展中心，內蒙古 鄂爾多斯 017400；2.河海大學水文水資源學院，江蘇 南京 210024；3.內蒙古自治區(qū)鄂爾多斯市河湖保護中心，內蒙古 鄂爾多斯 017400；4.內蒙古自治區(qū)鄂爾多斯市水利事業(yè)發(fā)展中心，內蒙古 鄂爾多斯 017400；5.江蘇省南京市長江保護與綠色發(fā)展研究院，江蘇 南京 210024；)

1 概述

地表蒸散發(fā)是土壤-植被-大氣系統(tǒng)中能量與水分轉換的主要過程，也是能量平衡與水量平衡的關鍵組成因素。因此，準確估算地表蒸散發(fā)對區(qū)域水循環(huán)和水資源綜合利用方面具有重要意義[1]。然而由于空間異質性和觀測站點稀疏等原因，傳統(tǒng)基于站點數據估算蒸散發(fā)的方法存在較多問題。因此，國內外研究者們開發(fā)了不同類型的大尺度地表蒸散發(fā)估算方法，總共可分為5類：經驗統(tǒng)計公式、遙感正蒸散模型、陸面水文模型、基于蒸散發(fā)互補方法地表蒸散發(fā)模型和基于水量平衡的蒸散發(fā)估算方法，并依托這些方法陸續(xù)生產和發(fā)布了不少全球和區(qū)域尺度的蒸散發(fā)數據產品[2]。然而，Zhang等[3]發(fā)現大多數蒸散發(fā)產品在計算冠層氣孔阻力時，未能將蒸散發(fā)和植被初級總生產力(Gross Primary Productivity，GPP)耦合起來。因此，不能充分考慮不同環(huán)境因子對植被冠層氣孔導度的影響，在流域尺度上的地表蒸散發(fā)估算存在較大的不確定性。但遺憾的是，目前為止很少有研究在計算蒸散發(fā)的過程中引入GPP過程。

與此同時各類蒸散發(fā)產品由于模型結構、輸入數據集、模型參數存在差異等原因，都不可避免地有一定的不確定性。因此，迫切需要對其在蒸散發(fā)估算中的精度和不確定性進行分析。水量平衡法能夠在流域尺度對蒸散發(fā)估算值進行較為精準的驗證，因此被業(yè)內研究者廣泛使用[4]。TCH方法則可以在沒有蒸散發(fā)觀測值的情況下，量化流域蒸散發(fā)柵格尺度上的不確定性[5]。目前，在區(qū)域尺度上對蒸散發(fā)產品評價的研究較多，但均未對引入GPP過程的遙感蒸散發(fā)模型和傳統(tǒng)蒸散發(fā)產品在不同時空尺度上采用不同方法進行綜合對比研究。本研究通過引入GPP過程以實現對遙感蒸散發(fā)模型的改進，并采用水量平衡方法和TCH方法，綜合對比評價改進SW模型和不同蒸散發(fā)產品在黃河流域不同時空尺度下的精度和不確定性，以期為更加準確估算黃河流域蒸散發(fā)提供有益參考。

2 數據來源與研究方法

2.1 研究區(qū)概況

黃河流域位于中國的北部地區(qū)，面積為7.95×105km2(如圖1所示)。多年平均氣溫為7.2℃，多年平均降水量為495.6mm。由于地理位置，地形，植被覆蓋等因素影響，該流域全年降水分布極不均勻，約60%～80%的降雨量集中在6—9月。流域主要植被類型為常綠針葉林(ENF)、落葉針葉林(DNF)、常綠闊葉林(EBF)、落葉闊葉林(DBF)、灌叢、草地和農田[6]。本研究根據黃河流域自然地理環(huán)境特征，選擇花園口水文站作為黃河流域的徑流控制站。該水文站的控制流域面積為7.3×105km2，約占黃河流域面積的92%。

圖1 研究區(qū)地理空間分布

2.2 數據收集與處理

本研究采用了兩種被廣泛使用的蒸散發(fā)產品，改進SW模型所需的氣象和遙感數據，GRACE及GRACE-FO重力探測衛(wèi)星的陸地水儲量數據以及實測徑流數據。各數據集具體信息詳見表1。依據表1，選定2003—2017年作為研究時段，并采用雙線性插值法將各數據集的空間分辨率統(tǒng)一重采樣至0.1°×0.1°。

表1 蒸散發(fā)產品特征信息

2.3 SW雙源蒸散發(fā)模型的改進及應用

SW模型是由Penman Monteith模型發(fā)展而來的雙源蒸散發(fā)模型。SW模型的詳細信息可以參考Shuttleworth[11]和Hu[12]等人的研究，模型可用以下公式表示：

λET=CsPMs+CcPMc

(1)

PMs=ΔR+[ρcpVPD-ΔracRs]/

(2)

PMc=ΔR+[ρcpVPD-ΔracRs]/

(3)

式中，PMs、PMc—土壤蒸發(fā)、冠層蒸騰，mm；Cs、Cc—土壤表面阻力系數、冠層阻力系數；ρ—氣體密度，kg/m3；—飽和蒸氣壓-溫度曲線的斜率，kPa/k；γ—恒定常數，kPa/k；VPD—蒸氣壓差，kPa；R、Rs—冠層上方、土壤表面的凈輻射通量；rac—葉片至冠層高度的空氣動力阻力，s/m；ras—土壤表面至冠層高度的阻力，s/m；raa—冠層

高度至參考高度的阻力，s/m；rss—土壤表面阻力，s/m；rsc—冠層氣孔阻力，s/m。

本研究采用Hu等人[12]的研究在SW模型基礎上引入Ball-Berry氣孔導度模型，以解決在區(qū)域尺度上SW模型氣孔導度參數難以獲取的關鍵問題。其中，改進SW模型的rsc計算公式如下：

(4)

式中，g0、a1—經驗參數，可以參考Hu等人的研究；hs—葉片表面空氣相對濕度，%，采用空氣相對濕度代替；cs—CO2含量，ppm；Pn—光合速率，μmolm2s-1，是估算rsc中的關鍵因素，采用GLASS GPP產品替代。在本研究中，改進SW模型計算得到的黃河流域月尺度蒸散發(fā)簡稱為ETSW_GLASS。

2.4 基于水量平衡法的多源蒸散發(fā)產品精度驗證

在缺少實測蒸散發(fā)觀測數據的流域尺度下，基于降水、徑流和總蓄水量變化的水量平衡方法通常被廣泛用于蒸散發(fā)產品檢驗[13]。因此，本研究將采用流域水量平衡方法計算的黃河流域逐月蒸散發(fā)量作為“真實值”用于驗證ETGLASS、ETCR、ETSW_GLASS在黃河流域的數據準確性。水量平衡方法計算蒸散發(fā)公式如下：

ETWB(i)=P(i)-R(i)-TWSC(i)

(5)

TWSC(i)=TWSA(i+1)-TWSA(i-1)/2Δi

(6)

式中，i—時間尺度，月；ETWB—水量平衡方法計算的蒸散發(fā)量，mm/月；P、R—降水量、徑流量，mm/月；TWSA—總蓄水量異常；TWSC—總蓄水量變化。

2.5 基于TCH法的多源蒸散發(fā)產品不確定性評估

與傳統(tǒng)的誤差估計方法相比，TCH方法可以在未知區(qū)域蒸散真實值的情況下評估3種或更多不同數據集的不確定性[14]。本研究選取改進SW模型得到的蒸散發(fā)估算值(ETSW_GLASS)和兩個被廣泛使用的蒸散發(fā)產品(ETGLASS、ETCR)，采用TCH方法對其在流域尺度上的不確定性進行評價。評估過程的具體細節(jié)描述如下：

Xi=Xt+εi，?i=1，2，3，…，N

(7)

Yi，M=Xi-XR=εi-εR，?i=1，2，…，N-1

(8)

S=cov(Y)

(9)

S=J·R·JT

(10)

(11)

式中，Xt—蒸散發(fā)產品的“真實值”；Xi—蒸散發(fā)產品的模擬值；εi—相應的隨機誤差值；Y—M×(N-1)的方差矩陣；M—時間樣本數量；N—不同蒸散發(fā)產品數量；S—Y的協(xié)方差矩陣；R—N×N的噪聲協(xié)方差矩陣；J—單位矩陣。

然而，由于未知數(N×(N+1)/2)大于方程數(N×(N-1)/2)，式(10)無法求解。根據Galindo和Palacio等人研究，可以基于Kuhn-Tucker定理來解決這個問題。根據Kuhn-Tucker定理，目標函數可表示為：

(12)

其約束函數為：

H(r1N，…，rNN)=-|R|/|S|·K<0

(13)

(14)

結合式(11)計算得到蒸散發(fā)產品模擬結果的方差，最后取標準差作為不確定性的結果。

2.6 多源蒸散發(fā)產品精度評價指標

本研究采用統(tǒng)計參數包括決定系數(R2)、均方根誤差(RMSE)、平均絕對誤差(MAE)作為蒸散發(fā)產品精度檢驗的依據。各項指標的計算公式如下：

(15)

(16)

(17)

3 結果

3.1 蒸散發(fā)產品時空格局對比分析

圖2—3分別展示了2003—2017年的多年平均ETGLASS、ETCR和ETSW_GLASS在黃河流域的空間分布格局和年際變化趨勢。

圖2 2003—2017年多年平均蒸散發(fā)空間分布

由圖2可以看出，2003—2017年黃河流域基于像元尺度的多年平均ETGLASS、ETCR和ETSW_GLASS總體上都呈現出從東南向西北逐漸減小的趨勢，但不同蒸散發(fā)產品在一定程度上也表現出了較為不同的空間分布特征：ETGLASS在中部、南部和東部地區(qū)的多年平均蒸散發(fā)要明顯高于北部和西部地區(qū)；ETCR的多年平均蒸散發(fā)則表現為在西南地區(qū)較高，而在北部地區(qū)較低；ETSW_GLASS的多年平均蒸散發(fā)表現出在東南部較高，西北部較低。由圖3可見，2003—2017年ETGLASS在黃河流域的平均年蒸散發(fā)高于ETCR和ETSW_GLASS的平均年蒸散發(fā)值。ETGLASS和ETSW_GLASS平均年蒸散發(fā)都呈現出不顯著的增加趨勢，線性擬合斜率分別為2.79、2.08mm/年。ETCR平均年蒸散發(fā)呈下降趨勢，斜率為-1.62mm/年。ETGLASS、ETCR和ETSW_GLASS的多年平均蒸散發(fā)分別為502.28、408.46、442.82、414.4mm。

圖3 2003—2017年平均年蒸散發(fā)變化

綜合表明，ETGLASS、ETCR和ETSW_GLASS估算的蒸散發(fā)在黃河流域存在時空特征差異，且ETGLASS相較于ETCR和ETSW_GLASS的差異最大，差異主要表現在：ETGLASS的逐年蒸散發(fā)與多年平均蒸散發(fā)要顯著高于ETCR和ETSW_GLASS，且ETGLASS的多年平均蒸散發(fā)在北部和東部大片地區(qū)存在高估。

3.2 蒸散發(fā)產品精度評估

圖4展示了2003—2017年黃河流域ETGLASS、ETCR和ETSW_GLASS月平均蒸散發(fā)和基于水量平衡方法計算的流域月平均蒸散發(fā)“真實值”(ETWB)之間的比較結果。圖5展示了ETGLASS、ETCR和ETSW_GLASS月平均蒸散發(fā)與ETWB月平均蒸散發(fā)之間的決定系數、均方根誤差、平均絕對誤差。

圖4 水量平衡蒸散發(fā)與蒸散發(fā)產品的比較

圖5 水量平衡蒸散發(fā)與蒸散發(fā)產品的精度驗證

從圖5可以看出，盡管ETGLASS、ETCR和ETSW_GLASS的月平均蒸散發(fā)與基于水量平衡方法計算的流域月平均蒸散發(fā)“真實值”(ETWB)的對比存在高估和低估的情況，但總體上它們的月均蒸散發(fā)與月平均ETWB表現出相似的變化規(guī)律。圖6中R2、RMSE和MAE的值也驗證了這一點。T檢驗證明，在0.01的顯著水平上(p<0.01)，月平均蒸散發(fā)ETGLASS、ETCR和ETSW_GLASS與月平均ETWB之間存在顯著相關性，且R2、RMSE和MAE值較高，表明ETGLASS、ETCR和ETSW_GLASS在黃河流域的月平均蒸散發(fā)與月ETWB具有相似的精度。由圖6可知，ETSW_GLASS表現較好，具有相對較高的R2值以及較低的RMSE和MAE值。綜合來看，ETSW_GLASS表現最好，其次是ETCR和ETGLASS。以上結果表明通過引入GPP過程來改進的SW模型可以提高估算黃河流域蒸散發(fā)的精度。

圖6 黃河流域蒸散發(fā)產品不確定性空間分布

3.3 蒸散發(fā)產品不確定性評估

圖6顯示了2003—2017年黃河流域ETGLASS、ETCR和ETSW_GLASS的逐像元的月平均不確定性空間分布特性。

如圖6所示，基于TCH方法確定的ETGLASS、ETCR和ETSW_GLASS在黃河流域的不確定性范圍為0.04mm/月-32.53mm/月，其中不確定性最低的是ETSW_GLASS(4.66mm/月)，不確定性范圍為0.2～11.58mm/月，ETGLASS的不確定性最高為(11.61mm/月)，不確定性范圍為1.42～32.53mm/月，其次為ETCR(9.96mm/月)，不確定性范圍為4.42～27.55mm/月。圖6也表明，ETGLASS、ETCR和ETSW_GLASS的柵格尺度月平均不確定性具有相似的空間變化特性，總體上都呈現出黃河流域西部的不確定性普遍低于東部的特征。然而，ETGLASS、ETCR和ETSW_GLASS在月均不確定性的空間分布上也都表現出一定程度的不同：ETGLASS在中部地區(qū)具有較高的不確定性，而在東部和南部地區(qū)具有較低的不確定性；ETCR在中部和西部地區(qū)的不確定性較低，東南部地區(qū)的不確定性較高；ETSW_GLASS東南部地區(qū)的不確定性較高，北部和西部地區(qū)的不確定性較低?？偟膩碚f，2003—2017年黃河流域ETGLASS和ETCR的月均不確定性高于ETSW_GLASS。以上分析也證明，通過引入GPP過程改進SW模型可以降低蒸散發(fā)估算值的不確定性。

4 結論

本文綜合對比評價了ETGLASS、ETCR和ETSW_GLASS在不同時空尺度下蒸散發(fā)估算值在黃河流域的精度和不確定性，并對比分析了2003—2017年3種ET產品的時空分布格局。研究發(fā)現：

(1)ETGLASS、ETCR和ETSW_GLASS的蒸散發(fā)估算值存在差異，這些差異反映在柵格尺度年平均ETGLASS、ETCR和ETSW_GLASS的時空格局分布以及變化趨勢上。

(2)就R2、RMSE、MAE而言，月ETGLASS、ETCR和ETSW_GLASS的精度與基于水量平衡方法計算的ETWB相似，其中ETSW_GLASS的表現優(yōu)于ETGLASS和ETCR。

(3)ETGLASS、ETCR和ETSW_GLASS柵格尺度月均不確定性呈現出不同的空間變化特征，其中ETSW_GLASS的月均不確定性低于ETGLASS和ETCR。

后續(xù)可以在3種蒸散發(fā)產品對比評估的基礎上，對不同蒸散發(fā)產品的不確定做進一步的歸因分析以提高模擬精度。