全天候區域地表蒸散發反演
——以黑河流域為例

廖前瑜 ，任 超 ，冷 佩 ，段四波，韓曉靜

（1. 桂林理工大學測繪地理信息學院，廣西桂林541004；2. 中國農業科學院農業資源與農業區劃研究所，北京100081）

0 引言

蒸散發是地表水分循環與能量平衡的重要組成分量，反映了大氣與地球表面的水分和能量交換［1］。由于全球氣候變化和人類對水資源的過度利用，導致了可利用水資源短缺和相關生態問題。因此，精確估算蒸散發在時空尺度上的變化，對農田灌溉管理、作物需水量估算、水資源優化配置等領域的研究具有重要意義。傳統的蒸散發估算大多基于氣象觀測的單點計算，盡管能得到相對準確的蒸散發量，但無法實現非均勻下墊面的蒸散發研究，而遙感對于區域蒸散發估算具有無可替代的優勢［2-3］。自20 世紀80 年代以來，眾多學者圍繞地表蒸散發遙感反演開展了大量工作，形成了許多切實可行的算法。這些算法大致可以分為基于能量平衡的方法［4-7］，基于水量平衡的方法［8］以及基于遙感地表溫度—植被指數特征空間的方法［9-13］。其中，基于遙感地表溫度—植被指數特征空間的方法近年來得到了最為廣泛的應用，主要得益于其簡單的原理和輸入以及它能夠避免用于計算水和熱傳遞的空氣動力學阻抗和表面阻抗的復雜參數化過程。應用地表溫度—植被指數特征空間方法反演地表蒸散發的一個重要前提是要確定特征空間的干濕邊。目前，確定地表溫度—植被指數特征空間干濕邊主要有基于傳統的空間統計的方法［14-15］和基于能量平衡理論干濕邊計算方法［13，16］。其中，基于傳統的空間統計的方法通常具有一定主觀性，且要求研究區氣象條件一致，導致其在區域尺度上難以應用。基于能量平衡的方法則需要輸入大量的參數，包括一些遙感難以獲取的參數如風速、空氣動力學阻抗及冠層阻抗等，中間參量過多會導致誤差的累積和傳遞。盡管基于能量平衡估算干濕邊的方法需要大量的輸入參數，但因其具有良好的普適性和物理基礎受到了研究學者們廣泛的關注。

最近Leng 等［16］提出了逐像元確定地表溫度—植被指數特征空間理論干濕邊的思想，避免了傳統特征空間方法在區域尺度應用的局限。由于在有云情況下無法獲取有效的植被指數和地表溫度參數，Leng 等［16］結合經典的Penman-Monteith 模型，發展了全天候地表蒸散發反演方法。然而，該研究僅在站點尺度進行了分析。在區域尺度上，個別站點觀測的氣象數據顯然無法很好地描述氣象因子的空間異質性。因此，文章在Leng等［16］研究的基礎上，基于目前常用的衛星陸表遙感產品和較高空間分辨率的格網氣象數據，深入探索區域尺度全天候蒸散發的反演。

1 數據

1.1 研究區與地面觀測數據

黑河流域（97°E～102°E，38°N～43°N）面積約為 14.29 萬 km2，位于我國西部甘肅省境內，流域氣候主要受到中高緯度西風帶環流和極地冷氣團的影響。該流域從上游到中、下游，以水為紐帶，形成了“冰雪／凍土—森林—草甸—人工/天然綠洲—荒漠—湖泊”的多元自然景觀。流域內寒區和干旱區并存，地表異質性顯著，尤其是山區冰凍圈和極端干旱的河流尾閭地區形成了鮮明對比。因此，準確獲取流域內復雜環境狀況下的地表蒸散發數據，對于深刻理解區域水循環和地氣之間能量交換，具有重要意義。此外，該文選擇黑河流域作為研究區的另一個重要原因是由中國科學院西北生態資源環境研究院主導開展一系列長期的生態水文實驗，積累了豐富的水文氣象觀測數據，為研究工作的開展奠定了良好的基礎。

