基于遙感的科爾沁沙地蒸散發時空動態

張 鵬, 張圣微, 徐 冉, 高 露, 高文龍, 杜銀龍

(1.內蒙古農業大學 水利與土木建筑工程學院, 呼和浩特 010018;2.內蒙古自治區水資源保護與利用重點實驗室, 呼和浩特 010018; 3.內蒙古自治區農牧業大數據研究與應用重點實驗室, 呼和浩特 010018; 4.內蒙古自治區鄂爾多斯水文勘測局, 內蒙古 鄂爾多斯 017020)

蒸散發作為連接地表能量平衡與水量平衡的關鍵紐帶，又是聯系生態環境與水文循環的重要因素[1]。蒸散發會消耗全球地表約60%的降水量，并以水汽的形式返回到空氣中[2]，干旱、半干旱地區(科爾沁沙地)占據我國國土面積50%以上，這部分地區降水較少，年降水量均在500 mm以下，而降水量要遠低于潛在蒸散量[3]，因此，想要深入了解區域生態環境問題、水循環與水文過程，準確的蒸散發是必不可少的，同時為地區水資源開發利用，荒漠化防治提供科學依據。

傳統的ET測量技術(渦度相關、閃爍儀、波文比等)適用于均勻覆蓋下墊面的單點ET計算，對于區域來講，由于下墊面情況比較復雜，遙感方法被認為是計算區域ET唯一可行的方法[4]。經過數十年的發展，遙感估算蒸散發模型在理論與應用中取得了不少的成果，遙感模型在理論建設與遙感數據種類有較大的發展，其中主要的理論方法有：溫度—植被指數特征空間法[5]、將遙感數據與傳統方法相結合的遙感模型[6-7]、陸面過程與數據同化[8]、地表能量平衡模型等[9-10]，應用較多的遙感數據有：Landsat系列數據、MODIS數據、AVHRR數據以及高分數據等。

到目前為止，應用最廣泛的地表能量余項法是Bastiaanssen等人提出的地表能量平衡算法(SEBAL)[11-13]。它利用遙感的地表溫度、地表發射率、NDVI和較少的地面氣象觀測數據(即日照時數和風速)來估計區域蒸散發。與其他使用遙感數據估算地表蒸散量模型相比，SEBAL模型具有以下優點[14]：(1) 避免了收集大量地面輔助數據；(2) 可以在每個區域進行內部定標過程來估算感熱通量；(3) 基于Monin-Obukhov相似理論的內部迭代計算，避免了氣溫空間插值帶來的不確定性。因此，SEBAL模型已被許多國家和地區成功地應用于地表熱通量和ET的估算，并且經過多年的發展許多學者在不同方面對SEBAL模型進行了優化，ALLEN等人[15]在SEBAL模型的基礎上加入了坡度與坡向函數，更加全面的考慮到高程的變化，進而開發出METRIC (Mapping EvapoTranspiration at high Resolution with Internalized Calibration)模型，因此METRIC模型在估算山地蒸散發中有著較高的精度；TANG等人[16]通過建立廣義解析模型分析了感興趣區以及遙感數據像元大小對于SEBAI模型結果的影響；之后TANG等[17]人通過對比SEBAL模型與Ts-Ⅵ三角模型(surface temperature-vegetation index triangle models)發現SEBAL模型能更好地估計感熱通量和潛熱通量。

