基于遙感的科爾沁沙地蒸散發(fā)時空動態(tài)

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

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

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

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

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

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

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

1 研究區(qū)概況及方法

1.1 研究區(qū)概況

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

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

表1 渦度相關(guān)塔測量項目及儀器

1.2 數(shù)據(jù)來源及處理

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

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

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

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

1.3 MEP模型

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

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

(1)

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

Rn=E+G+H

(2)

(3)

E=B(σ)H

(4)

(5)

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

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

(6)

(7)

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

1.4 MEP模型輸入?yún)?shù)計算方法

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

(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) 為劈窗算法的系數(shù);bk的取值隨著大氣水汽含量的變化而有所差異。使用劈窗協(xié)方差表示方差比法(Split Window Covariance-Variance Ratio,SWCVR)可以計算大氣水汽含量,計算公式如下:

W=a+bτj/τi

(9)

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

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

RH=e/es

(10)

(11)

(12)

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

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

(13)

1.5 SEBAL模型

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

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

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

(14)

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

土壤熱通量(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之間空氣動力學(xué)阻抗;Z1，Z2分別是指0.01 m和2 m高度處。

1.6 尺度擴展方法

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

(17)

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

2 結(jié)果與分析

2.1 MEP模型反演結(jié)果驗證

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

表2 地表參數(shù)及MEP模型結(jié)果估算誤差

表3 SEBAL模型估算誤差對比

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

圖1 MEP模型與SEBAL模型對比驗證結(jié)果

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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