基于Landsat8輻射模塊耦合SEBAL模型蒸散發(fā)估算

翟劭燚，王文種，劉九夫，2，王 歡，2，陸之昂，2

(1.南京水利科學(xué)研究院，江蘇南京210029；2.水利部應(yīng)對氣候變化研究中心，江蘇南京210029)

0 引 言

蒸散發(fā)是地表和大氣間水分、能量交換系統(tǒng)中的重要過程，是水文循環(huán)和能量循環(huán)的重要組成部分[1]；同時也是陸地表層水循環(huán)最大、最難估算的部分。傳統(tǒng)的估算蒸散發(fā)方法是基于點(diǎn)源的計(jì)算，研究方法包括：蒸滲儀法[2]、波文比能量平衡法[3]、渦度相關(guān)法[4]、蒸散速率與其他氣象要素的關(guān)系建立經(jīng)驗(yàn)公式推算法和水量平衡法[5]。傳統(tǒng)的氣候水文學(xué)方法對地表蒸散量的計(jì)算，僅限于離散點(diǎn)的觀測與小區(qū)域范圍內(nèi)估算。由于自然下墊面的異質(zhì)性，基于點(diǎn)尺度觀測數(shù)據(jù)在區(qū)域尺度上代表性較差，進(jìn)行空間尺度擴(kuò)展的準(zhǔn)確性未必能滿足需求。

隨著遙感技術(shù)的不斷發(fā)展，連續(xù)、動態(tài)地獲取大尺度區(qū)域地表信息的能力不斷提高，利用遙感手段對區(qū)域尺度與蒸散發(fā)相關(guān)的地表信息和遙感參數(shù)進(jìn)行監(jiān)測，利用遙感技術(shù)測量和估算區(qū)域尺度的蒸散發(fā)算法也成為新的技術(shù)手段，近年正蓬勃發(fā)展。遙感數(shù)據(jù)被引入各種估算模型，與氣候數(shù)據(jù)相結(jié)合，能夠?qū)崿F(xiàn)大空間尺度的蒸散量計(jì)算，對水文水資源的綜合管理以及氣候變化、農(nóng)業(yè)發(fā)展等都有重要意義。對于干旱、半干旱地區(qū)生態(tài)系統(tǒng)合理利用和分配水資源的研究，迫切需要深入了解不同植被覆蓋和土地利用條件下的蒸散發(fā)情況[6]。在水文、氣候、農(nóng)業(yè)和生態(tài)學(xué)中研究，對地表蒸散量的準(zhǔn)確估算將會產(chǎn)生巨大的經(jīng)濟(jì)、社會、生態(tài)和科學(xué)效益。

SEBAL模型是基于能量平衡原理通過遙感技術(shù)反演蒸散發(fā)的典型方法。模型假設(shè)研究區(qū)域存在冷熱極限像元，很好地考慮了干旱地區(qū)空間異質(zhì)；因此，SEBAL模型[7- 8](包括METRIC模型[9]和M-SEBAL模型[10])在干旱區(qū)遙感蒸散發(fā)估算中得到了廣泛的應(yīng)用[11]。然而，SEBAL模型計(jì)算中采用了許多經(jīng)驗(yàn)公式，特別是凈輻射和土壤熱通量等輸入?yún)?shù)的估算，具有區(qū)域適用性；這樣就會由于SEBAL模型實(shí)際應(yīng)用時的區(qū)域、時間、大氣條件、影像質(zhì)量等的不同使參數(shù)的普適性受到影響，最終會影響反演精度[12]。這些公式的地表參數(shù)在應(yīng)用到其他區(qū)域時，需要地面觀測數(shù)據(jù)核實(shí)驗(yàn)證。此外，為減少凈輻射估算的不準(zhǔn)確性，大多數(shù)以前的研究選用有限的遙感影像，這限制了時間尺度擴(kuò)展。因此，為使SEBAL模型更具有普遍適用性，需要研究機(jī)理性方法。

龍笛[13]研究表明，SEBAL模型中顯熱通量計(jì)算對于凈輻射(RN)不敏感，根據(jù)地表能量平衡方程的余項(xiàng)法計(jì)算得到的潛熱通量，應(yīng)該是敏感的。RN作為主計(jì)算潛熱通量的主要驅(qū)動力決定蒸散量(ET)的空間分布[17]。RN受許多因素影響，包括時間、地點(diǎn)、土地覆蓋的影響，大氣條件等，通常是額外地面太陽輻射、大氣輻射、大氣透射率函數(shù)、地表反照率、表面溫度和表面發(fā)射率等參數(shù)的函數(shù)[14-15]，SEBAL模型采用經(jīng)驗(yàn)公式簡化計(jì)算這些參數(shù)。例如，大氣透過率，這是由如氣溶膠因素約束的水汽、臭氧深度，潛熱通量對于大氣透過率敏感性較高，SEBAL模型中簡化為只有高度的函數(shù)[16]，這當(dāng)然會導(dǎo)致高估，如果不仔細(xì)考慮各種大氣因素操作。Bisht等人[14,17]提出的機(jī)械方案可能是一個合適的選擇，因?yàn)樗軌虿蹲降絉N空間分布且無需地面數(shù)據(jù)作為模型的輸入，也清楚地認(rèn)識到在空間上變化的輸入?yún)?shù)的要求。

