蔣 磊 蔡甲冰 張寶忠 許 迪 魏 征
(1.天津農(nóng)學(xué)院水利工程學(xué)院,天津 300392; 2.中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室,北京 100038;3.國家節(jié)水灌溉北京工程技術(shù)研究中心,北京 100048)
作物蒸散發(fā)(Evapotranspiration,ET)是農(nóng)業(yè)水管理的重要組成部分,合理估算區(qū)域作物蒸散發(fā)對作物用水效率評價和精量灌溉決策具有重要的意義[1-2]。隨著遙感技術(shù)的發(fā)展,基于遙感數(shù)據(jù)和地面觀測數(shù)據(jù)建立的蒸散發(fā)模型為科學(xué)、準(zhǔn)確獲取區(qū)域作物ET信息提供了一種有效的途徑和方法[3-4]。
在區(qū)域作物遙感蒸散發(fā)估算方法中,發(fā)展較為成熟的是基于地表垂直方向能量平衡的思路,即將蒸散發(fā)(潛熱)作為能量平衡中的余項求出[5-6]。根據(jù)不同的輻射通量和湍流交換方式,能量平衡余項法可分為單源模型和雙源模型。單源模型又稱“大葉模型”,即不區(qū)分土壤蒸發(fā)和植被蒸騰,將土壤和植被視為一個均一的整體來考慮潛熱通量,常見的單源模型包括SEBAL(Surface energy balance algorithm for land)模型[7]和SEBS(Surface energy balance system)模型[8]等。與單源模型不同,雙源模型分別計算土壤和冠層兩個組分的凈輻射、顯熱和潛熱,進而能夠區(qū)分土壤蒸發(fā)和植被蒸騰。NORMAN等[9]開發(fā)的TSEB(Two-source energy balance)模型是雙源模型的代表。YANG等[10]在前人研究的基礎(chǔ)上提出了混合雙源蒸散發(fā)模型(Hybrid dual-source scheme and trapezoid framework-based evapotranspiration model,HTEM),該模型已在國內(nèi)外多個地區(qū)進行驗證,并取得較高的精度[3,11],與TSEB、MOD16-ET產(chǎn)品的對比表明,其計算精度高于其他兩類模型[12]。
作為遙感蒸散發(fā)模型的主要輸入信息,遙感數(shù)據(jù)對模型的估算精度有重要影響。目前,常用于估算區(qū)域作物蒸散發(fā)的遙感數(shù)據(jù)主要有MODIS(Moderate resolution imaging spectroradiometer)數(shù)據(jù)[13-14]、Landsat ETM+數(shù)據(jù)[15]、Sentinel-2數(shù)據(jù)[16]等。國內(nèi)灌區(qū)常見地形、田塊和種植結(jié)構(gòu)復(fù)雜,往往需要高時間分辨率和高空間分辨率的遙感數(shù)據(jù),以實現(xiàn)農(nóng)田作物蒸散發(fā)的精準(zhǔn)、連續(xù)監(jiān)測和反演[17]。盡管可以利用一些數(shù)據(jù)融合算法調(diào)和MODIS數(shù)據(jù)和Landsat ETM+數(shù)據(jù)在空間分辨率和時間分辨率上的矛盾[18],但復(fù)雜的算法和步驟限制了其大范圍實際應(yīng)用和推廣。Sentinel-2數(shù)據(jù)具備高時間分辨率(5 d)和高空間分辨率(10~60 m),彌補了MODIS數(shù)據(jù)和Landsat ETM+數(shù)據(jù)在“時、空”分辨率上的缺陷。近年來,Sentinel-2數(shù)據(jù)在作物識別、植被指數(shù)反演等方面得到了廣泛應(yīng)用[19-22]。但是基于Sentinel-2數(shù)據(jù)對作物蒸散發(fā)量進行估算的研究還較為鮮見。
地表溫度(Land surface temperature, LST)作為雙源模型中區(qū)分植被溫度和土壤溫度的關(guān)鍵因子,是作物需水狀況和土壤水分信息的綜合反映。目前,對于遙感地面溫度數(shù)據(jù)同樣存在著“時、空”分辨率矛盾的問題,因此,高時空分辨率的LST數(shù)據(jù)對提高遙感作物蒸散發(fā)量估算精度具有重要價值?;谏鲜鰡栴},本研究利用混合雙源蒸散發(fā)模型(HTEM),結(jié)合Sentinel-2數(shù)據(jù)和地面農(nóng)田實際觀測數(shù)據(jù),對寧夏回族自治區(qū)中衛(wèi)市玉米主要生育期(5—8月)的實際蒸散發(fā)量進行計算,以探究利用Sentinel-2數(shù)據(jù)和地表作物冠層溫度實測值進行區(qū)域作物實際蒸散發(fā)量估算的可行性,為區(qū)域玉米精量灌溉決策提供理論基礎(chǔ)。
試驗區(qū)位于寧夏回族自治區(qū)中衛(wèi)市沙坡頭區(qū)常樂鎮(zhèn)農(nóng)業(yè)示范區(qū)(105°7′6.76″E, 37°28′5.7″N)(圖1)。研究區(qū)屬于溫帶大陸性季風(fēng)氣候,海拔1 225 m。受沙漠影響,日照充足,晝夜溫差大,年內(nèi)平均氣溫為7.3~9.5℃,年均降水量180~367 mm,年蒸發(fā)量為1 829 mm,夏季降水量占全年降水量的61%,全年日照時長2 800 h左右[23]。當(dāng)?shù)赝寥乐饕獮楣嘤偻痢Q芯繀^(qū)種植的農(nóng)作物以玉米為主,土地面積約為566.95 hm2。于2019年在玉米的主要生育期(5—8月)進行試驗觀測,2019年5月初播種,9月中旬收青貯玉米。研究區(qū)內(nèi)自然降雨時段與玉米整個生育期吻合;玉米田塊共灌溉5次,總灌溉量為380.81 mm,灌溉日期及灌水量如圖2所示。
根據(jù)研究區(qū)2019年玉米分布情況,選取2個代表性地塊,在每個地塊中間布置1套CTMS-On line型作物冠層溫度及環(huán)境因子測量系統(tǒng)(分別為H1和H2,圖1)。觀測系統(tǒng)田間實物如圖3所示,數(shù)據(jù)采集時間間隔均為30 min。本系統(tǒng)同步連續(xù)監(jiān)測的農(nóng)田信息包括:空氣溫度/濕度、風(fēng)速、太陽輻射、光合有效輻射、大氣壓強、作物冠層紅外溫度、根區(qū)土壤墑情等[24]。每7~10 d人工進行觀測和采集玉米株高、干物質(zhì)量、葉面積等作物信息。地下水位觀測值通過寧夏水利科學(xué)研究院在示范區(qū)安裝的觀測井獲得。其他氣象數(shù)據(jù)(如降水量等)從國家氣象科學(xué)數(shù)據(jù)中心下載得到。
Sentienl-2數(shù)據(jù)來源包括Sentienl-2A和Sentienl-2B衛(wèi)星,各波段信息如表1所示。空間分辨率為10~60 m,時間分辨率為5 d。Sentinel-2數(shù)據(jù)下載及處理流程主要包括原始數(shù)據(jù)下載(https:∥scihub.copernicus.eu)和校正、輻射定標(biāo)和裁剪、反演相關(guān)參數(shù)等3個關(guān)鍵步驟。數(shù)據(jù)下載完成后利用SNAP平臺對原始數(shù)據(jù)進行處理,包括幾何校正、輻射定標(biāo)、重采樣和裁剪等處理。并進一步反演獲得遙感蒸散發(fā)模型所需的反照率、歸一化植被指數(shù)(NDVI)、植被覆蓋度(fc)等關(guān)鍵參數(shù)。通過去除生育期云遮蓋等的影響,共選取10景數(shù)據(jù)進行研究,如表2所示。
表1 Sentinel-2數(shù)據(jù)各波段信息Tab.1 Statistical band information of Sentinel-2 data
表2 選用的Sentinel-2數(shù)據(jù)Tab.2 Selections of Sentinel-2 remote sensed data
基于雙源模型的原理,忽略光合作用耗能和水平方向能量交換的地表能量平衡方程為
Rn=H+LE+G
(1)
式中Rn——凈輻射通量,W/m2
H——感熱通量,W/m2
LE——潛熱通量,W/m2
G——土壤熱通量,W/m2
潛熱通量(LE)作為能量余項,可通過求得地表能量平衡方程中的Rn、H和G求得。
凈輻射通量采用ALLEN等[25]提出的基于遙感信息的地表凈輻射計算方法獲得。土壤熱通量采用BASTIAANSSEN[26]提出的半經(jīng)驗?zāi)P瓦M行估算,其表達式為
G=Rn(Lst-273.16)(0.003 8+0.007 4α)·
(1-0.98NDVI4)
(2)
其中
α=0.356ρ2+0.130ρ4+0.373ρ8+
0.085ρ11+0.072ρ12-0.001 8
(3)
(4)
式中Lst——地表溫度,K
α——地表寬波段反照率[27]
NDVI——歸一化植被指數(shù)
其中,ρi為各窄波段反射率。地表寬波段反照率需用Band 2、Band 4、Bnad 8、Band 11和Band 12波段的反射率;歸一化植被指數(shù)需用Band 4和Band 8波段的反射率。
HTEM模型采用層狀模型對凈輻射在植被和土壤組分間進行分配,采用塊狀模型計算地表的顯熱通量、潛熱通量及土壤熱通量。植被和土壤潛熱通量計算式為
(5)
(6)
式中LEc、LEs——植被、土壤潛熱通量,W/m2
Rnc、Rns——植被、土壤凈輻射,W/m2,可通過消光系數(shù)計算得到
ρair——空氣密度,kg/m3
Cp——空氣比熱容,J/(kg·K)
Tc、Ts——植被、土壤表面溫度,K,可由植被指數(shù)-地表溫度梯形特征空間來確定
Ta——空氣溫度,K
參數(shù)具體計算方法和步驟見文獻[10]。
由遙感蒸散發(fā)模型計算得到的是衛(wèi)星過境瞬時的潛熱通量,需要通過一定的方法將瞬時潛熱通量進行時間尺度擴展來獲得日蒸散發(fā)量及更長時段內(nèi)的蒸散發(fā)量。采用修正蒸發(fā)比法對瞬時潛熱通量進行時間尺度擴展得到衛(wèi)星過境日的蒸散發(fā)量,計算公式為
EF=λETi/(Rn-G)i
(7)
λETd=mEF(Rn-G)d
(8)
式中EF——蒸發(fā)比
m——修正系數(shù),取1.1
λ——水的汽化潛熱,取2.45 MJ/kg
下標(biāo)i、d表示瞬時值和日值。
對于缺乏遙感數(shù)據(jù)的時段,蒸散發(fā)量采用日參考蒸發(fā)比三次樣條插值法得到。
采用水量平衡方法對遙感反演作物蒸散發(fā)量進行驗證和評價,水量平衡公式為
ET=Pr+U+I-D-R-ΔW
(9)
式中ET——作物蒸散發(fā)量,mm
Pr——有效降雨量,mm
U——地下水補給量,mm
I——灌水量,mm
D——深層滲漏量,mm
R——徑流量,mm
ΔW——試驗初期和試驗?zāi)┢谕寥纼λ康淖兓?,mm
因試驗區(qū)屬于干旱半干旱地區(qū)且地勢平坦,R和D可忽略不計,則水量平衡公式可簡化為
ET=Pr+U+I-ΔW
(10)
有效降雨量采用有效降水系數(shù)和實際降水量的乘積得到,試驗區(qū)5—8月有效降水系數(shù)分別為69.87%、78.17%、81.15%、88.16%[26]。地下水補給量采用地下水給水度和研究時段內(nèi)地下水位變化量乘積得到,本研究中地下水給水度取0.045[28]。土壤儲水量的變化量采用試驗?zāi)┢诤驮囼灣跗谕寥纼λ坎钪登蟮茫渲型寥纼λ坑嬎愎綖?/p>
(11)
式中W——土壤儲水量,mm
θi——第i層土壤體積含水率,%
hi——第i層土壤厚度,cm
n——土層的層數(shù),本試驗設(shè)0~20 cm、20~40 cm和40~60 cm共計3個土層
選用相對誤差(Relative error,RE)、絕對誤差(Absolute error,AE)和均方根誤差(Root mean square error,RMSE)作為精度評價指標(biāo),其中相對誤差、絕對誤差和均方根誤差越小,說明估算的精度越高。采用OriginPro 9.1軟件制圖,統(tǒng)計參數(shù)分析在Excel 2017中進行。
遙感蒸散發(fā)計算的關(guān)鍵輸入數(shù)據(jù)包括寬波段反照率(α)、歸一化植被指數(shù)(NDVI)、株高、地表冠層溫度、氣象數(shù)據(jù)等。其中α和NDVI通過遙感數(shù)據(jù)反演得到,株高通過NDVI進行估算;地表溫度(玉米冠層溫度)和氣象數(shù)據(jù)通過田間實測獲得。
2.1.1反照率和NDVI反演結(jié)果
Sentinel-2數(shù)據(jù)反演研究區(qū)反照率和NDVI結(jié)果如圖4所示。以2019年6月24日遙感數(shù)據(jù)為例,研究區(qū)農(nóng)田寬波段反照率主要集中在0.1~0.2之間,農(nóng)田NDVI主要集中在0.5~0.8之間,其中玉米觀測站點南部有部分田塊NDVI明顯低于周圍田塊,主要是由于2019年該部分田塊種植作物為枸杞。枸杞屬灌木,NDVI小于玉米等農(nóng)作物,反演結(jié)果空間分布與實際種植情況相吻合。另外,從反照率和NDVI反演結(jié)果可以看出,Sentinel-2數(shù)據(jù)具有較高的空間分辨率(反照率空間分辨率20 m,NDVI空間分辨率10 m),適合用于研究區(qū)農(nóng)田地塊復(fù)雜、種植規(guī)模小的情況,降低了混合像元的存在,提高了反演精度。
2.1.2玉米NDVI時間序列擬合
研究區(qū)玉米試驗田塊NDVI時間序列擬合曲線選取DoseResp“S”形曲線,迭代算法采用Levenberg-Marquardt優(yōu)化算法。DoseResp曲線表達式為
(12)
式中Z——儒略日(DOY)
a、b、c、e——待擬合參數(shù),其中a為曲線下漸近線,b為曲線上漸近線
對研究區(qū)玉米NDVI時間序列進行擬合,結(jié)果如圖5和表3所示。可以看出,“S”形曲線能夠很好地描述研究區(qū)玉米NDVI時間序列,決定系數(shù)R2達到0.99以上。從玉米NDVI時間序列變化曲線來看,玉米播種后至出苗前期,NDVI約為0.2,表現(xiàn)為與裸土NDVI接近。隨著玉米生育期的增加,在第160天到第180天,NDVI由0.2迅速增至0.8,這一階段為玉米快速生長階段。此后至收獲青貯玉米之前,NDVI一直保持在0.8左右。由擬合參數(shù)可以看出,a和b表示擬合曲線理論上的最大值和最小值,玉米NDVI時間序列理論最大值和最小值分別為0.797和0.183。
表3 2019年研究區(qū)玉米NDVI時間序列曲線擬合結(jié)果Tab.3 Maize NDVI series of study area and logistic curve fitting results in 2019
2.1.3玉米株高與NDVI關(guān)系
遙感蒸散發(fā)模型中,空氣動力學(xué)阻抗的計算需要冠層高度作為輸入數(shù)據(jù),2019年共進行10次玉米株高觀測,分別為5月19日、6月4日、6月14日、6月24日、7月4日、7月18日、7月30日、8月9日、8月19日和8月29日。玉米株高與NDVI關(guān)系如圖6所示,并利用指數(shù)曲線對株高與NDVI關(guān)系進行擬合。由擬合結(jié)果可以看出,玉米株高與NDVI呈指數(shù)關(guān)系,決定系數(shù)為0.94。
2.2.1遙感凈輻射和蒸散發(fā)計算結(jié)果
2019年研究區(qū)玉米站點(H1和H2)的凈輻射觀測值和衛(wèi)星過境瞬時凈輻射計算結(jié)果如圖7所示??梢姡芯繒r段內(nèi)衛(wèi)星過境瞬時凈輻射變化范圍為590~710 W/m2,利用遙感信息獲取地表凈輻射具有較高的精度,計算值與觀測值決定系數(shù)為0.786,均方根誤差(RMSE)為36.256 W/m2。
利用HTEM模型計算得到的研究區(qū)內(nèi)2個玉米試驗田塊瞬時潛熱通量和采用修正的蒸發(fā)比法得到的日蒸散發(fā)量結(jié)果如表4(E為土壤蒸發(fā)量,T為植被蒸騰量)和圖8所示??梢?,2個觀測站點(H1和H2)在玉米生育期內(nèi)潛熱通量和日蒸散發(fā)量變化趨勢相似,生育初期凈輻射量約為600 W/m2,潛熱通量約為350 W/m2,潛熱通量主要為土壤潛熱,占總潛熱通量的95%以上。時間尺度擴展得到的日蒸散發(fā)量變化范圍為3.742~4.676 mm/d,水分消耗主要以土壤蒸發(fā)為主,范圍為3.470~4.303 mm/d,植被蒸騰很小,變化范圍為0.112~0.373 mm/d。在玉米快速生長階段,凈輻射逐漸增加,約為700 W/m2,潛熱通量也隨之增加,變化范圍為495.481~613.873 W/m2。日實際蒸散發(fā)量在這一階段也有明顯增加,變化范圍為5.573~5.999 mm/d,主要的水分消耗形式從土壤蒸發(fā)逐漸變?yōu)橹脖徽趄v,在第175天時,植被蒸騰量占日蒸散發(fā)量的95%以上,并且這一比例一直保持到青貯玉米收獲。玉米日實際蒸散發(fā)量在第185天達到峰值,2個玉米觀測站點的日實際蒸散發(fā)量分別為6.311 mm/d和6.547 mm/d。隨后日實際蒸散發(fā)量逐漸減小至4.106 mm/d和4.617 mm/d。
采用日參考蒸發(fā)比插值的方法獲取無衛(wèi)星過境日的蒸散發(fā)量,進而得到生育期內(nèi)連續(xù)的玉米日實際蒸散發(fā)量,采用三次樣條插值,并對插值后的參考蒸發(fā)比異常值進行修正,得到玉米生育期內(nèi)連續(xù)的日實際蒸散發(fā)量(圖9)。結(jié)果表明,逐日實際蒸散發(fā)量同樣呈現(xiàn)出單峰的變化趨勢,即生育初期日實際蒸散發(fā)量在4 mm/d左右浮動。隨著生育期的增加,日實際蒸散發(fā)量逐漸增大,在第180天到第220天達到峰值,期間個別時段日實際蒸散發(fā)量超過8 mm/d,這是由于這一階段降水較為集中,其中年內(nèi)兩次較大的降水均發(fā)生在該時段,并且進行了兩次灌溉,土壤含水率較高,水分脅迫較小,加之這一時段作物潛在蒸散發(fā)量較大,導(dǎo)致這一階段實際蒸散發(fā)量達到了峰值。至生育期末,日實際蒸散發(fā)量逐漸減小為2~4 mm。從統(tǒng)計結(jié)果可以看出(表5),研究區(qū)2個玉米試驗田塊5—8月實際蒸散發(fā)總量分別為525.114 mm和533.690 mm,日均實際蒸散發(fā)量為4.269 mm/d和4.339 mm/d,其中植被蒸騰總量分別為363.483 mm和358.196 mm,從整個生育期來看,植被蒸騰是作物水分消耗的主要形式。
表5 2019年研究區(qū)玉米實際蒸散發(fā)量Tab.5 Statistical results of maize evapotranspiration of study area in 2019
2.2.2ET計算結(jié)果驗證
采用水量平衡法對計算得到的玉米生育期內(nèi)蒸散發(fā)總量進行驗證和評價(表6)。結(jié)果表明,基于冠層溫度觀測系統(tǒng)和Sentinel-2數(shù)據(jù)反演的研究區(qū)玉米蒸散發(fā)總量與水量平衡法計算得到的蒸散發(fā)總量相比,精度較高,2個觀測站點的相對誤差分別為-2.512%和-1.469%,絕對誤差分別為13.533、7.774 mm??梢?,利用當(dāng)?shù)靥镩g多參數(shù)觀測系統(tǒng)和Sentinel-2數(shù)據(jù)并結(jié)合混合雙源蒸散發(fā)模型獲取作物蒸散發(fā)量比較可靠。
表6 2019年研究區(qū)玉米蒸散發(fā)量驗證結(jié)果Tab.6 Verification of maize evapotranspiration in study area in 2019
Sentinel-2數(shù)據(jù)在可見光、近紅外和短波紅外工作譜段范圍內(nèi)具有高時空分辨率,能夠很好地描述下墊面的高異質(zhì)性,大大降低了混合像元的存在。但是地表溫度作為遙感蒸散發(fā)模型的關(guān)鍵輸入數(shù)據(jù),目前仍存在著“時、空”上的矛盾,在很大程度上增加了蒸散發(fā)估算的不確定性。隨著遙感技術(shù)的發(fā)展以及地表溫度解譯方法的成熟,建立高時空分辨率的地表溫度算法和產(chǎn)品是提高區(qū)域蒸散發(fā)估算精度的重要課題。
在瞬時蒸散發(fā)到日蒸散發(fā)量的時間尺度擴展方面,不同方法均有其適用條件,如何根據(jù)不同的環(huán)境因素、衛(wèi)星過境時刻、下墊面類型等選取適宜的時間尺度擴展方法仍值得深入研究。此外,利用擴展得到的日蒸散發(fā)量獲取連續(xù)時段的蒸散發(fā)量的插值方法仍需進一步研究,特別是針對有降水前后由于云遮蔽的影響缺乏遙感數(shù)據(jù),同時又由于降水導(dǎo)致的土壤含水率突變的情況下,如何提高插值的精度也是下一步需繼續(xù)開展和完善的工作。
(1)Sentinel-2數(shù)據(jù)具有高時空分辨率,能夠與研究區(qū)復(fù)雜的種植地塊相匹配,降低了混合像元的存在,提高了玉米蒸散發(fā)量估算的精度。利用“S”形曲線擬合研究區(qū)玉米NDVI時間序列具有較高精度,決定系數(shù)達到0.99以上。同時,研究區(qū)玉米株高與NDVI的關(guān)系可以用指數(shù)曲線擬合,決定系數(shù)達到0.94。
(2)基于Sentinel-2數(shù)據(jù)反演得到的衛(wèi)星過境瞬時凈輻射具有較高精度,其變化范圍為590~710 W/m2,與地面觀測值間的決定系數(shù)為0.786,均方根誤差為36.256 W/m2。
(3)衛(wèi)星過境日實際蒸散發(fā)量呈單峰變化,在第185天左右達到峰值,日蒸散發(fā)量達到6 mm/d以上。生育初期水分主要以土壤蒸發(fā)形式消耗,從第150天之后逐漸過渡到以植被蒸騰為主,到第175天時,植被蒸騰量占日實際蒸散發(fā)量的95%以上,并且一直保持到青貯玉米收獲。由參考蒸發(fā)比三次樣條插值得到的逐日實際蒸散發(fā)量結(jié)果表明,研究區(qū)玉米5—8月實際蒸散發(fā)總量為529.402 mm,日均實際蒸散發(fā)量為4.304 mm/d,其中植被蒸騰總量為360.840 mm,從整個生育期來看,植被蒸騰是作物水分消耗的主要形式。采用水量平衡法對玉米生育期蒸散發(fā)總量進行驗證,兩個觀測站點的平均相對誤差和絕對誤差分別為-1.990 5%和10.653 5 mm。
(4)將地面冠層溫度觀測數(shù)據(jù)和高時空分辨率遙感數(shù)據(jù)相結(jié)合,采用混合雙源遙感蒸散發(fā)模型(HTEM)獲取研究區(qū)玉米生育期蒸散發(fā)量(蒸發(fā)量、蒸騰量)具有較高的精度,對精量灌溉決策的制定和區(qū)域農(nóng)業(yè)水管理具有重要的科學(xué)意義和應(yīng)用價值。