基于改進的METRIC模型的農田潛熱通量估算

余 健， 姚云軍， 趙少華， 賈 坤， 張曉通， 趙 祥， 孫 亮

(1.北京師范大學遙感科學國家重點實驗室，遙感科學與工程研究院，北京 100875； 2.環境保護部衛星環境應用中心，北京 100094； 3.美國農業部水文與遙感實驗室，貝茨維爾 MD20705)

0 引言

潛熱通量是指地表發生土壤水分蒸發、水體或植被截留蒸發和植物體內水分的蒸騰過程中下墊面與大氣水分熱交換能量的總稱，在農業作物估產、大面積干旱監測和水資源管理中具有重要作用。我國是干旱發生頻率較多和水資源匱乏的國家，特別是華北平原地區，多年來農業糧食作物因缺水減產的情況較為嚴重。為合理地管理農業水資源和防旱抗旱，有必要開展農田潛熱通量的估算研究。

遙感高度融合了地表空間異質信息，能夠實現區域和田塊尺度上的潛熱通量估算和水分監測[1-2]。在眾多的遙感潛熱通量估算方法中，基于熱紅外遙感信息潛熱通量估算算法受到研究者的廣泛關注。1973年，Brown等[3]根據能量平衡原理和作物阻抗原理建立了作物阻抗—蒸散發模型，成為熱紅外遙感溫度應用于作物潛熱通量模型的理論基礎； 1983年，Seguin等[4]利用熱紅外遙感數據反演地表溫度估算日蒸發量，并作了進一步的闡釋，該方法在大尺度潛熱通量遙感估算中被廣泛應用； 隨后，Shuttleworth等[5]提出了經典的雙層模型，后于1991年進行了修正[6]，融合了植被和土壤的蒸散作用； Bastiaanssen等[7-8]開發出地表能量平衡算法(surface energy balance algorithm for land，SEBAL)模型，結合Landsat TM遙感影像反演了土耳其的蓋笛茲灌溉盆地的顯熱和潛熱通量，該模型在繼后的150多個研究項目中得到了應用。

最近，基于高空間分辨率的農田潛熱通量模型受到了許多學者的青睞。Allen等[9-10]于2007年提出了基于高分辨率和內在校準的蒸散估算法(mapping evapotranspiration at high resolution and with internalized calibration, METRIC)估算月度和季節性蒸散與潛熱通量，并將其應用于美國愛達荷州所需地下水開發量的計算和水權調度及計劃； 國內何磊等[11]和連晉姣等[12]利用METRIC和地表能量平衡系統(surface energy balance system, SEBS)模型估算了黑河流域的蒸散，驗證結果表明SEBS-METRIC方法可以用來估算黑河流域中游作物農田蒸散量，得出了在作物生長季內黑河流域中游地區蒸散量空間分布差異較大，大致在398～709 mm之間變化的結論。盡管METRIC模型物理意義明確，應用廣泛，但在具體區域潛熱通量反演實踐中受限制于地表觀測數據和氣象數據的缺失，往往導致反演精度不夠。

本文在分析以上方法與模型優缺點的基礎上，利用Landsat熱紅外數據獲取的地表溫度以及可見光和近紅外波段獲得的歸一化植被指數(normalized difference vegetation index，NDVI)改進地表粗糙度，提出基于地表粗糙度改進的METRIC模型來估算農田潛熱通量，并利用海河流域懷來和密云2個農田通量觀測站的通量觀測數據驗證精度，為開展大范圍的農田水分監測和農業灌溉提供一種簡易可行的方法。

1 材料與方法

1.1 研究區概況及數據源

本文的研究區海河流域位于E112°～120°，N35°～43°之間，流域面積約為32萬km2(圖1)。

圖1 研究區通量觀測站點空間分布圖Fig.1 Location of 2 flux tower sites throughout the study area

研究區地勢為西北高東南低，大致分高原、山地及平原3種地貌類型，西部為黃土高原和太行山區，北部為蒙古高原和燕山山區。研究區屬于半濕潤半干旱地區，年平均降水量為540 mm，年陸表蒸散量為470 mm[13-14]。土地覆蓋類型有旱地、草地、林地、灌木林、水田、濕地、水體、灘地和裸地等。本文先選取懷來和密云站做站點驗證分析，再以懷來站為中心，選取約24 km×24 km的樣帶區域進行空間分析，樣帶區域包含海河流域主要的土地覆蓋類型，能較好地代表海河流域的耗水特征。另外在樣帶尺度進行分析可以避免遇到云覆蓋影像的概率，從而提高影像的時間分辨率。