Landsat8數(shù)據(jù)空間分辨率高、光譜信息豐富、定標(biāo)精確，通過適當(dāng)?shù)姆椒ň涂梢苑囱葚S富的地面信息，是反演區(qū)域蒸散發(fā)的理想數(shù)據(jù)源，為區(qū)域蒸散量的監(jiān)測提供支持。本文在總結(jié)文獻(xiàn)的基礎(chǔ)上，分析采用Landsat8數(shù)據(jù)反演區(qū)域蒸散發(fā)基本原理和過程并介紹SEBAL模型中所需參數(shù)的估算方法，耦合新的凈輻射能量估計(jì)方案改進(jìn)SEBAL模型，并在半干旱地區(qū)進(jìn)行了對比驗(yàn)證。對利用SEBAL模型采用Landsat8遙感影像估算區(qū)域蒸散發(fā)的研究及應(yīng)用有重要的指導(dǎo)意義。

1 模型基本原理

改進(jìn)的SEBAL模型與SEBAL模型都是基于能量平衡方程，對于研究區(qū)域，網(wǎng)格內(nèi)陸面能量平衡方程可以表達(dá)為

RN=G+H+λET

(1)

式中，RN是凈輻射量，W·m-2；G是土壤熱通量，W/m2，即地面與土壤間的熱量交換；H為感熱通量，W/m2；λ是水的汽化潛熱，W/(m2·mm2)。

利用晴朗天氣下的遙感數(shù)據(jù)反演相關(guān)地表參數(shù)(地表反照率、植被指數(shù)、地表發(fā)射率、地表溫度等)，結(jié)合相應(yīng)氣象及地表觀測資料(氣溫、風(fēng)速、平均高度等)，獲得區(qū)域范圍的凈輻射通量、土壤熱量，并根據(jù)遙感圖像中冷、熱像元點(diǎn)的選取確定地表溫度與溫度梯度差的線性關(guān)系，通過Monin-Obukhov相似理論迭代計(jì)算求得顯熱通量，從而根據(jù)余項(xiàng)法得到區(qū)域潛熱通量和瞬時蒸散發(fā)λET值，最后通過蒸發(fā)比不變原理[18]進(jìn)行時間尺度擴(kuò)展，求得區(qū)域日蒸散量或更長時間的蒸散總量。

圖1 改進(jìn)SEBAL模型流程

1.1 凈輻射計(jì)算

RN=(1-α)RS↓+RL↓-RL↑-(1-εs)RL↓

(2)

式中，α是地面反照率，W/m2；RS↓是下行的太陽短波輻射；RL↓是下行的長波輻射；RL↑是上行的長波輻射；是地面比輻射率。

太陽短波輻射能，它與太陽高度角、太陽輻射強(qiáng)度、大氣條件密切相關(guān)。Niemela[19]對比分析了6種不同晴空條件下的太陽短波輻射，認(rèn)為Zillman[20]提供的晴空條件下的短波輻射計(jì)算方法簡單，所需參數(shù)較少。公式為

(3)

式中，I0為太陽常數(shù)，取值1 367；θ為太陽天頂角，rad；1/R2為日地距離訂正因子，無量綱；τsw為短波大氣透射率，某特定地區(qū)上空大氣透過率受氣象因素影響很大，在較短時間內(nèi)，大氣的壓強(qiáng)、濕度、氣體密度可發(fā)生明顯的變化，透過率會因此而發(fā)生較大程度的改變。大氣單向透射率的值約為0.55～0.85，如果研究區(qū)的面積很小、海拔較低的時候可以用一個值代替，約取0.75，當(dāng)研究區(qū)面積較大時，必須將τsw轉(zhuǎn)化為研究區(qū)面上的計(jì)算，本文采用ASCE-EWRI推薦的公式[9]

(4)

式中，P為大氣壓強(qiáng)，kPa；K為無量綱校正系數(shù)，晴朗天氣下K=1.0；渾濁天氣下取K=0.5；w為大氣水汽含量，mm。王猛猛等[21]采用SWCVR算法應(yīng)用于Landsat8數(shù)據(jù)的水汽含量反演結(jié)果精度較高，得到Landsat8熱紅外大氣透過率比值和水汽含量的關(guān)系。即

(5)

(6)