近年來，表面通量建模涌現出許多新的思路和方法，其中引入非平衡熱力學系統中最大熵增理論(maximum Entropy Production,MEP)就是其中之一。最大熵原理(maxEnt)和貝葉斯概率論是最大熵增理論的基礎，由Wang和Bras經過系統的數學推導、理論分析和多次測試而發展起來的一種新型蒸散模型[18-19]。MEP模型與傳統蒸散模型有較大的差異：(1) MEP模型是基于非平衡熱力學理論，結合能量守恒、邊界湍流等物理理論以及概率論而建立的蒸散模型[20]；(2) MEP模型在求解表面通量時，僅需凈輻射、地表溫度、表面比濕作為已知條件，不通過風速和地表粗糙度等物理量來求解；(3) MEP模型的前提條件為地表能量閉合，以能量平衡方程作為地表通量解的數學約束條件[21]。MEP模型的創建經過了嚴密的推理與驗證，起初Wang和Bras等[22]在估算陸面蒸散發過程中將最大熵增理論和地表能量平衡原理引入其中，并證明該方法的可行性；然后Wang和Bras等[23]認為植被在發生蒸騰的過程中內部的水分是飽和的，葉面氣孔是一律張開的，此時植被蒸騰作用最強，由此建立植被蒸騰模型；最后Wang和Bras[19]對MEP模型進行了完善與驗證，結果表明顯熱通量與感熱通量二者具有良好的線性關系。自從MEP的建立，MEP模型已經得到越來越多學者的關注，大部分學者使用渦度數據對MEP模型進行驗證，且驗證地點多數為濕潤地區，對于MEP模型在干旱半干旱地區的應用較少。Huang等[24]創新性的將MEP模型與地表輻射溫度遙感數據結合，估算了時空分辨率分別為3 h和1°的全球蒸散和地表熱通量，并首次給出了全球洋面熱通量，以及極地地區的地表熱通量；XU等[25]使用MEP模型結合站點渦度數據估算了亞馬遜地區的蒸散發量，并使用MEP模型評估了MODIS ET產品，發現兩者之間平均值差異較小，但是在空間上表現出較為明顯的差異。隨著MEP模型的不斷完善與發展，MEP模型在水文學方面還存在較大的應用空間。

本文首先利用Landsat8與氣象數據，估算MEP模型所需參數：凈輻射(Rn)、相對濕度(RH)、地表溫度(Ts)，其次使用MEP模型與SEBAL模型估算研究區2016年5—10月的ET，使用渦度數據以及SEBAL模型進行驗證并評估MEP模型在研究區的適用性，最后結合研究區不同土地覆被類型分析不同生態類型的ET變化特征，以期揭示科爾沁沙地不同生態類型區的水熱條件的分布特征，從而為區域水資源的合理利用及荒漠化治理提供理論支持。

1 研究區概況及方法

1.1 研究區概況

研究區位于內蒙古自治區通遼市科爾沁左翼后旗阿古拉鎮，地處科爾沁沙地南緣，地理坐標為122°33′00″—122°41′00″E，43°18′48″—43°42′24″N，面積約55 km2。總體上研究區地形自西向東、由南向北緩慢傾斜，西高東低。境內海拔最高為232 m，最低為186 m。研究區為固定沙丘、半固定沙丘、流動沙丘相結合，坨甸相間分布，中間地帶的草甸，農田和牧場鑲嵌交錯分布，沙丘、草甸、湖泊、農田和村莊分別占研究區總面積的54.5%，26.6%，5.2%，10.4%，3.3%。屬于典型的半干旱荒漠化地區，溫帶半干旱大陸性季風氣候，Ф20 cm蒸發皿多年平均水面蒸發量為1 400.3 mm，多年平均降水量為379.8 mm。研究區植被種類豐富，主要有差巴嘎蒿、沙蓬、少花蒺藜草、草麻黃、蟲實、冷蒿、蘆葦、羊草、蒲公英等。土壤類型有：砂土、壤砂土、壤土、黑鈣土、栗鈣土等。土壤中砂粒占80%以上，有的甚至達到100%。

沙丘、草甸試驗站點分別設有高10 m的微氣象觀測塔，在塔的不同高度分別布設了風、溫、濕梯度傳感器，并布設了輻射及渦度相關系統。兩個站點所安裝的均為開路式渦度通量觀測系統，附加了溫濕度傳感器用于修正渦動相關系統所測定的空氣溫濕度。站點的渦動相關系統架高分別為2.62 m，4.61 m，采集頻率均為10 Hz。具體觀測項目見表1。

表1 渦度相關塔測量項目及儀器

1.2 數據來源及處理

本文利用的遙感數據為Landsat-8/OLI影像(2016年5—10月)，獲取自USGS (https:∥earthexplorer.usgs.gov/)。取清晰度高、無云或云量較少的6期影像(5—9月)，每月一幅，之后對影像進行幾何校正、輻射定標、大氣校正等預處理工作。

