南京地區(qū)蒸散發(fā)降尺度研究
——基于增強型時空自適應反射融合模型

尉毓姣,朱 琳,*,曹鑫宇,王文科,龔建師,余慧琳,孟 丹

1 首都師范大學 資源環(huán)境與旅游學院,北京 100048 2 城市環(huán)境過程與數字模擬國家重點實驗室培育基地,北京 100048 3 首都師范大學 水資源安全北京實驗室,北京 100048 4 長安大學 水利與環(huán)境學院,西安 710054 5 中國地質調查局南京地質調查中心,南京 210016 6 自然資源部流域生態(tài)地質過程重點實驗室,南京 210016

蒸散發(fā)是水文循環(huán)中自降水到達地面后由液態(tài)或固態(tài)轉化為水汽返回大氣的過程,包括水面蒸發(fā)、土壤蒸發(fā)和植被蒸騰,是陸地水分循環(huán)和能量交換的重要組成部分,能夠調節(jié)全球能量和水量平衡。精確估算蒸散發(fā)以及獲取高時空分辨率的蒸散發(fā)數據對區(qū)域乃至全球氣候變化、水資源評價和陸地碳水循環(huán)具有重要意義[1—2]。

針對大范圍、快速變化的地表信息監(jiān)測,單一傳感器無法滿足需求,許多學者針對這一問題,提出了切實可行的遙感數據時空融合技術[3]。2006年,Gao等人[4]首次提出STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)融合算法,開創(chuàng)了遙感數據融合的先河。Yi等人[5]使用STARFM模型對MODIS和ASTER蒸散發(fā)結果進行數據融合,用于獲取農田尺度上連續(xù)的日蒸散發(fā),發(fā)現在非均質農業(yè)區(qū)融合精度低于均質區(qū)域。Yang等人[6]使用STARFM模型融合基于Landsat和MODIS數據估算的蒸散發(fā),發(fā)現精度與僅使用Landsat估算的蒸散發(fā)相似或略低。STARFM算法進行遙感數據的融合可以較好的捕捉植被的物候變化,適用于地表覆蓋同質的像元,而對于復雜、異質性的地表覆蓋,預測效果不佳。Zhu等人[7]提出增強型時空自適應反射融合模型(Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model,ESTARFM),該模型考慮了像元之間的光譜和空間相似性,保留了空間細節(jié),適合復雜地物區(qū)域的數據預測。程筱茜等人[8]利用ESTARFM模型模擬了內蒙古紅堿淖1987—2018年間缺失的Landsat影像數據,用于反映水體的變化情況,該模型模擬精度達到0.95,彌補了Landsat數據在時間分辨率上的不足。Zhou等人[9]利用ESTARFM模型對歸一化植被指數(Normalized Difference Vegetation Index,NDVI)和地表溫度(Land Surface Temperature,LST)分別進行融合,用于提高干旱監(jiān)測的精度,發(fā)現預測結果與實際數據的相關系數分別為0.9和0.74。Heimhuber等人[10]分別利用STARFM和ESTARFM模型對Landsat影像和MODIS影像進行了8 d間隔、30 m分辨率的數據預測,發(fā)現ESTARFM結果優(yōu)于STARFM,在漫灘地區(qū)融合精度最佳。鄔明權等人[3]將STARFM、ESTARFM以及統計回歸等模型進行對比,發(fā)現ESTARFM模型的融合精度最高。

目前,對于蒸散發(fā)的研究以中國北部地區(qū)居多,且研究內容以蒸散發(fā)估算、時空特征分析為主,而對于南部地區(qū)蒸散發(fā)降尺度的研究較少。因此,本文選擇南京市為案例區(qū),結合Landsat- 8遙感影像數據和氣象數據,采用基于能量平衡原理的SEBS(Surface Energy Balance System)模型估算南京市的地表日蒸散量。在此基礎上,選取南京市部分區(qū)域(400像元×400像元,面積144 km2),采用ESTARFM模型將高空間分辨率(30 m)的Landsat蒸散發(fā)估算結果與高時間分辨率(8 d)的MOD16A2產品數據進行時空融合,并評價模型的融合精度。研究結果可為南京市水資源的合理規(guī)劃提供基礎數據。