該研究反演區域尺度全天候地表蒸散發，涉及的地面觀測數據包括自動氣象站氣象觀測數據和基于渦動的蒸散發觀測數據。由中國科學院西北生態資源環境研究院主導的“黑河流域生態—水文過程綜合遙感觀測聯合試驗”于2012 年5 月底至9 月中旬在黑河流域中游地區布置了多套渦動相關儀和自動氣象站，同步獲取蒸散發與氣象數據，為該文全天候蒸散的反演和驗證提供了優質數據源。這些觀測數據可從寒區旱區科學數據中心（http：//westdc.westgis.ac.cn/）免費下載。考慮到不同站點觀測時期的差異，該文選取研究區內2012 年6 月1 日至9 月15 日的數據開展研究。根據獲取的驗證數據選擇了4個具有代表性站點來探索區域地表全天候蒸散發反演，表1 給出了這4 個站點的編碼名稱、地表類型、高程、測量高度及經緯度。

圖1 研究區域的地表高程Fig.1 DEM in the study area

表1 黑河生態水文遙感試驗通量觀測站點信息Table 1 Information of Heihe ecological hydrological remote sensing test flux observation sites

1.2 遙感數據

與地面觀測數據對應，該文使用2012 年6 月1 日至9 月15 日的MODIS（Moderateresolution Imaging Spectroradiometer）陸表參數產品作為研究的遙感數據源。MODIS 是搭載在Terra 和Aqua 衛星上的一個重要傳感器，具有36 個波段的觀測通道，覆蓋中國全境，具有較高的空間分辨率和時間分辨率。目前，基于MODIS 觀測數據，已經發布了多種陸表和大氣產品數據，在全球范圍內得到了廣泛應用。該文使用的MODIS 產品及其參數如表2 所示。

表2 MODIS 數據詳情表Table 2 MODIS data detail

MODIS 產品數據來源于LAADS DAAC（Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center）（https：//ladsweb.modaps.eosdis.nasa.gov/search/）網站。該文使用MRT 和MRTSwath 工具對不同級別的MODIS 數據進行重投影等預處理。研究中MOD021KM、MOD03 和MOD05_L2 產品被用來反演地表短波凈輻射，其中MOD021KM 產品空間分辨率為250 m 和500 m 的大氣頂部輻射亮度和反射率數據被整合成空間分辨率為1 km 的數據；MOD35_L2 云掩膜產品被用來判斷MODIS 圖像的每個像元是否受到云的干擾；MOD09A1、MOD11A1、MOD13A2 和MOD15A2 作為全天候地表蒸散發反演的參數輸入，其中MOD09A1 和MOD15A2 在8 d 時間內被認為是不變的，MOD13A2 歸一化植被指數通過hants 濾波后插值得到日尺度的植被指數［17］。

1.3 格網氣象數據

研究所用的格網氣象數據為中國氣象局陸面數據同化系統CLDAS（CMA Land Data Assimilation System）的Verson 1.0 產品（簡稱CLDAS-v1.0），數據下載于中國氣象數據網（http：//data.cma.cn/）。CLDAS-v1.0 產品覆蓋范圍為東亞區域（0°～60°N，70°E～150°E），空間分辨率為0.0625°，時間分辨率為1 h。研究中用到的CLDAS-v1.0 氣象產品要素包括短波輻射、氣溫、風速和比濕。這些數據被重采樣到MODIS 數據1 km空間分辨率。以2012 年6 月1 日為例，圖2 顯示了北京時間上午11：00 黑河流域的短波輻射、氣溫、風速和比濕。

2 方法

該文在Leng 等［16］建立的全天候地表蒸散發反演方法基礎上，深入探索區域全天候地表蒸散發的反演。在該反演方法中，所有的像元根據是否存在有效的地表溫度數據被劃分為晴空和有云像元。其中，晴空像元的地表蒸散發利用逐像元地表溫度—植被指數特征空間方法進行反演，有云像元的地表蒸散發則基于Penman-Monteith 模型計算。基于MODIS 陸表產品和CLDAS 格網氣象數據，探索該全天候蒸散發反演方法在區域尺度上的可行性。圖3 為區域尺度全天候蒸散發反演的技術流程圖。

圖2 2012 年6 月1 日CLDAS 產品：（a）短波輻射；（b）氣溫；（c）風速；（d）比濕Fig.2 CLDAS products on June 1，2012

圖3 全天候蒸散發反演方法流程圖Fig.3 Flow chart of the all-sky evapotranspiration retrieval method

2.1 LST-NDVI 特征空間估算晴天像元蒸散發

