張 亮，崔林林，苑 躍

(1.高原與盆地暴雨旱澇災害四川省重點實驗室，成都 610072；2.四川省氣象災害防御技術中心，成都 610072； 3.成都信息工程大學資源環(huán)境學院，成都 610225)

前 言

地表蒸散(Evapotranspiration, ET)是全球水循環(huán)的重要組成部分，蒸散量每年約占陸地降水量的60%[1]，并在區(qū)域水平衡中發(fā)揮著重要作用。定量研究不同時空尺度下蒸散的分布，對于全球氣候變化、生態(tài)環(huán)境保護、水資源管理和農(nóng)業(yè)生產(chǎn)都具有重要意義。

目前，蒸散的估算方法主要分為傳統(tǒng)法和基于衛(wèi)星遙感的方法。傳統(tǒng)方法主要有蒸滲儀法、波文比法-能量平衡法、空氣動力學法、渦度相關儀法等，這些方法能夠在均勻的地表條件下提供較為精確的計算結(jié)果，但存在著理論模型復雜并且只能應用于小尺度區(qū)域等缺陷[2-3]。近年來隨著衛(wèi)星遙感技術的發(fā)展，基于衛(wèi)星遙感數(shù)據(jù)的蒸散估算方法可對大范圍區(qū)域的蒸散量進行計算[4～6]?；诘乇砟芰科胶饫碚摰腟EBAL(Surface Energy Balance Algorithm for Land)模型是當前應用較多的模型之一，在蒸散的計算過程中具有以下優(yōu)勢：(1)最小程度的使用地面觀測數(shù)據(jù)；(2)模型中近地表的空氣溫度的非強制性；(3)可通過識別干濕像元以及測定近地表空氣溫度差進行自校正[7～9]。研究表明，SEBAL模型蒸散計算結(jié)果在日尺度、季節(jié)尺度和年尺度上都具有較高的精度[9]。近年來國內(nèi)外學者利用SEBAL模型對不同生態(tài)系統(tǒng)的地表蒸散量進行了大量研究，如：Opoku-Duah et al.利用ASTER和MODIS數(shù)據(jù)對西非熱帶草原地區(qū)的蒸散量進行估算，并將結(jié)果與Landsat ETM+圖像中提取的ET進行了比較[10]；Teixeira et al.對SEBAL模型進行了進一步改進，并對巴西東北部舊金山河流域的蒸散量進行了估算[11]；程明瀚等對北京市日蒸散區(qū)域分布研究后發(fā)現(xiàn)，北京市蒸散整體呈現(xiàn)出西北山區(qū)高、東南平原低的趨勢[12]；周妍妍等對疏勒河流域蒸散量的時空變化研究后發(fā)現(xiàn)，且其空間分布異質(zhì)性十分明顯，整體上從東南向西北呈逐漸減小趨勢[13]；李寶富等估算了塔里木河干流區(qū)蒸散量，結(jié)果表明該區(qū)日蒸散量較高，其中水體蒸散最高，耕地和林地次之，居民和工業(yè)用地蒸散量最小[14]。

若爾蓋高原地處青藏高原東北部，區(qū)內(nèi)河網(wǎng)密集，是黃河上游重要水源補給區(qū)，被譽為“中國西部高原之腎”，同時也是全球氣候影響的敏感區(qū)。近幾十年，由于受人口增長、氣候變化、過度放牧等影響，該區(qū)域的土地利用/覆被發(fā)生劇烈變化，沙丘面積逐步擴大，進而導致區(qū)域水文過程(蒸散、徑流、截留等)和水循環(huán)過程的改變[15～17]。前人主要利用彭曼公式(Penman-Monteith formula)和MODIS的蒸散產(chǎn)品對若爾蓋高原的潛在蒸散量進行計算，但得到的結(jié)果與站點觀測值相差較大[18～20]。因此，本文利用SEBAL模型對近20年來若爾蓋高原的實際蒸散量進行估算，并對其時空變化特征進行分析，研究結(jié)果可為本地區(qū)的氣候變化、水資源管理和生態(tài)環(huán)境保護提供科學支撐。

1 研究地區(qū)與方法

1.1 研究區(qū)概況

若爾蓋高原位于青藏高原東北部邊緣，海拔介于3 400～3 900 m 之間，隸屬四川省若爾蓋縣、紅原縣和阿壩縣，以及甘肅省的瑪曲縣和碌曲縣(圖1)。若爾蓋高原為大陸性高原氣候，寒冷濕潤，長冬無夏，霜凍期長，日溫差大。研究區(qū)多年平均氣溫1℃左右，年日照時數(shù)均在2 000h以上，全年降水量 600～800mm，雨季分布在5～10月，約占全年降水量的 90%。本區(qū)植被類型以高山草甸和沼澤植被為主，區(qū)內(nèi)河流多屬于黃河水系，是黃河上游重要的水源補給區(qū)。

