張振宇,李小玉,孫 浩

1 浙江農(nóng)林大學(xué)林業(yè)與生物技術(shù)學(xué)院, 杭州 311300 2 中國科學(xué)院新疆生態(tài)與地理研究所荒漠與綠洲生態(tài)國家重點實驗室, 烏魯木齊 830011 3 中國科學(xué)院大學(xué),北京 100049

地表蒸散(Evapotranspiration,簡稱ET)是地球表面與大氣相互作用過程中的關(guān)鍵一環(huán),是地球水量平衡的重要組成部分,地表蒸散的過程伴隨著能量的吸收與釋放,也是維持地球能量平衡的重要環(huán)節(jié)。近年來,針對全球氣候變化的研究不斷深入,對地表蒸散的研究也越來越多,在這個過程中產(chǎn)生了較多的ET估算模型,如SEBAL模型[1]、SEBS模型[2]、ETwatch模型[3]等,其中SEBAL模型是應(yīng)用最為廣泛的ET估算模型之一。SEBAL模型是由Bastiaanssen提出的一種地表蒸散估算模型,研究者利用SEBAL模型結(jié)合不同種類遙感數(shù)據(jù)進行了不同區(qū)域的蒸散反演,并且取得了較好的反演結(jié)果。Santos[4]利用Landsat5 TM數(shù)據(jù)結(jié)合SEBAL模型進行了巴西布吉達河流域蒸散反演,Rawat等[5]基于Landsat7 TM數(shù)據(jù)和SEBAL模型對印度哈里亞邦區(qū)域的小麥農(nóng)田進行了蒸散反演,Alemayehu等[6]利用MODIS數(shù)據(jù)結(jié)合SEBAL模型進行了東非瑪拉盆地的地表蒸散估算。曾麗紅、杜嘉等[7-9]利用Landsat TM和MODIS數(shù)據(jù)結(jié)合SEBAL模型對中國東北地區(qū)的扎龍濕地、三江平原以及松嫩平原地區(qū)進行了蒸散量的反演與分析；張殿君等[10]利用SEBAL模型對中國黃土高原典型流域區(qū)進行了蒸散反演,并討論了土地利用對蒸散的影響；張發(fā)耀等[11]利用MODIS數(shù)據(jù)結(jié)合SEBAL模型對浙江省蒸散量進行了估算；王秋云[12]利用SEBAL模型進行了北京平原造林區(qū)的蒸散研究；周玲等[13]基于Landsat TM/ETM+數(shù)據(jù),應(yīng)用SEBAL模型對中國西南地區(qū)漓江流域進行了蒸散量變化分析；楊肖麗等[14]利用MODIS數(shù)據(jù)和SEBAL模型對中國半干旱區(qū)的老哈河流域進行了蒸散研究；李寶富、陳亞寧等[15]利用TM數(shù)據(jù),結(jié)合SEBAL模型對中國干旱區(qū)的塔里木河干流區(qū)進行了蒸散估算；張楠楠等[16]利用HJ- 1B數(shù)據(jù)結(jié)合SEBAL模型對淮河上游段進行了蒸散估算。以上研究都是利用SEBAL模型實現(xiàn)了地表蒸散的反演,而對SEBAL模型的參數(shù)計算方法影響討論較少,只有少數(shù)研究討論了SEBAL模型反演結(jié)果受氣象條件[17]和地表類型的影響。地表反照率,是指地表物體各個方向上發(fā)射的太陽輻射通量與到達該物體表面上的總輻射通量之比,是各個方向上反射率的積分。地表反照率是影響地球能量平衡的關(guān)鍵變量,相關(guān)研究[18-19]表明,地表反照率的增加促使地表凈輻射減少,造成區(qū)域降水和云量的降低,同時,云量的減少使得地表凈輻射通量進一步增大,促進降水和云的形成,這種反饋作用維持著區(qū)域氣候與熱量的穩(wěn)定與平衡。在遙感反演研究中,針對同一物理量的多種求解方法可能導(dǎo)致衍生結(jié)果的差異。地表反照率作為影響地表能量平衡的關(guān)鍵變量,是SEBAL模型的重要參數(shù),分析不同地表反照率計算方法造成的SEBAL模型反演結(jié)果差異,對提高SEBAL模型反演精度具有重要意義。本文擬在控制SEBAL模型其他參數(shù)一致的基礎(chǔ)上,基于兩種地表反照率計算方法,對研究區(qū)進行蒸散反演,通過反演值與實測值的擬合程度、偏離程度差異,實現(xiàn)不同地表反照率計算方法對SEBAL模型的影響分析。

