基于改進(jìn)的METRIC模型的農(nóng)田潛熱通量估算

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

(1.北京師范大學(xué)遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遙感科學(xué)與工程研究院，北京 100875； 2.環(huán)境保護(hù)部衛(wèi)星環(huán)境應(yīng)用中心，北京 100094； 3.美國(guó)農(nóng)業(yè)部水文與遙感實(shí)驗(yàn)室，貝茨維爾 MD20705)

0 引言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況及數(shù)據(jù)源

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

圖1 研究區(qū)通量觀測(cè)站點(diǎn)空間分布圖Fig.1 Location of 2 flux tower sites throughout the study area

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

遙感數(shù)據(jù)采用Landsat5和Landsat8的原始影像和地表反射率數(shù)據(jù)，其中Landsat5熱紅外波段空間分辨率為120 m，Landsat8熱紅外波段空間分辨率為100 m，其余波段空間分辨率為30 m,數(shù)據(jù)來(lái)自美國(guó)地質(zhì)調(diào)查局?jǐn)?shù)據(jù)中心(http: //earthexplorer.usgs.gov/)。研究時(shí)間段共獲取影像58景，其中懷來(lái)站研究時(shí)間為2013—2014年，獲取28景影像； 密云站研究時(shí)間為2008—2010年，選取30景影像。因?yàn)樵苹蛟朴皩?duì)地表溫度的反演結(jié)果有很大影響，只選取了晴空條件下的影像進(jìn)行遙感估算。

氣象和渦度相關(guān)觀測(cè)數(shù)據(jù)由寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心提供(http: //westdc.westgis.ac.cn/)，其中氣象數(shù)據(jù)主要包括氣溫、相對(duì)濕度、風(fēng)速、太陽(yáng)輻射和地表溫度等。

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

1.2 地表參數(shù)計(jì)算

1.2.1 地表溫度反演

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

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為中間變量，由大氣透射率和地表輻射率計(jì)算可得。

1.2.2 地表凈輻射估算

瞬時(shí)地表凈輻射Rn是指地表吸收的太陽(yáng)總輻射和地表有效輻射之差，其表達(dá)式為

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

(2)

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

(3)

(4)

(5)

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

1.3 METRIC模型

1.3.1 METRIC模型框架

傳統(tǒng)的METRIC模型[9]是基于能量余項(xiàng)法計(jì)算潛熱通量，即

LE=Rn-G-H,

(6)

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

(7)

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

(8)

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

METRIC模型假設(shè)dT與Ts成正比，即

dT=cTs+d,

(9)

式中c和d為經(jīng)驗(yàn)系數(shù)，由遙感圖像上的極端點(diǎn)(冷熱點(diǎn))的相關(guān)參數(shù)計(jì)算得到。

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

故熱點(diǎn)的顯熱通量和冷點(diǎn)的潛熱通量分別為

Hhot=Rn,hot-Ghot,

(10)

LEcold=Rn,cold-Gcold-Hcold,

(11)

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

METRIC算法通過(guò)迭代計(jì)算熱像元rah，結(jié)合冷像元信息求解影像內(nèi)系數(shù)c和d。初始假設(shè)大氣中性穩(wěn)定，rah計(jì)算公式為

(12)

式中： k為von Karman常數(shù)，本文取0.41；u*為摩擦風(fēng)速，m/s，在中性大氣條件下的計(jì)算公式為

(13)

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

1.3.2 地表粗糙度參數(shù)改進(jìn)

在對(duì)地表粗糙度的估算上，研究學(xué)者們發(fā)展了許多方案，如利用粗糙度與葉面積指數(shù)(leaf area index，LAI)之間的關(guān)系來(lái)計(jì)算地表粗糙度，即

zom=0.018LAI。

(14)

Moran等[17]研究發(fā)現(xiàn)，地表粗糙度和光譜反射有很大的相關(guān)性。因此，本研究通過(guò)建立zom與NDVI的經(jīng)驗(yàn)關(guān)系估算zom。本文采用粗糙度zom與NDVI之間的關(guān)系進(jìn)行粗糙度估算，關(guān)系表達(dá)式為

zom=exp(a1+b1NDVI)。

(15)

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

zom=exp(-6.57+7.33NDVI)。

(16)

研究區(qū)域的zom取值主要集中在0.01～0.03 m之間，均值在0.018 m左右相當(dāng)于植被株高在15 cm左右，這對(duì)于對(duì)應(yīng)測(cè)量的時(shí)間段樣本區(qū)下墊面植被類型是玉米來(lái)說(shuō)是合理的。

2 結(jié)果與討論

2.1 地表溫度和地表凈輻射驗(yàn)證

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

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

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

圖2地表溫度及地表凈輻射通量與站點(diǎn)觀測(cè)值散點(diǎn)圖

Fig.2Scatter-plotofretrievedversusobservedTsandestimatedversusobservedRn

2.2 潛熱通量驗(yàn)證與算法比較

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

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

(c) 改進(jìn)METRIC模型模擬地表通量 (d) 傳統(tǒng)METRIC模型模擬地表通量

圖3模擬值與站點(diǎn)觀測(cè)值散點(diǎn)圖

Fig.3Scatter-plotofestimatedandground-measuredvalues

圖3(a)所示，顯熱通量的相關(guān)系數(shù)平方為0.89，偏差為1.48 W/m2，模型模擬的結(jié)果顯著性明顯。最后利用能量余項(xiàng)法作差得到潛熱通量，如圖3(b)所示，模型計(jì)算的潛熱值相關(guān)系數(shù)平方為0.97(p<0.05)，偏差為2.31 W/m2，均方根誤差為31.06 W/m2，進(jìn)一步表明METRIC模型估算該區(qū)域的農(nóng)田潛熱通量具有較高的模擬精度。本文同時(shí)也利用懷來(lái)站的部分影像進(jìn)行傳統(tǒng)的METRIC模型估算，圖3(c)是地表粗糙度改進(jìn)的METRIC模型模擬的地表通量估算效果，圖3(d)是傳統(tǒng)METRIC模型模擬的地表通量估算效果?？梢钥闯?，傳統(tǒng)的METRIC模型模擬的顯熱通量和潛熱通量結(jié)果驗(yàn)證散點(diǎn)圖中值分散，估算結(jié)果與站點(diǎn)觀測(cè)值的誤差在50～100 W/m2左右，而地表粗糙度改進(jìn)的METRIC模型模擬的結(jié)果與觀測(cè)值的誤差基本在50 W/m2以下，并且經(jīng)過(guò)地表粗糙度的改進(jìn)后，相關(guān)系數(shù)平方從0.9提高到0.96，偏差和均方根誤差都顯著減小，這說(shuō)明地表粗糙度改進(jìn)的METRIC模型估算結(jié)果更為合理準(zhǔn)確，同時(shí)也更適用于海河流域的潛熱通量估算。

2.3 潛熱通量制圖與分析

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

由圖4(c)對(duì)比后發(fā)現(xiàn)，地表粗糙度改進(jìn)的METRIC模型模擬研究區(qū)的潛熱通量增大，這是由于傳統(tǒng)的METRIC模型地表粗糙度存在較大偏差導(dǎo)致潛熱通量偏小。同時(shí)發(fā)現(xiàn)，模型的改進(jìn)對(duì)農(nóng)田和的影響比較大而對(duì)水體的影響比較小，這是因?yàn)楸疚闹饕菍?duì)模型中地表粗糙度和地表溫度進(jìn)行改進(jìn)以及對(duì)坡度坡向進(jìn)行糾正。此外在第238日這天，站點(diǎn)的觀測(cè)值為487.7 W/m2，傳統(tǒng)METRIC估算值為403.3 W/m2，地表粗糙度改進(jìn)后的估算值為501.3 W/m2，粗糙度改進(jìn)后的模型估算值和觀測(cè)站觀測(cè)值之間誤差更小，說(shuō)明地表粗糙度改進(jìn)的METRIC模型模擬更加理想。

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

圖4改進(jìn)METRIC模型與傳統(tǒng)模型模擬潛熱通量分布及其差值分布

Fig.4SpatialdistributionofLEusingimprovedMETRIC,traditionalMETRICandtheirdifferences

3 結(jié)論

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

1)基于地表粗糙度改進(jìn)的METRIC模型與實(shí)測(cè)潛熱通量之間具有較好的相關(guān)性，相關(guān)系數(shù)平方可達(dá)0.97，經(jīng)統(tǒng)計(jì)檢驗(yàn)，達(dá)到顯著水平(p<0.05)。

2)盡管基于地表粗糙度改進(jìn)的METRIC模型和傳統(tǒng)的METRIC模型都能模擬農(nóng)田潛熱通量，但是基于地表粗糙度改進(jìn)的METRIC模型的模擬值與實(shí)測(cè)潛熱通量的相關(guān)性更好，表明了改進(jìn)模型的優(yōu)越性。

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