圖1 研究區(qū)域位置圖Fig.1 Location of research area

1.2 數(shù)據(jù)及方法

本研究所用若爾蓋高原及周邊299個國家及區(qū)域自動氣象站逐日數(shù)據(jù)，來源于四川省氣象局，包括氣溫(℃) 、平均風速(m/s) 、日照時數(shù)(h)、相對濕度(%) 和蒸散(mm)等，研究時段為2001～2017年，所有數(shù)據(jù)均經(jīng)過質(zhì)量控制，并通過回歸插值法對缺測數(shù)據(jù)進行補充。遙感數(shù)據(jù)為2001～2017年EOS/MODIS數(shù)據(jù)，主要包括1 km分辨率的地表溫度產(chǎn)品(MOD11A1)、1 km分辨率的植被指數(shù)產(chǎn)品(MOD13A2)、1 km分辨率的土地利用產(chǎn)品(MCD12Q1)、500 m分辨率的反照率產(chǎn)品(MCD 43A3)和1 km分辨率的反射率產(chǎn)品(MOD09GA)等，下載于美國國家航空和航天局的地球觀測系統(tǒng)數(shù)據(jù)和信息系統(tǒng)(http://reverb.echo.nasa.gov/reverb/)。所有的輸入數(shù)據(jù)均預處理為Albers等面積投影和1 km空間分辨率的數(shù)據(jù)集。

區(qū)域蒸散采用SEBAL模型進行計算，SEBAL模型基于地表能量平衡方程，通過衛(wèi)星遙感與氣象數(shù)據(jù)分別對各能量分量進行估算，最后得到區(qū)域的實際蒸散值。地表能量平衡方程為：

λET=Rn-G-H

(1)

公式(1)中λ為潛熱蒸發(fā)系數(shù)，取值2.45×106 MJ/kg，λET是潛熱通量(W/m2)，Rn是地表凈輻射通量(W/m2)， G是土壤熱通量(W/m2)，H是顯熱通量(W/m2)，Rn、G和H的計算方法介紹如下。

地表凈輻射通量Rn根據(jù)公式(2)進行計算。

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

(2)

公式(2)中： RS↓是下行到達地表的太陽短波輻射(W/m2)； RL↓是下行的長波輻射(W/m2)； RL↑是上行的長波輻射(W/m2)；α為地表反照率，εα為地表比輻射率，具體的參數(shù)計算過程如下。

RS↓=GSCcosθdrτSW

(3)

(4)

(5)

(6)

ε0=1.0094+0.0047lnNVDI

(7)

τSW=0.75+2×10-5z

(8)

公式(3)～(8)中: GSC為太陽常數(shù)，1 367 W/m2;θ為太陽天頂角; dr為日地距離因子，可通過儒略日求得; τsw為大氣單向透射率，根據(jù)研究區(qū)內(nèi)氣象觀測站的海拔獲得; z為高程; σ為5.67 × 10-8W/(m2·K4)；Ts為地表溫度(K); Ts_cold為“冷”點地表溫度(K)。α為地表反射率; αtoa為大氣層頂反射率; αpath_radiance為大氣層輻射值，NVDI為歸一化植被指數(shù)。

土壤熱通量G根據(jù)公式(9)進行計算。

(9)

顯熱通量H根據(jù)公式(10)進行計算。

(10)

(11)

(12)

公式(10)～(12)中：ρ為空氣密度(kg/m3)；Cp為空氣定壓比熱，取值1004J/(kg·K)；dT為地面Z1、Z2兩個高度的溫差，Z1和Z2是參考高度植被的零平面位移；rah為空氣動力學阻抗(s/m)。u*是摩擦速度(m/s); k是卡門常數(shù)，取0.4; u200是200 m 的摩擦風速(m/s)；zom是地表動量傳輸?shù)拇植诙乳L度(m)。其中，dT與rah為相互關聯(lián)的2個未知量，直接求解比較困難，SEBAL 模型引入Monin-Obukhov循環(huán)迭代計算，直至感熱通量穩(wěn)定。

由于公式(1)～(12)計算的是衛(wèi)星過境瞬時地面蒸散值(Λinst)。研究表明，一天內(nèi)蒸發(fā)比(Λ24)隨時間不發(fā)生變化，即Λ24=Λinst。因此利用蒸發(fā)比的這一特性對蒸散進行時間尺度的擴展，計算出一天的總蒸散量(式13和式14)。

(13)

ET24=86400Λ24(Rn24-G24)/λ

(14)

公式(13)(14)中：ET24為日蒸散量，Rn24和G24分別為日凈輻射通量和日土壤熱通量，其中日土壤熱通量G24由于白天和夜晚保持能量收支平衡，可以忽略不計。日凈輻射通量Rn24計算公式如下。