1 研究區(qū)概況

三工河流域(圖1),位于新疆阜康市境內(nèi),區(qū)域主要由綠洲與部分荒漠構(gòu)成,地勢南高北低。平原區(qū)北部為古爾班通古特沙漠,南部為天山山脈,境內(nèi)自西向東分布有水磨河、三工河、四工河3條河流,區(qū)域內(nèi)干旱少雨,年降水量150—240 mm,水資源主要依賴高山冰川和積雪融水。研究區(qū)氣候類型為典型的溫帶大陸性氣候,夏季日照時間長,熱量高,冬季較早,寒冷漫長,流域內(nèi)主要農(nóng)作物為小麥、玉米等經(jīng)濟作物。

圖1 三工河流域位置示意圖Fig.1 Location of Sangong river basin

2 研究數(shù)據(jù)及預(yù)處理

2.1 遙感和其他空間數(shù)據(jù)獲取與預(yù)處理

①Landsat 8 OLI/TIRS影像數(shù)據(jù)來源于美國地質(zhì)調(diào)查局網(wǎng)站(https://glovis.usgs.gov/),共六期,日期分別為2016年5月25日、2016年6月26日、2016年7月12日、2016年7月28日、2016年8月29日和2016年9月30日。影像數(shù)據(jù)整體云量小于0.6%,其中6月26日和7月12日影像云量較高,分別為1.1%和7.45%。對6期Landsat 8 OLI/TIRS影像數(shù)據(jù)進行輻射定標(biāo)、大氣校正獲取糾正后數(shù)據(jù)。

②ASTER高程數(shù)據(jù)來源于馬里蘭大學(xué)全球土地覆蓋數(shù)據(jù)庫(http://www.landcover.org/data/),空間分辨率為30 m,對高程數(shù)據(jù)進行拼接與裁剪得到所需研究區(qū)高程數(shù)據(jù)。

③三工河流域矢量邊界由中國科學(xué)院新疆遙感與地理信息系統(tǒng)重點實驗室提供。

2.2 波文比數(shù)據(jù)

波文比數(shù)據(jù)主要來自于中國科學(xué)院新疆遙感與地理信息系統(tǒng)重點實驗室設(shè)立的2臺波文比觀測站,經(jīng)緯度坐標(biāo)分別為(44.33°N,87.89°E)和(44.36°N,87.87°E),波文比觀測站放置于均一的農(nóng)田內(nèi)部,基于波文比-能量平衡法(BREB)[20]獲得瞬時蒸散,進一步獲得每小時蒸散和日蒸散。大量研究證明[21-22]利用波文比觀測儀器測量蒸散精度較高。波文比數(shù)據(jù)監(jiān)測間隔為1h。

3 模型及參數(shù)

SEBAL模型[23]的理論基礎(chǔ)是地表能量平衡方程,即地表的凈輻射通量由從下墊面到大氣的感熱通量、從下墊面到大氣的潛熱通量和土壤熱通量組成：

Rn=G+H+λET

(1)

式中,Rn代表凈輻射通量(W/m2),G代表土壤熱通量(W/m2),H代表感熱通量(W/m2),λET代表潛熱通量(W/m2),λ代表水的汽化潛熱,通常為2.49×106W m-2mm-1；ET表示蒸散量,kg m-1s-1。

3.1 凈輻射通量(Rn)計算

凈輻射通量計算公式為：

Rn=(1-a)Rsd+Rld-Rlu-(1-k)Rld

(2)

式中,a為地表反照率,Rsd為太陽下行短波輻射(W/m2),Rld為太陽下行長波輻射(W/m2),Rlu為地表反射的長波輻射(W/m2),k為地表比輻射率(發(fā)射率)。

3.2 土壤熱通量(G)的計算

土壤熱通量是單位時間、單位面積上的土壤熱交換量,一般存儲在土壤和植被中,這部分能量在整個熱量平衡過程中占比重較小,而且計算較為困難,這里使用曾麗紅[24]推導(dǎo)的公式。

(3)

式中,Ts為地表溫度,NDVI為歸一化植被指數(shù),c表示衛(wèi)星過境時間對G的影響,過境時間在地方時12：00點以前取0.9；在12：00點到14：00點之間取1.0；在14：00點到16：00點之間取1.1。

3.3 感熱通量(H)的計算

感熱通量表征的是下墊面與大氣間的湍流形式的熱量交換,與大氣穩(wěn)定度、風(fēng)速、表面粗糙度有關(guān)。通用計算公式如下：

(4)

式中,Pair為空氣密度(kg/m3)；Cp是空氣定壓比熱(1004 J kg-1K-1)；dT是距離地面高Z1和Z2處的溫差(一般Z1為0.1 m,Z2為2 m),rah為空氣動力學(xué)阻抗。

3.4 小時蒸散量與日蒸散量計算

小時蒸散量的計算公式[25]為：

(5)

λ=2.501-0.002361×(Ts-273.15)

(6)

式中,λ為水的汽化潛熱,ET為蒸散量,ETinst為小時蒸散量。

通過小時蒸散量結(jié)合日照時長計算日蒸散量ETday：

(7)

式中,N為區(qū)域內(nèi)平均日照時長,t為日出到衛(wèi)星過境時間間隔(h)。

4 地表反照率計算

4.1 基于各波段反射率的大氣表觀反射率計算

求取大氣表觀反射率,這一過程在輻射定標(biāo)中即可完成：

(8)

式中,ρ代表表觀反射率,DN為影像灰度值,μ、φ分別為landsat8 OLI/TIRS影像的反射率增量和修正值,θ為太陽高度角。

4.2 求取地表反照率

①Liang地表反照率計算方法(下稱Liang方法)[26]:

(9)

式中,ρ2、ρ4、ρ5、ρ6、ρ7為各波段大氣表觀反射率,ap為大氣程輻射值影響,一般取0.03,tsw為大氣單向透過率。

大氣單向透過率計算公式如下[27]：

tsw=0.75+2×10-5×Z

(10)

式中,Z為海拔高度(m),從DEM數(shù)據(jù)中獲取。

②Smith地表反照率計算方法(下稱Smith方法)[28]:

(11)

式中,ρi為波段反射率(對于Landsat 8 OLI/TIRS數(shù)據(jù)使用2、4、5、6、7波段),Δλi為各波段范圍,Si為1—7波段的平均太陽輻照度(表1),ap為大氣程輻射值影響,一般取0.03,tsw為大氣單向透過率。

表1 Landsat 8 OLI數(shù)據(jù)波段平均太陽輻照度[28]

Si,平均太陽輻照度,Average solar irradiance

5 結(jié)果與分析

5.1 日蒸散量比較

基于SEBAL模型得到兩種地表反照率計算方法下反演的地表日均蒸散量(表2)及其空間分布圖(圖2和圖3),對比6個研究日區(qū)域日蒸散量的分布,結(jié)果顯示兩種地表反照率方法下最終得到的日蒸散量空間和時間分布規(guī)律保持一致。時間上,日蒸散量在6月—7月末明顯高于其他月份；空間上,流域綠洲的日蒸散量明顯高于其他類型區(qū)域(水體除外)；另一方面,兩種地表反照率方法下基于SEBAL模型得到的日蒸散量分布曲線(圖4)則有較大不同。5月25日,兩種地表反照率方法下的日蒸散量分布曲線在形態(tài)上保持一致,但是Liang方法下的日蒸散量分布曲線在6 mm/d處開始接近橫坐標(biāo)軸,而Smith方法下日蒸散量分布曲線在日蒸散量7 mm/d處才開始接近橫坐標(biāo)軸。6月26日、7月12日、7月28日三天同一地表反照率計算方法下的日蒸散量分布曲線形態(tài)上較為一致,但是不同地表反照率計算方法下得到的曲線則差異較大：Liang方法下得到的3期日蒸散量分布曲線波峰出現(xiàn)在日蒸散量4 mm/d附近,在日蒸散量6 mm/d處開始接近橫坐標(biāo)軸,而Smith方法下的3期日蒸散量分布曲線均在日蒸散量6—6.5 mm/d處達到峰值,在日蒸散量7.5 mm/d處才開始接近橫坐標(biāo)軸。8月29日和9月30日兩期日蒸散量分布曲線則呈現(xiàn)出與前4期較為不同的特征：曲線波峰開始向低日蒸散量區(qū)間移動,Liang方法下得到的日蒸散量分布曲線波峰分布在日蒸散量1—1.5 mm/d附近,而Smith方法下得到的日蒸散量分布曲線在日蒸散量1.5—2.5 mm/d處達到峰值。兩種方法下最終得到的多條日蒸散量分布曲線在日蒸散量區(qū)間末端均出現(xiàn)了一個較為明顯的波峰,是由區(qū)域內(nèi)的水庫水體所引起?？偟膩碚f,Smith方法下基于SEBAL模型反演的區(qū)域日蒸散量分布曲線峰值出現(xiàn)更晚。

此外,兩種地表反照率計算方法下最終求得的日均蒸散量也有較大差異,Liang方法下最終得到的區(qū)域日均蒸散量在所選取的6個研究日內(nèi)均小于Smith方法下最終得到的結(jié)果,其中,春季(5月25日)、秋季(8月29日和9月30日)的日均蒸散量差距較小,平均差異0.2 mm/d,而夏季(6月26日、7月12日、7月28日)的日均蒸散量差距較大,平均差異0.64 mm/d。

表2 兩種地表反照率計算方法下得到的區(qū)域日均蒸散量/(mm/d)

圖2 Liang地表反照率計算方法下日蒸散量分布Fig.2 Daily evapotranspiration distribution using Liang surface albedo formula

圖3 Smith地表反照率計算方法下日蒸散量分布Fig.3 Daily evapotranspiration distribution using Smith surface albedo formula

圖4 Liang反照率計算方法和Smith反照率計算方法下的日蒸散量分布曲線Fig.4 Distribution curve of daily evapotranspiration using Liang and Smith surface albedo method

5.2 兩種地表反照率計算方法下的日蒸散反演精度比較

將SEBAL模型反演值與實測值進行擬合,并使用決定系數(shù)R2、均方根誤差RMSE、平均偏差BIAS和平均絕對偏差MAE[29]對兩種地表反照率方法下基于SEBAL模型反演的日蒸散量進行比較與評價。評價指標(biāo)計算公式如下:

(10)

(11)

(12)

式中,N為觀測個數(shù),Pi為反演值,Qi為實測值。得到的結(jié)果如表3所示。

表3 兩種地表反照率計算方法下SEBAL模型反演值精度比較

R2,決定系數(shù),Coefficient of determination；RMSE,均方根誤差,Root mean square error；BIAS,平均偏差,Bias；MAE,平均絕對偏差,Mean absolute error

兩種地表反照率計算方法下基于SEBAL模型反演的日蒸散量與實測值散點圖及擬合方程(圖5)表明,Smith方法下最終得到的日蒸散量與實測值的線性擬合程度更高,決定系數(shù)達到了0.85,相比之下,Liang方法下得到的模型反演日蒸散量與實測值的線性擬合決定系數(shù)僅為0.76；F檢驗顯示,Smith方法下最終得到反演值與實測值線性擬合P值<0.001,而Liang反照率計算方法下得到的擬合直線P值為0.0016,大于0.001,可見相同時間空間條件下,Smith方法下基于SEBAL模型得到的日蒸散量與實測值更為接近,可信程度更高。另外通過RMSE、BIAS、MAE 3種指標(biāo)對反演值與實測值進行偏離程度分析,Smith方法下最終得到的日蒸散量與實測值的偏離程度小于使用Liang方法下得到的結(jié)果,多項結(jié)果表明：在應(yīng)用SEBAL模型對干旱區(qū)地表進行蒸散反演過程中,使用Smith方法結(jié)果作為SEBAL模型地表反照率參數(shù)輸入更優(yōu)。

圖5 Liang反照率計算方法和Smith反照率計算方法下的日蒸散量擬合Fig.5 Daily evapotranspiration fitting using Liang and Smith surface albedo method

6 討論

6.1 兩種方法計算下的地表反照率比較

2016年6個研究日的流域地表反照率分布曲線(圖6)表明,兩種方法得到的地表反照率分布曲線類似,都呈較明顯的正態(tài)分布,然而不同方法下的曲線波峰處反照率存在明顯差異,基于Liang方法得到的地表反照率分布曲線波峰處反照率在6個時期均高于Smith方法得到的曲線波峰反照率,波峰處地表反照率相差0.02—0.06,盡管這種差異較小,但是這種差異使得整個區(qū)域的反照率分布出現(xiàn)較大不同,Smith方法下得到的地表反照率更低,兩種方法計算出的地表反照率均值(表4)可進一步印證這一點,地表反照率的差異使得由此計算出的流域凈輻射通量、顯熱通量產(chǎn)生了更大的差異,進一步的造成了兩種地表反照率計算方法下基于SEBAL模型反演結(jié)果在分布上的差異。

圖6 兩種方法計算下反照率像元分布曲線Fig.6 Albedo pixel distribution curves under two methods

地表反照率計算方法Surface albedo calculation method5月25日May 25th6月26日June 26th7月12日July 12th7月28日July 28th8月29日August 29th9月30日September 30thLiang0.36380.34690.32850.31830.30130.2916Smith0.29440.26840.25520.25580.25200.2520

6.2 兩種方法計算下地表反照率與MCD43A3產(chǎn)品的比較

圖7 地表反照率檢測點位置Fig.7 Location of surface albedo detection points

為進一步比較兩種方法計算出的地表反照率在研究區(qū)的可靠性,在缺乏地表反照率測量儀器的情況下,參考He等[30]的方法,使用MODIS地表反照率產(chǎn)品MCD43A3與兩種方法計算出的地表反照率進行比較,由于研究區(qū)所選Landsat 8數(shù)據(jù)云量均較少(晴空條件),因此提取MCD43A3數(shù)據(jù)中的短波白空反照率(WSA)作為參考標(biāo)準(zhǔn)。在研究區(qū)內(nèi)部隨機布置了36個監(jiān)測點(圖7),通過6個研究日共得到216個數(shù)據(jù)。數(shù)據(jù)比對結(jié)果(圖8)表明：①兩種方法下得到的地表反照率均偏高于MCD43A3 WSA,Liang方法計算出的地表反照率與MCD43A3 WSA擬合程度(R2=0.35)略高于Smith方法計算結(jié)果(R2=0.2),說明Liang方法計算出的地表反照率一致性更好。②RMSE比較上,Smith方法計算出的地表反照率與MCD43A3 WSA誤差相對較小(Smith方法RMSE=0.075；Liang方法RMSE=0.13)。

圖8 兩種方法計算下地表反照率與MCD43A3短波白空反照率對比Fig.8 Comparison between surface albedo calculated by two methods and MCD43A3 shortwave white sky albedo

同時由于Landsat 8數(shù)據(jù)部分參數(shù)官方尚未明確公布,計算Landsat 8數(shù)據(jù)地表反照率時仍采用針對上一代傳感器的方法,這也可能會對結(jié)果造成一定影響；基于此,Smith提出了針對Landsat 8數(shù)據(jù)地表反照率的計算方法,并在阿拉斯加、佛羅里達等地區(qū)進行了實驗,近期部分學(xué)者[30]提出利用大氣頂反照率和站點測量數(shù)據(jù)直接估計Landsat 8地表反照率,實驗效果較好,但由于中國西北地區(qū)缺少地表反照率國際測量站點,因此,各種方法在中國西北地區(qū)的適用性仍有待進一步檢驗。

另一方面,由于USGS(美國地質(zhì)調(diào)查局)至今未公布Landsat 8數(shù)據(jù)各波段平均太陽輻照度,因此在使用Smith方法過程中使用了其提供的各波段平均太陽輻照度。在理想狀態(tài)下,平均太陽輻照度可通過大氣光譜曲線和各波段響應(yīng)函數(shù)獲得,分別使用WRC太陽光譜曲線和Thuillier 太陽光譜曲線結(jié)合Landsat 8數(shù)據(jù)1—7波段響應(yīng)函數(shù)進行平均太陽輻照度估計[31-32],結(jié)果(表5)顯示,兩種譜線對應(yīng)的波段平均太陽輻照度存在差異,其中1和4波段差異較大,其他波段結(jié)果相近,同時Smith方法的1、2、3波段平均太陽輻照度均小于兩種譜線對應(yīng)估計值,而在第6波段(短波紅外)上高于譜線估計值,這使得Smith方法在計算地表反照率時更側(cè)重于短波紅外波段,而由此得出的地表反照率與相關(guān)產(chǎn)品的偏差較小,這說明Landsat 8地表反照率計算中短波紅外波段的權(quán)重增加可能改善了研究區(qū)地表反照率反演質(zhì)量,部分研究[33]也表達了相似的結(jié)果。

表5 利用Thuillier和WRC太陽光譜曲線估算的Landsat 8 OLI數(shù)據(jù)各波段平均太陽輻照度

ThuillierSi,利用Thuillier太陽光譜曲線估算的波段平均太陽輻照度,Band average solar irradiance based on Thuillier Solar spectral curve；WRCSi,利用WRC太陽光譜曲線估算的波段平均太陽輻照度,Band average solar irradiance based on WRC Solar spectral curve

7 結(jié)論

研究基于兩種地表反照率計算方法結(jié)合SEBAL模型,對西北干旱區(qū)典型流域日蒸散進行反演并進行比較,結(jié)果表明：(1)Simth方法下最終得到的日均蒸散量均高于Liang方法下最終得到的日均蒸散量,這種差異在夏季最明顯。同時,Smith方法下基于SEBAL模型得到的反演結(jié)果與實測值的擬合程度和偏離程度略優(yōu)于Liang方法下得到的結(jié)果。(2)兩種方法計算的地表反照率差異是導(dǎo)致SEBAL模型反演日蒸散量差異的主要原因,相同日期內(nèi)Liang方法得到的地表反照率均高于Smtih方法結(jié)果(高約0.04—0.08)。(3)與MCD43A3產(chǎn)品的比較結(jié)果顯示,Liang方法得到的地表反照率與MCD43A3產(chǎn)品一致性較好,但是其偏離程度也較高。綜合考慮之下,Smith方法在研究區(qū)的適用性略好。