氣象數據(氣溫、標準等壓面的氣壓、露點溫度)源自中國氣象數據網(http:∥data.cma.cn/)，使用通遼國家氣象站實測數據(2016年5—10月)。其他氣象數據(氣溫、風速、降水、土壤溫度等)分別源自研究區試驗站2016年(5—10月)沙丘與草甸兩個觀測站的實測數據。

渦度相關系統安裝的理想條件是下墊面水平均一，氣流穩定等，而在實際測量中，很難達到理想狀態，因此渦度數據在使用前需要預處理，處理步驟主要有：異常值及野點剔除、坐標旋轉修正、數據插補和求日平均等。

DEM數據為SRTM(ShuttleRadarTopographyMission)數據，經過簡單裁剪拼接等處理工作后即可使用，其空間分辨率為30 m，由于空氣溫度隨高程變化，因此需要DEM數據進行修正，數據來源于地理空間數據云(http:∥www.gscloud.cn/)。

1.3 MEP模型

Wang和Bras在估算蒸散發的過程中，引入最大熵增原理與地表能量余項法，這一方法使得估算蒸散發更加簡便[19-28]。本次研究使用他們的研究成果，并給出基于最大熵增原理的蒸散發模型的推理過程和公式。

對于無植被覆蓋的地表(裸地)，熵增函數D(H,E,G)的計算公式如下：

(1)

式中:H是顯熱通量(W/m2);Is為H的熱慣性參數(W·m2·K·s1/2);E是潛熱通量(W/m2);Ia為E的熱慣性參數(W·m2·K·s1/2);G為土壤熱通量(W/m2);Ie為G的熱慣性參數(W·m2·K·s1/2)。H,E和G是利用最大熵增原理,在Rn確定的前提下,運用拉格朗日乘數法極值化D(H,E，G),得到的計算公式:

Rn=E+G+H

(2)

(3)

E=B(σ)H

(4)

(5)

式中:Rn是凈輻射(W/m2),計算方法見公式(14);qs為表面比濕(g/kg);Ts為地表溫度(℃);σ為無量綱參數,可通過qs與Ts得出;I0為與顯熱通量獨立的簡潔表達式,在計算H的過程中需要利用二分法或牛頓迭代法,其計算精度需要提前確定并作為約束條件。

在MEP模型中認為，植被覆蓋的地表為單一封閉冠層，這時G相對于H和E可忽略不計，即Is=0，冠層上的能量平衡可以表示為：RnEv+H，Ev和H的計算公式如下:

(6)

(7)

式中:Ev為植被蒸騰作用的潛熱通量(W/m2),此時的Ts為葉面溫度;qs為葉面的表面比濕。

1.4 MEP模型輸入參數計算方法

MEP模型的優點在于輸入參數少，只需要地表溫度(Ts)、凈輻射(Rn)以及表面比濕(qs)，地表溫度的計算使用Ren等[29-30]應用通用劈窗算法針對Landsat8數據開發的模型，公式如下：

(8)

式中:Ti是第10波段(Band 10：10.6～11.19 μm)的亮度溫度;Tj是第11波段(Band 11：11.5～12.51 μm)的亮度溫度;ε為第10波段與第11波段的比輻射率平均值;Δε為第10波段與第11波段比輻射率的差值;bk(k=0,1,…，7) 為劈窗算法的系數;bk的取值隨著大氣水汽含量的變化而有所差異。使用劈窗協方差表示方差比法(Split Window Covariance-Variance Ratio,SWCVR)可以計算大氣水汽含量,計算公式如下:

W=a+bτj/τi

(9)

式中:W為大氣水汽含量(g/cm2);a和b為系數從模擬數據中獲得;τi和τj分別為i和j波段的大氣透過率。

為了方便后續模型的計算以及數據驗證，將表面比濕換算為相對濕度(RH)，本文將Peng等[31]使用MODIS數據估算大氣水汽含量方法移植到Landsat8數據，公式如下：

RH=e/es

(10)

(11)

(12)

式中:e為地面水汽壓(hpa);es為飽和水汽壓(hpa);Ta為大氣溫度(℃);Pa為大氣壓(hpa);W與qs之間有著良好的線性關系,使用通遼站的實測探空數據,得到W與qs的回歸方程,如式(13):

qs=0.001(-0.682W2+6.677W-1.1123)

(13)

1.5 SEBAL模型