1 研究區(qū)概況

南京市(31°14′—32°37′N,118°22′—119°14′E)位于中國東部、長江下游中部地區(qū),是長三角輻射帶動中西部地區(qū)發(fā)展的國家重要門戶城市,總面積6587 km2。區(qū)內四季分明,雨水充沛,多年平均溫度為16.6 ℃(2000—2019年),年平均降雨量為1150.8 mm,屬亞熱帶季風濕潤氣候。長江穿城而過,沿江岸線總長近200 km,境內共有大小河道120條,水域面積占總面積的11%。南京市是中國重要的農業(yè)和商品糧基地之一,植物資源豐富、種類繁多。全市森林覆蓋率27%,建成區(qū)綠化覆蓋率45%,人均公共綠地面積13.7 m2。研究區(qū)示意圖如圖1所示。

圖1 研究區(qū)示意圖Fig.1 Location map of study area

2 數據與研究方法

2.1 數據來源與預處理

本文采用的數據包括遙感影像數據和專題數據,其中專題數據包括：氣象數據、土地利用類型數據、數字高程模型以及MOD16A2蒸散發(fā)產品數據。

2.1.1遙感影像數據

由于時空數據融合需要盡量選擇日期相同或相近的兩種數據源,且融合需要三個時期的數據,本次經過數據篩選后,發(fā)現Landsat- 8遙感影像在2017年7月21日、10月9日和12月12日數據質量較好(云量小于2%),且能夠與MOD16A2蒸散發(fā)產品保持日期一致。故選擇覆蓋南京市的2017年7月21日、10月9日和12月12日三期Landsat- 8遙感影像(地理空間數據云：http://www.gscloud.cn/),空間分辨率為30 m,時間分辨率為16 d。在ENVI 5.3中對影像數據進行輻射定標、大氣校正、圖像裁剪等預處理,通過波段運算工具計算歸一化植被指數NDVI、地表比輻射率和反照率參數,采用大氣校正法計算地表溫度[11]。

2.1.2專題數據

(1)氣象數據

氣象數據來源于中國氣象數據網(http://data.cma.cn/),選取南京市及周邊城市共計11個站點的日風速、氣溫、氣壓、相對濕度等氣象要素,作為SEBS模型的輸入數據。在ArcGIS 10.5中采用反距離加權插值方法將站點數據進行空間插值,得到30 m分辨率的柵格數據。

(2)土地利用類型數據

土地利用類型數據為2018年1 km空間分辨率的柵格數據,來源于中科院資源環(huán)境科學與數據中心網站(http://www.resdc.cn/Default.aspx)。原數據包括耕地、林地、草地、水域、居民地和未利用土地6個一級類型和25個二級類型。在ArcGIS 10.5中對二級類型進行合并處理,并將合并結果重分類為耕地、林地、水域、草地和其他五類(表1),通過裁剪得到南京市土地利用類型圖(圖2),其中其他類別主要包括建設用地和未利用土地。

表1 南京市土地利用類型分類表

圖2 研究區(qū)土地利用類型圖Fig.2 Land use type map of the study area

(3)數字高程模型(DEM)

數字高程模型(Digital Elevation Model,DEM)為30 m空間分辨率的SRTM產品,下載于Google Earth Engine平臺(https://earthengine.google.com/),并通過裁剪處理得到南京市DEM圖。

(4)MOD16A2蒸散發(fā)產品數據

MOD16A2蒸散發(fā)產品數據來源于美國國家航空航天局(NASA)網站(https://ladsweb.modaps.eosdis.nasa.gov),空間分辨率為500 m,時間分辨率為8 d。本次通過MRT(MODIS Reprojection Tool)工具對MOD16A2產品數據進行格式轉換、投影轉換等預處理,并重采樣至30 m格網大小,便于后續(xù)與基于Landsat-遙感影像反演的蒸散發(fā)結果進行融合計算。此外,該產品數據對建設用地、水域、裸地等區(qū)域賦值無意義數值,為無效值(https://lpdaac.usgs.gov/products/mod16a2v006/)。為了減少對融合結果的影響,本次采用鄰域均值代替無效值,并選取無效值較少的區(qū)域(400像元×400像元,面積144 km2)作為后續(xù)的融合計算區(qū)。

2.2 研究方法

2.2.1SEBS模型

2002年,荷蘭學者蘇中波[12]提出基于能量平衡原理的單層模型—SEBS模型,該模型具有較高的估算精度和更好的實用價值,目前廣泛應用于干旱監(jiān)測、水資源評價等方面[13—16]。該模型假設每個像元的顯熱通量都介于該像元的干限和濕限顯熱通量之間,一定程度上避免了氣象數據空間插值不確定性帶來的誤差[17—18]。SEBS模型的地表能量平衡方程[19]為：

Rn=G0+H+λE

(1)

式中,Rn為地表凈輻射通量(W/m2)；G0為土壤熱通量(W/m2)；H為感熱通量(W/m2)；λE為潛熱通量(W/m2),其中E為蒸散量(mm),λ為水的汽化潛熱(Wm-2mm-1),與地表溫度Ts(K)之間的關系式為[20]：

λ=[2.501-0.002361×(Ts-273.15)]×106

(2)

(1)地表凈輻射通量Rn

地表凈輻射(輻射平衡)指地表向上和向下的總輻射之差,是地面所獲取的凈輻射能量[21]。計算公式為：

(3)

ε=0.004×Pv+0.986

(4)

式中,Pv為植被覆蓋度,公式計算如下：

(5)

式中,NDVI為歸一化植被指數；NDVISoil為裸土或無植被覆蓋區(qū)的NDVI值；NDVIVeg為植被完全覆蓋區(qū)的NDVI值。本次取累計概率為5%和95%的NDVI值作為NDVISoil和NDVIVeg,對植被覆蓋度進行計算[23]。

(2)土壤熱通量G0

土壤熱通量是單位面積土壤在單位時間內的熱交換量[24],是由于傳導導致的儲存到土壤和植被中的熱量的比率[20],經驗統計公式為：

(6)

(3)感熱通量H

感熱通量也叫做顯熱通量,是指由于溫度變化而引起的大氣與下墊面之間發(fā)生的湍流形式的熱交換[25],是由地表溫度和參考高度的氣象參數獲取的,與溫度差值成正比[26],計算公式為：

H=ρCp(Ts-Ta)/γa

(7)

式中,ρ=1.29(kg/m3)為空氣密度；Cp=1004(Jkg-1K-1)為空氣定壓比熱容；Ta為參考高度的溫度(K)；γa為空氣動力學阻抗(s/m)。空氣動力學阻抗γa與風速、粗糙度和大氣層結構有關,根據大氣邊界層相似理論,存在以下關系：

(8)

(9)

(10)

(4)蒸發(fā)比

基于能量平衡公式計算的蒸散發(fā)是衛(wèi)星過境時刻的瞬時值,本文采用恒定蒸發(fā)比法將瞬時蒸散發(fā)擴展到日蒸散發(fā)[29]。該時間尺度擴展方法簡單易用,適用于晴朗或云量保持恒定不變的白天[30]。計算公式為：

(11)

(12)

式中,ET24為日實際蒸散發(fā)量(mm)；EF為蒸發(fā)比；Rn24為日凈輻射通量(W/m2)；G024為日土壤熱通量(W/m2)；ρw為水的密度。

2.2.2ESTARFM時空融合模型

ESTARFM時空融合模型是對STARFM模型的改進,其原理是考慮同一時間同一區(qū)域內像元空間和光譜的相似性,依據預測日期前后至少兩對高低分辨率影像,以及預測日期當天的一幅低分辨率影像,計算對應像元之間的權重和轉換系數,從而模擬出預測日期的高時空分辨率數據[31]。本次利用待預測時期(2017年10月9日)前后兩個時期(2017年7月21日、12月12日)基于Landsat-8遙感影像反演的蒸散發(fā)數據和MOD16A2產品數據以及待預測時期的MOD16A2產品數據,模擬出待預測時期30 m的蒸散發(fā)數據[32]。計算公式為[5]：

ET(xw/2,yw/2,tp)=TmETm(xw/2,yw/2,tp)+TnETn(xw/2,yw/2,tp)

(13)

(14)

式中,ET為最終預測時期的蒸散發(fā)；tp為預測時期；ETm為Tm時期預測的蒸散發(fā)；ETn為Tn時期預測的蒸散發(fā)；ETL為基于Landsat-8遙感影像反演的蒸散發(fā)；ETM為MODI6A2蒸散發(fā)；ETk為Tk時期預測的蒸散發(fā)；(xw/2,yw/2)為中心像元的位置；w為相似像元搜索窗口；(xi,yi)為第i個相似像元的位置；N′為相似像元的個數；Vi為轉換系數；Wi為綜合權重因子,計算公式如下：

(15)

Di=(1-Ri)di

(16)

(17)

式中,di為距離權重；Ri為光譜相似權重。

Tk為Tm、Tn時期的時間權重因子,計算公式為：

(18)

3 結果與分析

3.1 SEBS模型結果驗證

以往研究利用通量站數據、蒸發(fā)皿觀測數據或P-M公式計算的潛在蒸散發(fā)量與蒸散發(fā)模型估算的結果進行對比驗證,以證明結果的合理性[33—35]。因未收集到通量站實測數據,本文采用南京市氣象站點(31°55′N,118°54′E)在3個研究時期的蒸發(fā)皿折算數據以及MOD16A2產品數據進行SEBS模型30 m Landsat蒸散發(fā)估算結果的驗證。由于蒸發(fā)皿實測數據為水面蒸發(fā),本次參考相關文獻確定該研究區(qū)水面蒸發(fā)的折算系數[36—37],將水面蒸發(fā)轉化為實際蒸發(fā)后與SEBS模型結果進行對比。統計結果如表2所示,模型估算結果與蒸發(fā)皿折算數據的平均相對誤差為0.14 mm/d,與MOD16A2產品數據的平均相對誤差為0.22 mm/d。三者變化趨勢相同,并與已有研究[38—39]的蒸散發(fā)驗證結果的表現規(guī)律一致,可認為SEBS模型的估算結果合理。

表2 SEBS模型計算結果與蒸發(fā)皿折算數據、MOD16A2產品數據結果對比

3.2 SEBS模型估算的蒸散發(fā)時空特征

研究區(qū)3個不同時期的日蒸散量空間分布如圖3所示,7月21日平均蒸散量為7 mm,最大值為10 mm；10月9日平均蒸散量為3.9 mm,最大值為5.6 mm；12月12日平均蒸散量為1.1 mm,最大值僅為2.8 mm。由此得出,蒸散量在季節(jié)上差異較大,呈現夏季>秋季>冬季的變化規(guī)律。7月21日為植被的生長季中期,氣溫高,降水充沛,日照時間長,植被生長比較旺盛,尤其是農田水分灌溉充足,植被覆蓋度高,蒸騰量增加等都是導致該時間蒸散量較高的原因。10月9日為植被的生長季末期,氣溫降低,降水減少,氣候已發(fā)生明顯變化,并且農田水分灌溉減少,植被覆蓋度降低,因而蒸散量有所降低。12月12日正值冬季,氣溫較低從而導致蒸散量很低。

圖3 研究區(qū)日蒸散量空間分布圖Fig.3 Spatial distribution of daily evapotranspiration in the study area

在空間上,南京市各區(qū)在7月21日的平均蒸散量差異較大,最大為六合區(qū)的7.7 mm,最小為秦淮區(qū)的5.4 mm。六合區(qū)耕地面積626km2,水域面積202km2,植被蒸騰量和水面蒸發(fā)量在夏季都比較大,所以該地區(qū)蒸散量最大。秦淮區(qū)屬于城市中心區(qū)域,只有少部分自然河、人工河錯落和極少的城市綠地分布,蒸散量最小。

為了進一步分析研究區(qū)不同土地利用類型的蒸散發(fā)差異,選取特征最為明顯的7月21日的蒸散發(fā)結果,將土地利用與Landsat蒸散發(fā)數據進行疊加并分區(qū)統計,得到研究區(qū)不同土地利用類型的蒸散量。由表3可知,南京市不同土地利用類型的蒸散量大小為：水域>林地>耕地>草地>其他。水域是蒸散量最高的地物類型,因其水量充足,水面蒸發(fā)量最大；林地具有涵養(yǎng)水源的能力及較高的植被蒸騰,雖然面積占整個南京市面積的比例不大,但蒸散量較高；耕地面積最大,約占南京市總面積的51%,并且南部高淳區(qū)分布大量水田,夏季較高的氣溫、充足的灌溉條件等使得耕地蒸散量較高。

表3 研究區(qū)7月21日各土地利用類型平均蒸散量統計

3.3 ESTARFM模型結果及精度評價

本次輸入ESTARFM模型的初始數據為2017年7月21日、12月12日的兩組Landsat、MOD16A2數據以及預測時期10月9日的MOD16A2數據,以此融合預測10月9日的Landsat蒸散發(fā)數據,融合區(qū)域范圍見圖1。

圖4為融合區(qū)10月9日基于ESTARFM模型的融合結果與基于Landsat-8遙感影像反演的蒸散發(fā)結果,圖中亮白色代表蒸散發(fā)較高,為植被覆蓋度較高的區(qū)域；黑灰色代表蒸散發(fā)較低,主要為其他類型的地物。經統計,融合結果最大值為5.8 mm,最小值為0 mm,平均值為4 mm；Landsat蒸散發(fā)結果最大值為5.5 mm,最小值為0 mm,平均值為3.6 mm。整體上看,融合結果數值偏高,但二者在空間分布上具有一定的相似性。為了進一步度量時空融合模型的預測結果精度,本次選擇融合區(qū)域內均勻分布的10000個驗證點,并計算二者的相關系數。圖中(圖5)散點分布在1：1線附近,相關系數為0.74,表明基于ESTARFM模型融合的高分辨率蒸散發(fā)結果可靠。

圖4 融合結果與Landsat蒸散發(fā)對比Fig.4 Comparison of fusion result and Landsat evapotranspiration

圖5 Landsat蒸散發(fā)與融合結果相關性分析 Fig.5 Correlation analysis between Landsat evapotranspiration and fusion result

4 討論與結論

4.1 討論

4.1.1SEBS模型的不確定性分析

本文利用Landsat- 8遙感影像數據進行地表參數的反演,采用SEBS模型估算地表日蒸散量。研究結果經過多種方法的驗證,表明SEBS模型估算的蒸散發(fā)具有合理性。但是模型本身會受到輸入氣象要素與地表參數精度的影響,帶來模型反演結果的不確定性。因此,保證模型輸入參數的準確性對于提高模型的精度極其重要。同時,受地面通量站實測數據的限制,精度驗證的結果受到了一定的制約,獲取充足的蒸散發(fā)實測驗證數據是需要進一步解決的問題。此外,基于地表能量平衡的蒸散發(fā)模型分為單層模型和雙層模型,其中單層模型又稱為大葉模型,把地表視為一張大葉與外界進行水分和能量的交換,在下墊面均勻的陸面反演蒸散發(fā)的精度較高。本文采用的SEBS模型是單層模型的典型代表,雖然應用較為廣泛,但實際上研究區(qū)下墊面并非單一均勻組成,這在一定程度上也會對反演結果造成誤差。而雙層模型是單層模型的延伸,能夠定量描述陸面非均勻性對地表通量的影響,分別獲得土壤蒸發(fā)和植被蒸騰[40—41]。故在后續(xù)研究中將考慮下墊面異質性問題,采用特征空間的雙層蒸散模型進行區(qū)域蒸散發(fā)反演方面的研究。

4.1.2ESTARFM模型的不確定性分析

本文采用ESTARFM模型進行蒸散發(fā)產品數據的時空融合,研究結果表明,該模型可以有效地進行數據時空降尺度。為了進一步討論融合結果的質量,將融合區(qū)土地利用分別與ESTARFM模型融合結果和Landsat蒸散發(fā)估算結果進行疊加并分區(qū)統計,分析二者在土地利用類型上的蒸散發(fā)差異。融合區(qū)不同土地利用類型在10月9日的平均蒸散量統計結果如表4所示,由于融合區(qū)域范圍小,沒有涉及草地這一地物類型,故沒有做相應的統計。由表可以看出,二者均呈現水域>林地>耕地>其他的特點。ESTARFM模型融合結果與基于Landsat-8遙感影像反演的蒸散發(fā)結果的土地利用蒸散量差異在0.2—0.4 mm之間,其中,差異最大的土地利用類型為其他類型和水域,這可能與MOD16A2產品本身在建設用地、水域、裸地等區(qū)域存在無效值有關。由此可見,輸入ESTARFM模型的遙感數據質量對融合結果的質量會產生直接影響。

本次采用的Landsat- 8遙感影像為2017年7月21日、10月9日和12月12日三期,長時間間隔對融合結果會產生一定的影響。此外,融合結果會受到模型本身參數的影響,導致相似像元的選取會產生一定的誤差,使模型存在一定的不確定性。本次只針對小部分區(qū)域進行了時空數據融合降尺度研究,在后續(xù)的工作中可以嘗試將該模型應用到更大的尺度上,并且盡量選擇同一季節(jié)、質量更好的產品數據進行融合,為獲取長時間序列的蒸散發(fā)數據集提供數據保證,這也是未來進一步研究的方向。

表4 融合區(qū)10月9日土地利用類型平均蒸散量統計

4.2 結論

本文基于Landsat- 8遙感影像數據和氣象數據采用SEBS模型進行南京市地表日蒸散發(fā)的估算,分析蒸散發(fā)的時空分布特征及不同地物類型的差異。在此基礎上,采用ESTARFM模型融合Landsat蒸散發(fā)數據與MOD16A2蒸散發(fā)產品數據,進行時空降尺度研究。主要結論如下：

(1)采用蒸發(fā)皿折算后的數據與MOD16A2產品數據進行SEBS模型估算結果的精度評價,平均相對誤差分別為0.14 mm/d和0.22 mm/d,表明使用該模型進行該研究區(qū)蒸散發(fā)估算的結果合理。

(2)南京市蒸散量季節(jié)差異明顯,夏季蒸散量最大,冬季最小；各區(qū)在夏季的日平均蒸散量差異較大,六合區(qū)蒸散量最大,秦淮區(qū)最小；從土地利用類型上看,日平均蒸散量差異也較大,呈現水域>林地>耕地>草地>其他的特點。

(3)基于ESTARFM模型融合的蒸散發(fā)與基于Landsat-8遙感影像反演的蒸散發(fā)相比,數值偏高,但二者在空間分布上具有相似性,相關系數為0.74,融合結果可靠。對于不同的土地利用類型,受產品數據本身無效值的影響,融合效果差異最大的為其他類型和水域。