區別于傳統基于空間信息的“地表溫度—植被指數”特征空間，Leng 等［16］發展了獨立于氣象和下墊面條件的逐像元兩階段“地表溫度—植被指數”特征空間新模式。圖4 為示意圖。該模式假設對于每一個晴空像元（如圖1 中的像元C），都對應存在一個虛擬的特征空間，該虛擬特征空間由 4 個極端條件（Ts，max，Ts，min，Tv，max，Tv，min）對應的地表溫度構成。其中，Ts，max和Tv，min兩點構成的斜線把梯形特征空間劃分為兩個三角形，分別為由 Ts，max，Ts，min，Tv，min3 點構成的下三角形，由 Ts，max，Tv，max，Tv，min3 點構成的上三角形。在實際應用中，只要已知這4 個極端地表情形（即干燥裸土、受水分脅迫的全植被覆蓋、供水良好的全植被覆蓋和濕潤裸土）對應的地表溫度，便可確定像元C 的虛擬特征空間，從而能夠根據像元C 在該虛擬特征空間中所處位置，估算蒸散發。該特征空間新模式的一個顯著優點是在實際應用中能夠避免傳統基于空間信息確定干邊和濕邊時對氣象和下墊面條件的苛刻要求，從而能夠直接反演大范圍蒸散發。

圖4 LST-NDVI 特征空間示意圖Fig.4 Schematic diagram of LST-NDVI feature space

LST-NDVI 特征空間假設NDVI 為某一值時，蒸散發隨LST 增大而線性遞增。當某像元的LST 和NDVI 落在下三角形中，土壤根層水分充足，植被蒸騰接近于潛在最大值，而表層土壤水分蒸散由像元處于三角形中的位置決定。在該情況下，土壤蒸發（ETs）和植被蒸騰（ETv）可以表示為：

式（1）～（2）中，Es，p和Ev，p分別為最大土壤蒸發量和最大植被蒸騰量。AB和AC為長度（圖4）。

類似地，當某像元的LST 和NDVI 落在上三角形中，則表層土壤已干，土壤蒸發為0，植被開始受到脅迫，植被蒸騰量取決于該像元處于三角形中的位置。在該情況下，土壤蒸發（ETs）和植被蒸騰（ETv）可以表示為：

其中 Ts，max，Ts，min，Tv，max和 Tv，min分別是干燥裸土，飽和裸土，水分脅迫全植被和充分灌溉的全植被的理論地表溫度（K），如圖4 所示。Moran 等［18］基于能量平衡原理，假設4 種極端條件給出了4 個理論地表溫度的計算公式：

式（7）中，Rn是地表凈輻射（W/m2），rs是表面阻抗。ra是空氣動力學阻抗，γ為干濕表常數（kPa/℃），Δ 為飽和水汽壓隨空氣溫度變化的斜率（kPa/℃），VPD為參考高度的空氣飽和水汽壓差（kPa），G是土壤熱通量（W/m2），Cp為空氣定壓比熱（J/（kg℃））。其中，rs和Rn參數的計算可參考 Leng 等［16］和 Todorovic 等［20］的文獻。

式（5）中，Rn，s和Rn，v分別為土壤表面凈輻射和植被表面凈輻射（W/m2），G是土壤熱通量（W/m2），Ta是空氣溫度（K）。Cv為空氣定容比熱（J/（℃m3）），γ為干濕表常數（kPa/℃），Δ 為飽和水汽壓隨空氣溫度變化的斜率（kPa/℃），VPD為參考高度的空氣飽和水汽壓差（kPa）。rcx是葉面氣孔幾乎完全關閉情況的冠層阻抗，設定為2 000 s/m；rcp是充足水分條件下的植被冠層阻抗，為100 與葉面積指數的比值；ra是空氣動力學阻抗。

地表蒸散發可以表示為植被蒸騰和土壤蒸發的和：

2.2 Penman-Monteith 公式估算有云條件像元蒸散發

在有云條件下，無法有效獲取LST 和NDVI 參數，因此使用經典的Penman-Monteith模型［19］估算地表蒸散發：

3 結果與分析

3.1 CLDAS 氣象數據分析

該文使用所選4 個站點自動氣象站實測的短波輻射、風速、氣溫、比濕對CLDAS 格網數據產品重采樣為1 km 空間分辨率后的格網數據進行分析，結果如圖5 所示。

圖5 CLDAS 數據產品與實測數據散點圖：（a）短波輻射；（b）風速；（c）氣溫；（d）比濕Fig.5 Comparison of CLDAS products and measured data