SEBAL模型[11-13]在世界上得到了廣泛應用，在不同地區、氣候類型以及不同土地覆被類型下，有大量學者對SEBAL模型進行了驗證，相比其他遙感模型其精度相對較高，因此本文使用SEBAL模型進行對比驗證，計算過程如下：

Rn的計算公式如式(14)：

Rn=(1-α)RS↓+RL↓-RL↑-(1-ε0)RL↓

(14)

式中:Rs↓為入射短波輻射(W/m2);α為地表反照率;RL↓為入射長波輻射(W/m2);RL↑為出射長波輻射(W/m2);εo為寬波段表面熱發射。

土壤熱通量(G)計算公式如下：

(15)

顯熱通量(H)計算公式如下：

(16)

式中:ρair為空氣密度(kg/m3);Cp為空氣定壓比熱(取1 004.07 J/kg/K);Ts為地表溫度或冠層的表面溫度(K);Ta是2 m處的氣溫(K);dT為Z1與Z2兩個高度之間的溫差(K);rah是Z1與Z2之間空氣動力學阻抗;Z1，Z2分別是指0.01 m和2 m高度處。

1.6 尺度擴展方法

遙感模型中計算出的ET是衛星過境時刻的瞬時ET，在實際應用中，往往需要的是日ET，在本次研究中通過式(17)進行日尺度擴展[29]。

(17)

式中:ET24為日蒸散量(mm);λ為水分汽化潛熱指數;ETrF為衛星過境瞬時參考蒸散比(mm);Rn_24為衛星過境日的凈輻射通量(W/m2);G24為全天的土壤熱通量(W/m2)。

2 結果與分析

2.1 MEP模型反演結果驗證

利用2016年5—10月沙丘和草甸兩個觀測站的渦度相關數據以及氣象觀測數據與MEP模型估算結果、SEBAL模型估算結果、地表參數進行對比。為了更好地評價兩種模型，本文選用統計學中均方根誤差(RMSE)、相對誤差(MRE)與擬合系數(R2)作為評價指標，對比結果見表2與表3，圖1。

表2 地表參數及MEP模型結果估算誤差

表3 SEBAL模型估算誤差對比

由表2與表3得到，MEP模型反演G值的均方根誤差與相對誤差最大分別為61.50 W/m2，34%，G的擬合系數最小為0.57；SEBAL模型反演E值的均方根誤差最大為47.93 W/m2，G的相對誤差最大為15%，擬合系數最低為0.32；地表參數驗證結果相對較好，Rn的模擬結果誤差較小，其擬合系數最高為0.98，Ts的平均相對誤差最小為2%。由圖1看出：SEBAL模型與MEP模型之間反演結果整體誤差較小，其中MEP模型的H值偏高，G的誤差相對較小，個別日期的E值有著較大的差異。

圖1 MEP模型與SEBAL模型對比驗證結果

利用沙丘和草甸觀測數據，通過MEP模型分別反演了兩點的日ET，使用實測ET數據、遙感估算的ET數據以及Penmam-Monteith模擬ET之間進行對比，對比結果見圖2。

注：圖中ET模擬為Penmam-Monteith模擬的ET值。

從圖中看出MEP模型估算ET和實測ET變化趨勢基本一致，個別日期有所差異，遙感估算ET與實測ET以遙感反演ET值與觀測ET以及Penmam-Monteith模擬ET之間誤差較小，結果表明MEP模型能在不同的土地覆被類型模擬出精度較高的ET值。研究區ET總體波動較為平緩，表現為5—6月呈上升趨勢，7—8月相對較高，9—10月處于下降趨勢。沙丘的ET明顯低于草甸，隨著降水量的變化，沙丘的ET波動較為劇烈。降水對ET有著較大的影響，當有降雨事件發生的日期ET會隨之降低，出現一個波谷，降雨事件結束之后ET會有明顯升高，出現一個波峰。總體來說：在不同土地覆被類型下，MEP模型可以估算出相對精確的ET以及較為合理的變化趨勢。

2.2 MEP和SEBAL模擬ET的時空分布特征

本文通過MEP模型與SEBAL模型分別反演出研究區瞬時潛熱通量，通過日ET尺度擴展方法的到日ET值，日ET時空分布見附圖13，附圖14。