式中，τ10、211為第10、11波段的大氣透過率；z為地面高程(m)。到達(dá)地面的長波輻射能量可由斯蒂芬-玻耳茲曼(Stefan-Boltzmann)方程計(jì)算。即

(7)

(8)

式中，σ為常數(shù)；εa為空氣比輻射率，Ta為空氣溫度；Ts為地表溫度。

在地表反演過程中，水汽是大氣透過率估算所必需考慮的主要因素。通常的做法是通過MODT-RAN、6S等大氣模型軟件模擬大氣透過率與水汽含量的關(guān)系。對于Landsat 8數(shù)據(jù)而言，大氣水汽含量很難從影像反演得到。宋挺等[22]使用與Landsat8影像同期過境時刻最為接近的采用多通道反演的MOD05-L2產(chǎn)品[23]，采用3次卷積內(nèi)插法對MOD05-L2水汽產(chǎn)品數(shù)據(jù)進(jìn)行重采樣處理成與Landsat8數(shù)據(jù)分辨率一致，并用MODIS MOD11-L2地表溫度產(chǎn)品與基于Landsat8反演的地表溫度結(jié)果進(jìn)行對比驗(yàn)證及精度分析，表明Juan C.Jiménez-Muoz劈窗算法(SW2)對輸入?yún)?shù)敏感性最低。

Ts=T10+1.378(T10-T11)+0.183(T10-T11)2-

0.268+(54.30-2.238ω)(1-ε)-

(129.20-16.40ω)ε

(9)

式中，T10、T11為第10、11波段的亮度溫度，K；ε為平均地表比輻射率。

1.2 土壤熱通量(G)

土壤熱通量是指進(jìn)入土壤和植被內(nèi)部的熱交換能量，表征表層土壤與下層土壤之間的熱交量換，是地表能量平衡中能量再分配的重要因素，是所有能量平衡方程分析中要涉及的關(guān)鍵因素。土壤熱通量與凈輻射熱通量、顯熱能量相比，土壤熱通量的值較小，其變化趨勢與凈輻射通量大致相同，采用參數(shù)化公式計(jì)算。即

G=0.3Rn(1-0.98NDVI4)

(10)

式中，NDVI為歸一化植被指數(shù)。

1.3 感熱通量(H)

感熱通量(Sensible Heat Flux)是指以對流和傳導(dǎo)形式進(jìn)入空氣中的那部分用于加熱空氣的熱能量，是由于溫度的梯度差異作用而造成的，其計(jì)算公式為

(11)

式中，ρ是空氣密度，1.21 kg/m3；Cp表示空氣定壓比熱，取值1004 J/(kg·K)；dT是零平面位移以上不同高度間的溫差，K；rah為空氣動力學(xué)的阻抗，s/m。rah是未知量，它們與不同梯度的溫度、地表粗糙度和風(fēng)速相關(guān)。

空氣動力阻抗是摩擦風(fēng)速、大氣穩(wěn)定度和地表動量粗糙度的函數(shù)。運(yùn)用風(fēng)速廓線理論計(jì)算

(12)

摩擦風(fēng)速

(13)

式中，k為卡爾曼常數(shù)；z1定義為地表面植被覆蓋處輻射能力轉(zhuǎn)變?yōu)楦袩嵬康母叨龋话闳≈?.01 m；z2為不再受地表面粗糙影響的邊界層高度，一般為2.00 m；u200為200 m高度處的風(fēng)速，m/s；ψm為動量傳輸修正度；ψh為熱量傳輸修正度；zzom為地表動量粗糙度，m。

由于表面加熱導(dǎo)致低層大氣內(nèi)部的不穩(wěn)定性，SEBAL根據(jù)Monin-Obukhov相似理論計(jì)算Monin-Obukhov相似長度及相應(yīng)大氣穩(wěn)定度修正因子ψm和ψh，對rah進(jìn)行校正，根據(jù)校正后的rah可重新求得dT，如此經(jīng)過若干次迭代，直至H達(dá)到穩(wěn)定，從而得到dT、Ts的經(jīng)驗(yàn)回歸系數(shù)a和b的最終值。根據(jù)區(qū)域內(nèi)TS的空間分布求得dT的空間分布，進(jìn)而求得每個像元的H，再根據(jù)余項(xiàng)法求解LE，通過時間尺度擴(kuò)展便得到一天的蒸散量。

2 研究應(yīng)用