Rn 24= (1-α)Ra 24τSW- 110τSW

(15)

(16)

ωs=arccos(-tanφtanδ)

(17)

公式(15)～(17)中：φ為緯度，δ為太陽赤緯。

2 結(jié)果分析

2.1 SEBAL模型估算結(jié)果的對比驗證

本文選取2015年的若爾蓋、馬爾康、紅原、阿壩4個站點的蒸散實測數(shù)據(jù)(E-601B型蒸發(fā)器)與SEBAL模型估算出的蒸散量進行對比驗證。如表1所示，SEBAL模型的估算結(jié)果與四個站點的實測值的相關性較高，相關系數(shù)在0.53～0.75之間，相對誤差為25.05%～33.53%之間。由此可見，SEBAL模型在本研究區(qū)有較好的適用性。

表1 站點蒸散實測值與SEBAL模型估算結(jié)果對比Tab.1 The results comparation between SEBAL model and measured values by meteorological station

2.2 研究區(qū)蒸散量的時空變化特征

利用SEBAL模型估算研究區(qū)2001～2017年平均單位面積蒸散量，并統(tǒng)計出不同年份平均蒸散量的方差和變異系數(shù)于表2中。可以看出，17年間研究區(qū)平均單位面積蒸散量介于546.32～627.37mm，并且單位面積的蒸散量變異性較小，位于4.82%～6.67%之間。對于不同行政區(qū)域而言，瑪曲縣(5.0%)和紅原縣(4.78%)的單位面積蒸散量變異性最大、阿壩縣居中(3.95%)、碌曲縣(2.96%)和若爾蓋縣(2.68%)最小。擬合2001～2017年蒸散量的多年變化趨勢后發(fā)現(xiàn)，研究區(qū)單位面積蒸散量總體上呈微弱波動上升趨勢(圖2)。

為了分析2001～2017年研究區(qū)蒸散量的空間分布格局，將單位面積蒸散量按照400～500mm、500～600mm、600～700mm和700～800mm劃分為4個等級(圖3)，從圖3中可知，同一年份不同區(qū)域的蒸散量具有顯著的差異性，且不同年份蒸散量的空間分布也具有差異性?？傮w上看，研究區(qū)單位面積蒸散量從西北向東南呈現(xiàn)出逐漸增加趨勢。對于不同行政區(qū)域而言，17年平均單位蒸散量從大到小依次為紅原縣(622.11mm)、瑪曲縣(587.49mm)阿壩縣(580.29mm)、碌曲縣(572.70mm)、和若爾蓋縣(565.50 mm)。

表2 2001～2017年研究區(qū)單位面積蒸散量年均值統(tǒng)計Tab.2 Statistics of average annual evapotranspiration per unit area of Ruoergai Plateau

圖2 若爾蓋高原單位面積年均蒸散量時間變化趨勢Fig.2 Change trend of average annual evapotranspiration per unit area of Ruoergai Plateau

圖3 2001～2017年研究區(qū)蒸散量空間分布格局 Fig.3 Distribution pattern of evapotranspiration of Ruoergai Plateau from 2001 to 2017

2.3 研究區(qū)蒸散量時間變化趨勢的空間分布特征

基于2001～2017年蒸散數(shù)據(jù)，逐像元計算蒸散量變化的斜率，并分為7個相應的蒸散變化等級：顯著減少(slope<-10)、中等減少(-1010)(圖4)。從圖4中可以看出，研究區(qū)大部分都屬于輕微增加等級，部分中等增加等級的斑塊點綴其間，靠近東南邊界有少量輕微減少等級的區(qū)域。輕微減少、輕微增加和中等增加區(qū)域所占的面積比例分別為2.11%、97.20%、0.68%，其余等級區(qū)域所占面積幾乎為0。從行政區(qū)域上講，研究區(qū)內(nèi)五縣均屬于輕微增加趨勢。

圖4 研究區(qū)2001～2017蒸散量年際變化趨勢Fig.4 Spatial distribution map of inter-annual evapotranspiration trend in Ruoergai Plateau from 2001 to 2017

2.4 不同植被類型下的蒸散量變化特征

本研究選擇了研究區(qū)內(nèi)三種主要的植被類型(濕地、草地和林地)作為分析對象，基于土地利用數(shù)據(jù)和蒸散數(shù)據(jù)，統(tǒng)計2001～2017年單位面積的平均蒸散量(表3)。從表3中可以看出，不同植被類型單位面積蒸散量具有差異，且同種植被在不同年份的單位面積蒸散量也變化較大，具體表現(xiàn)為濕地>草地>林地，17年單位面積蒸散量分別為598.17mm、587.46mm和549.83mm。分析三種植被類型單位面積蒸散量的年際變化趨勢后發(fā)現(xiàn)(圖5)，濕地、草地和林地單位面積蒸散量均呈現(xiàn)微弱增加趨勢，且三者表現(xiàn)出同步變化規(guī)律。