從圖5 可以看出CLDAS 短波輻射、風速、氣溫、比濕數據誤差較大，尤其是短波輻射和風速產品。這可能是由于CLDAS 氣象數據在空間上受到氣象站點疏密的影響，導致了在缺少氣象站點的西北地區數據精度下降。Vinukollu 等［21］使用相同的一套驅動數據對比驗證了P-M 模型［22］、P-T 模型［23］、表面能量平衡模型［5］，發現不同模型之間蒸散發反演的差異主要取決于凈輻射數據的選擇。考慮到短波輻射是地表蒸散發的能量來源，其精度高低不僅限制了地表蒸散發的反演結果，還是地表蒸散發的主要影響因子。基于此，該文采用Tang 等［24］的方法利用MODIS 數據估算短波凈輻射轉換為短波輻射來代替CLDAS 中的短波輻射，作為全天候蒸散發反演的輸入參數。

為了驗證MODIS 估算的短波凈輻射在黑河地區的可靠性，該文選取了黑河流域4個站點的2012 年6 月1 日至9 月15 日的地面有效測量數據進行對比。其中，地面觀測的短波凈輻射采樣頻率為10 min，MODIS 數據反演的短波凈輻射是分辨率為1 km 的瞬時值。另外，根據MOD35_L2 云產品判斷每個像元是否會受到云的影響。MOD35_L2 中認為是clear 和probably clear 的像元，判斷為晴空的像元；MOD35_L2 中認為是uncertain和cloudy 的像元，判斷為有云的像元。圖6 是該4 個站點的地表短波凈輻射觀測值和MODIS 估算的短波凈輻射散點圖，從圖6 可以看出，晴空條件下，基于MODIS 反演的短波凈輻射的均方根誤差小于30 W/m2；有云時，均方根誤差在70 W/m2左右。相比于利用CLDAS 數據獲取的短波輻射精度，基于MODIS 反演的短波凈輻射顯然具有更高的精度。

圖6 估算的地表短波凈輻射與實測的地表短波凈輻射散點圖（a）1 號站；（b）7 號站；（c）10 號站；（d）17 號站Fig.6 Comparison of estimated net surface short-wave radiation and measured data

3.2 全天候區域地表蒸散發反演與分析

為了驗證Leng 等［16］方法的可靠性，該文首先使用單點的實測氣象數據和MODIS數據反演地表蒸散發，并利用黑河流域4 個站點的2012 年6 月1 日至9 月15 日的地面有效測量蒸散發數據進行對比。驗證結果如圖7 所示。使用單點的氣象數據和MODIS數據反演的地表蒸散發RMSE 值在70 W/m2以內，具有較高的精度。這表明了利用Leng 等［16］的方法反演全天候地表蒸散發是可行的。在此基礎之上，利用CLDAS 短波輻射、風速、比濕、氣溫的產品與MODIS 估算了黑河流域的全天候地表蒸散發，并將反演結果與地面實測數據進行了比較，4 個站點的平均偏差值為30.8 W/m2，平均均方根誤差值為111.9 W/m2，誤差相對較大。采用Tang 等［24］的方法轉換得到短波輻射作為全天候蒸散發反演的輸入參數后，精度有了明顯的提高，其RMSE 在80 W/m2以內。

圖7 不同驅動數據估算的蒸散發和實測的蒸散發散點圖（a）1 號站、（b）7 號站、（c）10 號站、（d）17 號站Fig.7 Comparison of estimated evapotranspiration and measured data

基于能量平衡逐像元計算干濕邊的方法具有良好的普適性和物理基礎。然而大量的輸入參數造成了誤差的積累和傳遞。通過提高關鍵輸入參數的精度，能夠顯著地提高蒸散發的反演精度。

4 結論

該文以黑河流域為研究區，基于遙感數據和較高空間分辨率的格網氣象數據，實現了區域尺度地表蒸散發反演。在驗證分析了CLDAS 數據精度的基礎之上，利用MODIS數據估算的短波輻射代替CLDAS 短波輻射數據作為蒸散發反演模型的輸入參數，最后利用2012 年6 月1 日至9 月15 日的CLDAS 數據和MODIS 數據反演了黑河流域的全天候地表蒸散發，并對結果進行精度檢驗。研究表明，基于MODIS 遙感數據估算得到精度較高的短波輻射替代CLDAS 中的短波輻射產品，能顯著提高黑河地區蒸散發的反演精度。

雖然通過該文方法得到了精度較高的區域尺度全天候地表蒸散發，但研究仍然存在以下不足。首先是驗證站點較少，且相對空間距離較近，不利于對地表蒸散發的全面評價；其次是由于蒸散發反演和驗證過程中涉及不同尺度的數據源，該研究沒有充分考慮尺度效應對反演結果可能造成的影響。