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

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

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

0 引 言

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

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

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

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

Landsat8數據空間分辨率高、光譜信息豐富、定標精確，通過適當的方法就可以反演豐富的地面信息，是反演區域蒸散發的理想數據源，為區域蒸散量的監測提供支持。本文在總結文獻的基礎上，分析采用Landsat8數據反演區域蒸散發基本原理和過程并介紹SEBAL模型中所需參數的估算方法，耦合新的凈輻射能量估計方案改進SEBAL模型，并在半干旱地區進行了對比驗證。對利用SEBAL模型采用Landsat8遙感影像估算區域蒸散發的研究及應用有重要的指導意義。

1 模型基本原理

改進的SEBAL模型與SEBAL模型都是基于能量平衡方程，對于研究區域，網格內陸面能量平衡方程可以表達為

RN=G+H+λET

(1)

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

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

圖1 改進SEBAL模型流程

1.1 凈輻射計算

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

(2)

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

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

(3)

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

(4)

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

(5)

(6)

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

(7)

(8)

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

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

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)

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

G=0.3Rn(1-0.98NDVI4)

(10)

式中，NDVI為歸一化植被指數。

1.3 感熱通量(H)

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

(11)

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

空氣動力阻抗是摩擦風速、大氣穩定度和地表動量粗糙度的函數。運用風速廓線理論計算

(12)

摩擦風速

(13)

式中，k為卡爾曼常數；z1定義為地表面植被覆蓋處輻射能力轉變為感熱通量的高度，一般取值0.01 m；z2為不再受地表面粗糙影響的邊界層高度，一般為2.00 m；u200為200 m高度處的風速，m/s；ψm為動量傳輸修正度；ψh為熱量傳輸修正度；zzom為地表動量粗糙度，m。

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

2 研究應用

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

圖2 Walnut Gulch流域示意

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

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

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

表1 反演結果對比

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

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

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

3 結 語

SEBAL模型利用很少的輸入參數就能反演區域蒸散量，反演的過程中， 模型中所需參數的估算是影響反演精度的主要因素 ，本文嘗試利用Landsat8數據計算水汽含量及大氣透射率代替經驗公式來提高凈輻射能量的估算精度及空間分布的準確性、可靠性，從能量平衡的角度闡述了計算凈輻射耦合SEBAL模型的計算過程、步驟所需參數的估算方法，為模型在蒸散發反演中的應用提供幫助。