表3 2001～2017年不同植被類型下單位面積蒸散量Tab.3 Statistics of average annual evapotranspiration per unit area for different vegetation types (mm)

圖5 2001～2017年三種植被類型下 的單位面積蒸散量變化趨勢Fig.5 Variation tendency of average annual evapotranspiration per unit area for different vegetation types from 2001 to 2017

3 討 論

本文利用SEBAL模型計算得到2001～2017年若爾蓋高原平均單位面積實際蒸散量介于546.32～627.37mm之間，小于劉佳等[18]和王建兵等[19]利用傳統(tǒng)FAO56 P-M方法計算的若爾蓋高原蒸散量，分別為1 200～1 600 mm和700～820 mm。這是由于FAO56 P-M法計算的為潛在蒸散量(ET0)，不能代表實際蒸散量(ET)。李志威等[17]將各類下墊面面積作為權(quán)重因子，利用FAO56 P-M 法加權(quán)計算若爾蓋高原的實際蒸散量，結(jié)果表明1961～2011年若爾蓋高原平均單位面積實際蒸散量介于452.6～535.2 mm之間。該結(jié)果略小于本文計算結(jié)果，但量級一致，造成這種差異的原因可能為研究方法及時空尺度不同。本文2001～2017年的多年平均蒸散量為587.06mm，與葉紅等[20]利用改進的MODIS ET算法得到的2000～2014年多年平均蒸散量538.61mm類似，這可能是由于所用數(shù)據(jù)一致和研究年份相近的緣故。

研究區(qū)2001～2017年蒸散量整體上呈波動微弱增加趨勢，其中2001～2008年蒸散量呈降低趨勢，2009～2017年蒸散量呈波動增加趨勢(圖2)，這可能與研究區(qū)實施生態(tài)保護密切相關[21-23]。同時與葉紅等[20]得到的結(jié)果類似：研究區(qū)蒸散量在2002、2008、2010和2013年處于低谷，且時間段整體和分段趨勢一致。針對不同的土地利用類型，蒸散量從大到小依次為濕地>草地>林地，這也與文獻[24-25]的研究結(jié)論相似。在模型參數(shù)設置方面，相比于傳統(tǒng)FAO56 P-M方法，SEBAL模型能更好的對不同植被類型下的實際蒸散量進行計算。此外，本文結(jié)果從側(cè)面印證了若爾蓋高原濕地出現(xiàn)了明顯的暖干化趨勢。由于若爾蓋高原是黃河上游重要的生態(tài)屏障，同時也是氣候變化的敏感區(qū)和脆弱帶，本文側(cè)重于若爾蓋高原實際蒸散量時空變化特征分析，后續(xù)工作將對影響若爾蓋高原實際蒸散量的氣象因子進行研究。

4 結(jié) 論

本文基于遙感和氣象臺站數(shù)據(jù)，利用SEBAL模型反演了若爾蓋高原的蒸散量，并采用統(tǒng)計分析方法定量評價了區(qū)域蒸散的時空變化特征，得出主要結(jié)論如下。

4.1 SEBAL模型的反演結(jié)果與地面實測數(shù)據(jù)的相關系數(shù)介于0.53～0.75之間，計算誤差小于35%，說明SEBAL模型在若爾蓋高原具有較好的適用性。

4.2 2001～2017年研究區(qū)單位面積年均蒸散量介于546.32～627.37mm，變異性較小(<7%)。從時間上看，單位面積年均蒸散量呈波動上升趨勢，其中輕微減少、輕微增加和中等增加等級區(qū)域所占的面積比例分別為2.11%、97.20%、0.68%，其余等級區(qū)域所占面積幾乎為0%。從行政區(qū)域上講，研究區(qū)內(nèi)五縣均屬于輕微增加趨勢。

4.3 從空間上看，研究區(qū)單位面積蒸散量呈現(xiàn)出從西北向東南逐漸增加趨勢。不同行政區(qū)域而言，17年平均單位面積蒸散量從大到小依次為紅原縣(622.11mm)>瑪曲縣(587.49mm)>阿壩縣(580.29mm)>碌曲縣(572.70mm)>若爾蓋縣(565.50 mm)。

4.4 研究區(qū)不同植被類型(草地、林地和濕地)單位面積蒸散量均呈現(xiàn)微弱增加趨勢，且三者表現(xiàn)出同步變化規(guī)律，17年平均單位面積蒸散量表現(xiàn)為濕地(598.17mm)>草地(587.46mm)>林地(549.83mm)。