遙感數據采用Landsat5和Landsat8的原始影像和地表反射率數據，其中Landsat5熱紅外波段空間分辨率為120 m，Landsat8熱紅外波段空間分辨率為100 m，其余波段空間分辨率為30 m,數據來自美國地質調查局數據中心(http: //earthexplorer.usgs.gov/)。研究時間段共獲取影像58景，其中懷來站研究時間為2013—2014年，獲取28景影像； 密云站研究時間為2008—2010年，選取30景影像。因為云或云影對地表溫度的反演結果有很大影響，只選取了晴空條件下的影像進行遙感估算。

氣象和渦度相關觀測數據由寒區旱區科學數據中心提供(http: //westdc.westgis.ac.cn/)，其中氣象數據主要包括氣溫、相對濕度、風速、太陽輻射和地表溫度等。

本研究中在計算相關地表參數時需要對影像進行預處理，包括輻射定標、幾何糾正和圖像裁剪等。另外需注意的是要將影像的熱紅外波段重采樣成30 m×30 m，與可見光和近紅外波段保持一致。

1.2 地表參數計算

1.2.1 地表溫度反演

在經典METRIC模型中，對地表溫度的計算是采用熱紅外輻亮度通過對比輻射率的校正所得，而本文采用覃志豪的單窗算法[15-16]計算得到，即

Ts={a(1-C-D)+[(b-1)(1-C-D)+1]Tb-DTa}/C,

(1)

式中：Ts為地表溫度，K； a,b為常量；Tb為亮度溫度，K；Ta為大氣平均作用溫度，K；C和D為中間變量，由大氣透射率和地表輻射率計算可得。

1.2.2 地表凈輻射估算

瞬時地表凈輻射Rn是指地表吸收的太陽總輻射和地表有效輻射之差，其表達式為

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

(2)

式中：α為地表反照率；RS↓為下行短波輻射，W/m2；RL↓為下行長波輻射，W/m2；RL↑為上行長波輻射，W/m2；ε為地表比輻射率。各輻射通量的計算方法為

(3)

(4)

(5)

式中： Gsc為太陽常數(本文取1 367 W/m2)；θrel為太陽入射角，rad；τsw為大氣透射率；d為相對日地距離；εa為大氣發射率； σ為Stefan-Boltzmann常數，本文取5.67×10-8W/(m2K4)；Td為近地表空氣溫度，K。由于大氣長波輻射是來自整個大氣層，區域內差異較小，而近地表空氣溫度不好獲取，因此取水分含量較高的地表溫度代替近地表空氣溫度，或取冷點的地表溫度作替代。

1.3 METRIC模型

1.3.1 METRIC模型框架

傳統的METRIC模型[9]是基于能量余項法計算潛熱通量，即

LE=Rn-G-H,

(6)

式中：LE為瞬時潛熱通量，W/m2；G為瞬時土壤熱通量，W/m2；H為瞬時顯熱通量，W/m2。其中H和G分別為

(7)

G=(Ts-273.15)(0.003 8+0.007 4α)(1-0.98NDVI4)Rn,

(8)

式中：NDVI為歸一化植被指數；ρα為空氣密度，kg/m3；Cp為空氣定壓比熱，J/(kgK)； dT為地氣溫差，通常指z1和z2高度之間的溫差(一般取z1=0.01 m，為地表高度，z2=2 m,為常規氣象站高度)；rah為熱量傳輸空氣動力學阻抗，s/m。

METRIC模型假設dT與Ts成正比，即

dT=cTs+d,

(9)

式中c和d為經驗系數，由遙感圖像上的極端點(冷熱點)的相關參數計算得到。

在METRIC模型中，假設研究區土壤水分存在2種情況： ①極端干燥環境下，沒有水分的蒸發，潛熱通量基本為0，而顯熱通量達到最大，即熱點，本文取溫度較高(圖像最高溫度的0.95倍)且沒有植被覆蓋(NDVI<0.3)的像元； ②濕潤環境下，此時水分蒸發達到最大，而顯熱通量達到最小值，即冷點，本文取溫度較低(一般約取1.05倍的圖像中的最低溫度)且植被茂密(NDVI>0.7)的像元。

故熱點的顯熱通量和冷點的潛熱通量分別為

Hhot=Rn,hot-Ghot,

(10)

LEcold=Rn,cold-Gcold-Hcold,

(11)

式中：Hhot和Hcold分別為熱點和冷點的顯熱通量，W/m2；Rn,hot和Rn,cold分別為熱點和冷點的瞬時凈輻射通量，W/m2；Ghot和Gcold分別為熱點和冷點的土壤熱通量，W/m2；LEcold為冷點的潛熱通量，W/m2。

METRIC算法通過迭代計算熱像元rah，結合冷像元信息求解影像內系數c和d。初始假設大氣中性穩定，rah計算公式為

(12)

式中： k為von Karman常數，本文取0.41；u*為摩擦風速，m/s，在中性大氣條件下的計算公式為

(13)

式中：u200為200 m高空處風速，m/s；zom為地表粗糙度。

1.3.2 地表粗糙度參數改進

在對地表粗糙度的估算上，研究學者們發展了許多方案，如利用粗糙度與葉面積指數(leaf area index，LAI)之間的關系來計算地表粗糙度，即

zom=0.018LAI。

(14)

Moran等[17]研究發現，地表粗糙度和光譜反射有很大的相關性。因此，本研究通過建立zom與NDVI的經驗關系估算zom。本文采用粗糙度zom與NDVI之間的關系進行粗糙度估算，關系表達式為

zom=exp(a1+b1NDVI)。

(15)

確定回歸參數a1和b1包括以下3個步驟： ①在研究區域選取一定數量的樣本區域，測量出樣本區域的平均株高h，根據公式zom=0.123h算出樣本區域的粗糙度[17]； ②在計算得到的NDVI影像上找出對應樣本區域的像元，找出像元的NDVI值； ③利用非線性擬合的方法算出回歸參數。本文擬合得到的粗糙度和NDVI關系為

zom=exp(-6.57+7.33NDVI)。

(16)

研究區域的zom取值主要集中在0.01～0.03 m之間，均值在0.018 m左右相當于植被株高在15 cm左右，這對于對應測量的時間段樣本區下墊面植被類型是玉米來說是合理的。

2 結果與討論

2.1 地表溫度和地表凈輻射驗證

為了驗證Ts反演的精度，將其與站點觀測值進行相關分析。從圖2(a)可以看出，反演得到的Ts值的平均誤差是-1.65 K，均方根誤差為3.4 K，相關系數平方為0.93，經統計檢驗，達到顯著水平(p<0.05)。理想狀況下遙感地表溫度反演誤差為1 K之內，然而由于受到觀測數據誤差、復雜地表、尺度差異等多種因素導致誤差在3 K左右，可見單窗法反演地表溫度的精度能滿足METRIC模型的要求。

基于地表輻射平衡公式計算的凈輻射與通量觀測數據進行對比，以此檢驗算法估算的有效性。由圖2(b)，估算的凈輻射值的偏差是8.28 W/m2，均方根誤差為32.94 W/m2，相關系數平方為0.95(p<0.05)，由此可知遙感估算凈輻射通量結果是可以接受的。

(a) 地表溫度(b) 地表凈輻射通量

圖2地表溫度及地表凈輻射通量與站點觀測值散點圖

Fig.2Scatter-plotofretrievedversusobservedTsandestimatedversusobservedRn

2.2 潛熱通量驗證與算法比較

根據地表粗糙度改進的METRIC模型先估算顯熱通量H，進而推算潛熱通量LE,模擬值與觀測值的關系如圖3所示。

(a) 改進METRIC模型模擬顯熱通量 (b) 改進METRIC模型模擬潛熱通量

(c) 改進METRIC模型模擬地表通量 (d) 傳統METRIC模型模擬地表通量

圖3模擬值與站點觀測值散點圖

Fig.3Scatter-plotofestimatedandground-measuredvalues

圖3(a)所示，顯熱通量的相關系數平方為0.89，偏差為1.48 W/m2，模型模擬的結果顯著性明顯。最后利用能量余項法作差得到潛熱通量，如圖3(b)所示，模型計算的潛熱值相關系數平方為0.97(p<0.05)，偏差為2.31 W/m2，均方根誤差為31.06 W/m2，進一步表明METRIC模型估算該區域的農田潛熱通量具有較高的模擬精度。本文同時也利用懷來站的部分影像進行傳統的METRIC模型估算，圖3(c)是地表粗糙度改進的METRIC模型模擬的地表通量估算效果，圖3(d)是傳統METRIC模型模擬的地表通量估算效果。可以看出，傳統的METRIC模型模擬的顯熱通量和潛熱通量結果驗證散點圖中值分散，估算結果與站點觀測值的誤差在50～100 W/m2左右，而地表粗糙度改進的METRIC模型模擬的結果與觀測值的誤差基本在50 W/m2以下，并且經過地表粗糙度的改進后，相關系數平方從0.9提高到0.96，偏差和均方根誤差都顯著減小，這說明地表粗糙度改進的METRIC模型估算結果更為合理準確，同時也更適用于海河流域的潛熱通量估算。

2.3 潛熱通量制圖與分析

選取懷來觀測區域2014年第238日遙感影像對瞬時潛熱通量進行制圖與分析，由于水體潛熱通量估算影響因素較多，本研究采取Priestley-Taylor模式對水體潛熱通量進行了單獨計算。從圖4(a)可以看出，潛熱通量的差異明顯，變化范圍大致在30～830 W/m2，均值約為470 W/m2。水體的潛熱通量較大約在700～750 W/m2之間，而處于河邊的灘地由于土壤含水量較高，潛熱通量也較高，主要集中在550 W/m2左右； 農田的潛熱通量呈偏峰型，眾數在450 W/m2左右，此時處于夏季，玉米地處于生長季； 林地日潛熱通量眾數在550 W/m2左右，這主要是由于夏季植被覆蓋率較高； 而處于農田旁邊的居民建設用地潛熱通量也較高(約為480 W/m2)，是由于居民往往在居住地附近種植瓜果蔬菜； 另外官廳水庫邊草地的潛熱通量主要集中在350 W/m2。圖4(b)與圖4(a)空間格局相同，但整體模擬值偏低。

由圖4(c)對比后發現，地表粗糙度改進的METRIC模型模擬研究區的潛熱通量增大，這是由于傳統的METRIC模型地表粗糙度存在較大偏差導致潛熱通量偏小。同時發現，模型的改進對農田和的影響比較大而對水體的影響比較小，這是因為本文主要是對模型中地表粗糙度和地表溫度進行改進以及對坡度坡向進行糾正。此外在第238日這天，站點的觀測值為487.7 W/m2，傳統METRIC估算值為403.3 W/m2，地表粗糙度改進后的估算值為501.3 W/m2，粗糙度改進后的模型估算值和觀測站觀測值之間誤差更小，說明地表粗糙度改進的METRIC模型模擬更加理想。

(a) 改進METRIC模型模擬潛熱通量空間分布 (b) 傳統METRIC模型模擬潛熱通量空間分布(c) 本文模型與傳統模型模擬潛熱通量差值

圖4改進METRIC模型與傳統模型模擬潛熱通量分布及其差值分布

Fig.4SpatialdistributionofLEusingimprovedMETRIC,traditionalMETRICandtheirdifferences

3 結論

利用Landsat熱紅外數據獲取的地表溫度以及可見光-近紅外數據獲得的歸一化植被指數(NDVI)，改進地表粗糙度，提出基于地表粗糙度改進的METRIC模型來估算農田潛熱通量，并利用海河流域懷來和密云2個農田通量觀測站的通量觀測數據驗證該算法的精度，得出的結論如下：

1)基于地表粗糙度改進的METRIC模型與實測潛熱通量之間具有較好的相關性，相關系數平方可達0.97，經統計檢驗，達到顯著水平(p<0.05)。

2)盡管基于地表粗糙度改進的METRIC模型和傳統的METRIC模型都能模擬農田潛熱通量，但是基于地表粗糙度改進的METRIC模型的模擬值與實測潛熱通量的相關性更好，表明了改進模型的優越性。

3)從潛熱通量空間分布來看，地表粗糙度改進的METRIC模型模擬研究區的潛熱通量值比傳統METRIC模型模擬精度要高，模擬效果更加可靠。然而，由于數據獲取的局限性，本文只采用了海河流域2個站點通量數據進行這些模型的驗證與比較，在其他區域的驗證與比較仍需要進一步研究。