Walnut Gulch流域(WG)位于美國東南部亞利桑那州的東南部((31°43′N， 110°41′W) ，是美國農(nóng)業(yè)研究中心(USDA-ARS)和陸地水文研究(NASA)的長期研究區(qū)域，面積為150 km2。研究區(qū)屬于圣佩德羅河上游流域(跨越索諾拉州、墨西哥州和亞利桑那州)，位于索諾蘭沙漠和齊瓦瓦沙漠之間的過渡地帶，屬于半干旱地區(qū)，流域的海拔高度為 1 250～1 885 m之間(見圖2)。

圖2 Walnut Gulch流域示意

WG被認(rèn)為是世界上觀測數(shù)據(jù)最為密集的半干旱區(qū)域，提供豐富的數(shù)據(jù)資料可供下載(詳見網(wǎng)址http:∥www.tucson.ars.ag.gov/dap/)。其包含豐富氣象水文數(shù)據(jù)[26]、地理信息[27-28]、通量[29]、遙感數(shù)據(jù)[28]等數(shù)據(jù)資料。這些資料為基于遙感衛(wèi)星的反演區(qū)域蒸發(fā)模型研究提供了先驗(yàn)知識和驗(yàn)證資料。該實(shí)驗(yàn)流域共設(shè)有三個自動觀測氣象站點(diǎn)，氣象數(shù)據(jù)主要有：氣壓、氣溫、風(fēng)速、相對濕度。區(qū)域內(nèi)的兩個地表通量觀測站點(diǎn)記錄的通量數(shù)據(jù)用于對比驗(yàn)證。

Landsat系列衛(wèi)星為世界上應(yīng)用范圍最廣應(yīng)用次數(shù)最多的衛(wèi)星之一，遙感衛(wèi)星影像的空間分辨率和光譜分辨率等也在不斷的提高。Landsat 8上攜帶有兩個主要載荷：OLI(Operational Land Imager，陸地成像儀)和TIRS(Thermal Infrared Sensor，熱紅外傳感器)。OLI包括了ETM+傳感器所有的波段，對波段進(jìn)行了重新調(diào)整，排除水汽吸收特征；OLI全色波段Band8波段范圍較窄，這種方式可以在全色圖像上更好區(qū)分植被和無植被特征；Landsat8上攜帶的TIRS載荷，將是有史以來最先進(jìn)，性能最好的TIRS。TIRS將收集地球兩個熱區(qū)地帶的熱量流失，目標(biāo)是了解所觀測地帶水分消耗，特別是美國西部干旱地區(qū)。數(shù)據(jù)可以通過http:∥glovis.usgs.gov/下載，本研究以2015337影像為例反演該區(qū)域地表蒸散發(fā)。

由計(jì)算結(jié)果提取草地和稀疏灌木站點(diǎn)的地表溫度和通量數(shù)據(jù)，稀疏灌木和草地兩個站點(diǎn)的landsat8估算值與地而觀測值的對比見表1。從表1中可以看出，草地和稀疏灌木兩個站點(diǎn)利用SEBAL模型估算的凈輻射都存在高估的問題，與觀測值相比偏差為175.7 W/m2和210.6 W/m2，而本文采用的計(jì)算結(jié)果偏差為7.4 W/m2和33.2 W/m2，計(jì)算的結(jié)果總體上要比SEBAL模型估算的結(jié)果更靠近1∶1線(見圖3)。圖3顯示，凈輻射計(jì)算值比SEBAL模型采用經(jīng)驗(yàn)公式估算的結(jié)果更接近真實(shí)值。

表1 反演結(jié)果對比

圖3 兩種模型估算的凈輻射能量反演結(jié)果對比

草地和稀疏灌木兩個站點(diǎn)的觀測潛熱通量分別為24.0 W/m2和30.4 W/m2，殘余法計(jì)算的潛熱通量分別為109.2 W/m2和0.35 W/m2，SEBAL模型相對應(yīng)于觀測的計(jì)算結(jié)果偏差為30.5 W/m2和18.6 W/m2，相對應(yīng)于殘余法的計(jì)算結(jié)果偏差為-54.7 W/m2和48.6 W/m2；而本文采用的耦合模型計(jì)算結(jié)果對應(yīng)于觀測的偏差為15.6 W/m2和-26.3 W/m2，相對應(yīng)于殘余法的計(jì)算結(jié)果偏差為-69.6 W/m2和3.7 W/m2，模型算法的可靠性還需要進(jìn)一步驗(yàn)證。2015337期反演區(qū)域地表蒸散發(fā)如圖4所示。

圖4 反演潛熱通量結(jié)果(單位：W/m2)

3 結(jié) 語

SEBAL模型利用很少的輸入?yún)?shù)就能反演區(qū)域蒸散量，反演的過程中， 模型中所需參數(shù)的估算是影響反演精度的主要因素 ，本文嘗試?yán)肔andsat8數(shù)據(jù)計(jì)算水汽含量及大氣透射率代替經(jīng)驗(yàn)公式來提高凈輻射能量的估算精度及空間分布的準(zhǔn)確性、可靠性，從能量平衡的角度闡述了計(jì)算凈輻射耦合SEBAL模型的計(jì)算過程、步驟所需參數(shù)的估算方法，為模型在蒸散發(fā)反演中的應(yīng)用提供幫助。