地表反照率不同計算方法對干旱區流域蒸散反演結果的影響
——以新疆三工河流域為例

張振宇,李小玉,孫 浩

1 浙江農林大學林業與生物技術學院, 杭州 311300 2 中國科學院新疆生態與地理研究所荒漠與綠洲生態國家重點實驗室, 烏魯木齊 830011 3 中國科學院大學,北京 100049

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

1 研究區概況

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

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

2 研究數據及預處理

2.1 遙感和其他空間數據獲取與預處理

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

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

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

2.2 波文比數據

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

3 模型及參數

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

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為地表比輻射率(發射率)。

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

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

(3)

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

3.3 感熱通量(H)的計算

感熱通量表征的是下墊面與大氣間的湍流形式的熱量交換,與大氣穩定度、風速、表面粗糙度有關。通用計算公式如下：

(4)

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

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

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

(5)

λ=2.501-0.002361×(Ts-273.15)

(6)

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

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

(7)

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

4 地表反照率計算

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

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

(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數據中獲取。

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

(11)

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

表1 Landsat 8 OLI數據波段平均太陽輻照度[28]

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

5 結果與分析

5.1 日蒸散量比較

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

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

表2 兩種地表反照率計算方法下得到的區域日均蒸散量/(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模型反演值與實測值進行擬合,并使用決定系數R2、均方根誤差RMSE、平均偏差BIAS和平均絕對偏差MAE[29]對兩種地表反照率方法下基于SEBAL模型反演的日蒸散量進行比較與評價。評價指標計算公式如下:

(10)

(11)

(12)

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

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

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

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

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

6 討論

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

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

圖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產品的比較

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

為進一步比較兩種方法計算出的地表反照率在研究區的可靠性,在缺乏地表反照率測量儀器的情況下,參考He等[30]的方法,使用MODIS地表反照率產品MCD43A3與兩種方法計算出的地表反照率進行比較,由于研究區所選Landsat 8數據云量均較少(晴空條件),因此提取MCD43A3數據中的短波白空反照率(WSA)作為參考標準。在研究區內部隨機布置了36個監測點(圖7),通過6個研究日共得到216個數據。數據比對結果(圖8)表明：①兩種方法下得到的地表反照率均偏高于MCD43A3 WSA,Liang方法計算出的地表反照率與MCD43A3 WSA擬合程度(R2=0.35)略高于Smith方法計算結果(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數據部分參數官方尚未明確公布,計算Landsat 8數據地表反照率時仍采用針對上一代傳感器的方法,這也可能會對結果造成一定影響；基于此,Smith提出了針對Landsat 8數據地表反照率的計算方法,并在阿拉斯加、佛羅里達等地區進行了實驗,近期部分學者[30]提出利用大氣頂反照率和站點測量數據直接估計Landsat 8地表反照率,實驗效果較好,但由于中國西北地區缺少地表反照率國際測量站點,因此,各種方法在中國西北地區的適用性仍有待進一步檢驗。

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

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

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 結論

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