研究區日ET最大可以達到7 mm以上，主要集中在湖泊，湖泊相對較為穩定，日ET值在4～7 mm；農田與草甸由于植被與農作物的影響日ET變化較為明顯，尤其在6月、7月、8月，個別區域會有最大值出現，日ET保持在3～7 mm；沙丘日ET一直處于較低水平，保持在1～3 mm。

SEBAL模型與MEP模型得到的ET空間分布圖在5—7月差異較小，8—10月差異較大，二者差異較大的區域在半固定沙丘，誤差在0.5～1 mm；農田與草甸差異較小，誤差在0.2～0.6 mm。主要原因有：(1) MEP模型以能量閉合為邊界條件，而SEBAL模型會出現能量不閉合的情況，(2) 輸入參數不同，MEP模型所需參數少，只需要3個地表參數，而SEBAL模型輸入參數多，參數的不同引起結果差異，(3) 對于植被處理方法的不同：對于植被覆蓋的地表MEP模型認為土壤熱通量為零，直接計算顯熱通量，而SEBAL模型通過選取冷熱點的方式估算顯熱通量，5—7月研究區處于春季，植被較少，8—10月研究區處于夏季植被生長茂盛，因此導致8—10月差異較大。

為進一步分析MEP模型與SEBAL模型的差異以及ET時空變化，按照不同土地覆被類型對ET時空分布圖的ET平均值進行統計，得到各類土地覆被類型ET平均值變化曲線，見圖3。

圖3 2016年生長季不同土地覆被類型平均ET統計對比

隨著時間的變化不同土地覆被類型下的ET值波動情況有所差異，但是整體呈現先增長后降低的變化趨勢。5月所有土地覆被類型ET值均保持在較低水平，ET在2.5～4 mm；6月—8月大部分區域的ET值均有較為明顯的升高，而沙丘地由于水分補給不充足反而降低；在9—10月ET呈現明顯的下降趨勢，10月除湖泊仍保持在相對較高的水平，其他區域ET值較低，ET空間分布情況與5月相似。湖泊、草甸、農田整體呈較高水平，受月份影響較為明顯；沙丘波動較為明顯，受月份影響相對較小，沙丘類型不同，ET分布有較大差異，其中固定沙丘和半固定沙丘相似，而流動沙丘ET處于較低水平，尤其在6—8月差異較為明顯；村莊受人類活動影響較大，ET值分布相對較為集中并處于較低水平。

影響研究區ET時空分布的主要因素有：土壤類型不同，草甸地主要為黑鈣土、栗鈣土，沙丘地以砂土、壤砂土、壤土為主要的土壤類型，不同土壤類型下地表反照率有所差異，對太陽輻射的吸收也有差別[32]，導致ET的不同；植被類型不同，隨著生長季的到來植被的蒸騰作用所占蒸散發總量的比重越來越大，不同植被類型蒸騰作用有著較大差異，導致ET的不同[33]；土壤含水量以及供水條件不同，農田與草甸地的土壤含水量較高且水分供應充足，而沙丘地土壤含水量低，地下水埋深較深，除天然降水外，沒有其他水源供給，因此沙丘ET要低于農田與草甸[34]。

3 結 論

(1) 各地表參數以及ET估算值與地表觀測數據一致性較高，MEP模型與SEBAL模型對于日蒸散量的估算較為準確與誤差較小，表明兩個模型結合遙感數據均可為研究區提供合理的地表參數與ET模擬值。

(2) 通過對兩個模型的對比得出：兩個模型對于ET的估算較為一致，MEP模型估算的ET值要略低于SEBAL模型的估算值，個別日期兩個模型在沙丘地表現出較為明顯的差異，誤差在1.5 mm左右。導致差異的原因主要是由于邊界條件、輸入參數的不同以及對于植被處理方法的不同所導致。

(3) 研究區ET時空變化規律為：時間上ET呈現出5月較低，6—7月升高，8月達到最大，9月—10月呈下降趨勢，且降低幅度較大的變化規律；空間分布表現為：湖泊一直保持較高態勢，農田與草甸次之，沙丘最低的空間分布情況。影響ET主要因素有：土地類型以及土壤含水量不同、植被類型以及植被數量不同、人類活動。