基于能量平衡算法的精河流域綠洲蒸散發時空變化模擬

于 輝 ，代鵬超，張金燕，毋兆鵬,2

(1.新疆師范大學地理科學與旅游學院，新疆 烏魯木齊 830054；2.新疆干旱區湖泊環境與資源實驗室，新疆 烏魯木齊 830054)

水資源利用是干旱區資源環境和生態系統的核心問題，其中蒸散是干旱內陸水循環中水分消耗的最終途徑，也是區域水量平衡和能量平衡的最活躍因子。精確估算地表蒸散發的時空變化，對于了解水文循環和能量平衡的整個過程，評價區域水循環和水平衡的功能，實現水資源、生態、社會和諧發展具有十分重要意義[1-2]。

1802年，Dalton提出關于蒸散發的計算公式，開啟了國內外研究蒸散發的歷史性時代。隨后發展了一系列蒸散發觀測和計算方法[3]。應用遙感技術是蒸散發估算領域新的一大亮點，為區域空間尺度的蒸散發擴展計算找到了解決方法，許多適用不同區域的和改進的蒸散發計算模型就是由此而生，如陸地表面能量平衡算法(Surface energy balance algorithm for land ,SEBAL)[4]、地表能量平衡指數(Surface energy balance index, SEBI)[5]、簡化的地表能量平衡指數(S-SEBAL)[6]、地表能量平衡系統(Surface energy balance system, SEBS)[7]、MODIS蒸發比模型[8]和雙層蒸散模型TTME[9]等。

研究區屬天山北坡地區，位于全國“兩橫三縱”城市化戰略格局中陸橋通道的西端，是全國主體功能區規劃[10]確定的國家層面重點開發區域。研究區是北疆重要的優質棉基地之一，但隨著近幾十年區域人口增加及城市化進程加劇，加之干旱氣候條件、水資源分布空間差異的綜合作用，在綠洲開發過程中出現了一系列生態危機。2015年3月，新疆被確定為“絲綢之路”經濟帶核心區，使研究區從對外開放的“末梢”變成了開放前沿,也更使得這里成為了新時期我國開展生態環境國際合作和生態文明建設的重點區域。目前針對該區域蒸散發的研究雖已有部分成果，但數量較少且難以反映在“山地-綠洲-荒漠”系統背景下蒸散發的空間異質性，對全區未來蒸散發變化趨勢的定量研究還相對比較薄弱。因此，本文利用基于能量平衡開發的Surfance energy balance alogrithm for land (SEBAL)模型，利用Landsat_5 TM遙感數據和地面氣象數據，對精河流域綠洲地表通量及日蒸散發進行了估算，同時分析了蒸散發的時空變化規律和未來的變化趨勢，以期為該區域水資源管理規劃及其可持續發展提供科學依據。

1 研究區概況

研究區精河流域綠洲，位于新疆維吾爾自治區西北部準噶爾盆地西南邊緣，天山支脈婆羅科努山北麓，介于82°32′50″E-83°17′49″E、44°36′30″N-45°08′23″N之間，包括整個艾比湖湖區及精河流域平原綠洲，總面積3 503.41 km2(圖1)。該地區地處亞歐大陸腹地，距海洋較遠，地勢南高北低，屬于北溫帶干旱荒漠型大陸性氣候，日照充足，冬冷夏熱，干燥少雨，蒸發量大，冬夏長，春秋時間較短。年平均太陽總輻射量冬季和夏季相差較大，最大可達到36.7 kcal·cm-2。降水量由南部山區向中、北部平原逐漸減少，年均降水102 mm左右，其中山區降水較為豐富，最大時可達約700 mm。研究區全年的蒸發主要集中在4～9月，多年平均蒸發量約為1 625 mm，且蒸發量大于降水量。

圖1 研究區位置圖Fig.1 Location of the study area

2 數據來源與研究方法

2.1 數據來源

研究所使用的1998、2007年和2011年遙感影像數據為Landsat-5 TM影像，拍攝日期分別為9月25日、9月18日、9月13日，三期數據含云量均低于5%，影像清晰度良好，經過幾何精校正、輻射定標、Flaash大氣校正后適合本次建模的基礎條件。

土地利用數據參考中國科學院土地資源遙感調查與監測技術規程[11]，結合研究區的實際情況，將研究區土地利用類型劃分為耕地、林地、草地、水域、建設用地、未利用土地6個類型，基于TM影像結合地形圖、專題圖等相關資料自行解譯完成，并通過調查驗證，誤差控制在1個像元內，最終在ARC GIS平臺下生成研究區土地利用數據庫。

實際蒸散量主要根據以下5個氣象站點：阿拉山口氣象站點(站點編號：51232)、托里氣象站點(站點編號：51241)、溫泉氣象站點(站點編號：51330)、精河氣象站點(站點編號：51334)以及烏蘇氣象站點(站點編號：51346)的日值數據，利用聯合國糧農組織推薦的Penman-Monteith公式，并借助表征作物在不同生長時期中需水量的生物學特性的作物系數Kc計算獲得。

2.2 基于能量平衡算法的SEBAL模型

Bastiaanssen博士在1998年提出SEBAL模型，已經在多個國家和地區得到了廣泛地應用。該模型依據陸地表面能量平衡原理，即地球表面所獲得的凈輻射量等于土壤熱通量、感熱通量和潛熱通量之和。SEBAL模型的計算公式為：

Rn=H+LE+G

(1)

其中,Rn為地表凈輻射通量(W·m-2)，H為感熱通量(W·m-2),LE為潛熱通量(W·m-2),G為土壤熱通量(W·m-2)。

2.2.1 地面凈輻射(Rn)的計算 地面凈輻射是地表的主要能量來源，是地表通過太陽輻射的短波輻射、長波輻射過程得到的凈能量(圖2a)，計算公式如下：

Rn=(1-α)Rs+(Lin-Lout)-(1-ε)Lin

(2)

其中：Rn為地表凈輻射(W·m-2)；Rs為太陽總輻射(W·m-2)；α為地面反射率；ε為地面比輻射率；Lin為大氣長波輻射；Lout為地面長波輻射。

圖2 研究區模型關鍵參數運算結果Fig.2 Diagram of operational results of key parameter of study area model

2.2.2 土壤熱通量(G)的計算 土壤熱通量指存儲在土壤或植被中的能量，是一個相對較小的量，直接計算較困難，一般通過G與地表溫度(Ts)、Rn、α、植被指數(NDVI)的統計關系求得，同時根據衛星過境時間進行適當的校正(圖2b)，計算公式如下：

(1-0.978NDVI4)×Rn

(3)

其中：c11為衛星過境時間，在地方時12∶00之前取0.9，12∶00-14∶00取1.0，14∶00-16∶00取1.1。

2.2.3 感熱通量(H)的計算 感熱通量是大氣穩定度、風速和表面粗糙度的函數(圖2c)，計算公式如下：

(4)

式中,ρair是空氣密度(kg·m-2)；Cp為空氣定壓比熱(≈1004J·kg-1·K-1)；dT(K)為零平面以上高度z1和z2處的溫差(T1-T2)；rah為熱量傳輸的空氣動力學阻抗(s·m-1)。

上式中，dT、rah均為未知量且彼此相關，因此借用Monin-Obukhov理論，通過循環遞歸算法進行求解(圖3)。圖中U為風速，U*為摩擦風速，Zom為動力粗糙度；L是Monin-Obukhov長度，反映地面層湍流特性的關鍵參量；Ψh、Ψm是感熱通量穩定度修訂函數，可以根據Paulson[12]提出的大氣穩定度修正公式計算得到。1998年rah值在33.698～35.794之間震蕩，取第14次值；2007年rah值在32.567～32.872之間震蕩，取第12次值；2011年rah值在34.152～34.083之間震蕩，取第13次值。

圖3 Monin-Obukhov循環遞歸流程圖Fig.3 The flow chart of Monin-Obukhov circular recursion

2.2.4 潛熱通量(LE)的計算 通過模型分別反演出了凈輻射、土壤熱通量和感熱通量后，運用公式5即可反演出衛星過鏡時刻的潛熱通量。

LE=Rn-H-G

(5)

公式中通量單位均為W·m-2，其中Rn、G、H都是基于像元得到的柵格數據。

2.2.5 全天蒸散發估算 運用瞬時潛熱通量計算出蒸發比，將瞬時值延伸至全天的蒸散發值。SEBAL模型中，假設24 h之內蒸發比EF是相對穩定的，即

(6)

衛星過境時刻的瞬時蒸散發量計算公式為：

ETinst=3 600×(Rn-H-G)/I

(7)

其中，I為汽化潛熱(J·kg-1),3 600是秒到小時的轉換系數，所以單位是mm·h-1。則24 h蒸散發ET24為：

(8)

式中,Rn24、G24為日平均凈輻射通量和日平均土壤熱通量。

2.3 蒸散發周期預測

小波是一種特殊的長度有限、平均值為0的波形，傾向于不規則和不對稱，能分析出時間序列周期性變化的局部特性，同時可以更清楚地看出各周期隨時間的變化情況[14]。本文采用Morlet小波研究ET序列的特征尺度和周期性，其小波函數為：

(9)

在此基礎上，結合Mann-Kendall(M-K)突變檢驗進行分析，序列大于0表明參數呈上升趨勢，反之呈下降趨勢，當其超過0.05顯著性水平檢驗臨界值(1.96)時，則上升或下降趨勢顯著[15]。

3 結果與分析

3.1 SEBAL模型精度驗證

研究區內缺乏實際地表通量觀測數據，因此，參考朱明承等[16]、劉文娟[17]和姜紅[18]的方法，先基于研究區及周邊站點氣象數據，采用聯合國糧農組織(FAO)推薦的P-M法計算參考作物蒸散發，再結合作物系數得出主要土地利用類型的實際蒸散發，并將此結果與遙感估算值相對比。結果表明，SEBAL模型計算結果相對誤差均在18%之內，能夠滿足對反演結果的精度要求(表1)。

表1 SEBAL模型精度驗證

Table 1 Precision validation of the SEBAL model

土地利用類型LandusetapeKC系數KCcoefficientFAOPM蒸散發均值/mmFAOPMevapotranspirationaverage遙感反演蒸散發均值/mmEvapotranspirationaveragebyremotesensing絕對誤差Absoluteerror相對誤差/%Relativeerror耕地Cultivatedland1．155．926．150．233．89林地Woodland1．085．566．430．8715．65草地Grassland0．854．385．150．7717．57

3.2 日蒸散量時空變化分析

3.2.1 日蒸散量空間分異特征 就空間尺度而言，研究區位于中國西北部典型干旱區，海拔落差較小，土地利用類型結構比較單一。在氣候、降水、溫度和土壤含水量等因素的綜合影響下，該區域的日蒸散發模型反演結果有顯著的空間差異性。圖4反映了研究區不同年份夏季的日蒸散發分布情況，結果表明中部艾比湖湖面為極高值范圍，這是由于湖面上方水分相對于其他地類比較充足，湖水比熱容大，湖區溫度上升和下降的幅度小。西南區除托里鄉旦達蓋沙漠外為蒸散發次高值區域，因為這部分區域地形較為平坦、海拔低，研究區耕地主要分布于此，基礎灌溉設施完善，棉花此時也正處于花鈴期末、蒸騰作用強烈，蒸散發值高。研究區的西北部是著名的風區阿拉山口，導致植被遭到破壞或覆蓋度低，土地沙化現象較為嚴重，植物的蒸發、蒸騰量較少。艾比湖的北部、東部梭梭、紅柳等灌木由于數年前的連年樵采而受到不同程度的破壞，導致地表覆蓋較差出現了蒸散量低值區。研究區東南部開始進入天山山區，地形較高造成小范圍水分的重新分配，林地面積占比高，故也造成該區域日蒸散發量較高。

從時間上看，1998-2007年，日蒸散發艾比湖湖區的高值范圍變化不大，但西南區的次高值區域明顯減少，取而代之的是低值區域的增加。2007-2011年湖區的日蒸散高值范圍明顯減少，同時西南區的次高值區域則明顯擴大，低值區域明顯減少。這一方面是由于該區域的耕地面積呈增加趨勢，且增加趨勢在加強；另一方面則是由于艾比湖濕地自然保護區在2007年晉升為國家級自然保護區，周圍生態環境得到保護而導致這一結果。

圖4 研究區日蒸散發空間分布Fig.4 Spatial distribution diagram of daily evapotranspiration in the study area

3.2.2 日蒸散量時間分異特征 利用精河氣象站近40年數據和P-M公式，結合Mann-Kendall 檢驗對1970-2015年研究區年實際蒸散量序列進行分析(圖5)，結果表明，1970-2015年正序列曲線(UF)總體為持續下降趨勢，UF和反序列曲線(UB)于1973年在臨界線區間內有1個明顯的交點，表明研究區蒸散發在1973年發生突變，突變后精河流域平均年蒸散量較突變前減少218.318 mm，降低幅度為22.019%。進一步根據1970-2015年精河流域實際蒸散量數據，繪制Morlet小波變換系數的實部等值線圖(圖6)。結果表明，研究區蒸散量主要存在26～30 a的周期變化。

為探究研究區未來的蒸散量變化情況，本文基于實際蒸散量小波方差(圖7)，對28 a處的小波系數值[f(t)]與年份t建立回歸方程：f(t)=357sin(0.337 6t-115.7)，R2=0.98，P<0.05，并選擇2018-2030年為預測期對研究區蒸散量變化趨勢進行預測(圖8)。結果表明，在現有條件不發生大的改變前提下，預測其整體蒸散量呈波動變化。其中，自當前到2022年蒸散量不斷增加，并將于2022年發生突變，進入下降期。預測到2030年，蒸散量將再次進入上升周期。

圖5 研究區蒸散量Mann-Kendall突變檢驗Fig.5 Mann-Kendall tests of the evapotranspiration rate in the study area for sudden changes

圖6 研究區蒸散量Morlet小波系數實部等值線圖Fig.6 Contour map isogram of real part of Morlet wavelet coefficients of ET in the study area

圖7 研究區實際蒸散發量小波方差Fig.7 Wavelet variances of actual ET in the study area

圖8 研究區實際蒸散發量小波系數(28 a)Fig.8 Wavelet coefficients of actual ET in the study area (28 a)

3.3 不同土地利用類型日蒸散發變化

不同土地利用類型使地表土壤的濕度和地表溫度狀況發生改變，影響陸面能量平衡，因此對區域日蒸散發的影響非常顯著。利用ARC GIS空間分析技術對各土地利用類型平均日蒸散發量的研究表明(圖9)，三期的土地利用類型的日平均蒸散發量均為：水域﹥林地﹥耕地﹥草地﹥建設用地﹥未利用地。

研究區的水域主要是艾比湖和精河綠洲的河流，由于水體位于北部和東北部地表比較平坦的區域，接受的太陽能輻射充足，因此實際蒸散發大于等于該區域的潛在蒸散發，也成為整個研究區日蒸散發最大區域，達到7.3～9.32 mm之間。研究區的林地主要分布在耕地邊緣和河流附近等海拔相對較低但水分較為充足的地帶，由于對地下水源的涵養能力比人工灌溉的耕地強，因此其日蒸散發量相比其他土地利用類型較高，均值達到6.43 mm。研究區耕地是主要土地利用類型且區域比較集中，由于主要作物棉花在9月初即將步入采摘期，棉農會停止澆水、施肥等田間管理活動，因此實際蒸散發情況相對于生長旺盛期而言較低，均值為6.15 mm。研究區草地自1998年到2011年呈現縮減的趨勢，多數轉變為耕地(表2)，加之草地的水源主要來自降水，因此由于本區域降水的缺乏，限制了該土地利用類型的蒸散發活動過程，使得該地區的實際蒸散發量不是很大，均值為5.37 mm。研究區建設用地在精河綠洲的占地極少，呈零星分布，主要由交通建設用地和建設用地等組成，不透水面積大，不利于水分擴散，導致潛熱通量變小，使得整體蒸散發較低，日蒸散發值主要在4.49～5.34 mm之間，在研究區土地利用類型的蒸散發量圖上表現出較低的區域。未利用土地在研究區為裸露巖石、戈壁、荒漠和沙地，是整個研究區分布最多的地類，由于地表覆蓋度極低致使該類型的日蒸散發量成為了整個研究區的極低值區域。

圖9 不同土地利用類型日平均蒸散量Fig.9 The daily mean ET of different land types in the study area

表2 1998-2011年研究區土地利用類型轉移矩陣/km2Table 2 Land use transit matrix of study area from 1998 to 2011

4 結論與討論

本文在蒸散發遙感反演的基礎上，綜合分析了研究區日蒸散發特征，得出以下結論：

1)研究區實際蒸散量近20 a間總體呈下降態勢，雖然1995-2004年出現小幅增長趨勢，但研究區蒸散量存在26～30 a的周期變化規律，故可預見2022年將再次轉入下降，而2030年蒸散量則再次進入上升周期。

2)研究區中部艾比湖湖面為極高值范圍，雖然東南部開始進入天山山區日蒸散發量較高，但西南區除托里鄉旦達蓋沙漠外的綠洲耕作區，實為較高值蒸散發的主體分布區域，西北部、北部、東部則為蒸散發低值區。在不同土地利用類型中，未利用地日蒸散發最小，其次是建設用地，除水體外林地和耕地蒸散發最高。

3)在自然界中，潛熱(LE)和感熱通量(H)都是伴隨著地面凈輻射通量(Rn)的變化而變化，同時LE和H的相互比例決定了能量的轉換耗散方式。研究區三期潛熱通量均值為450.6 W·m-2，感熱通量均值為392.6 W·m-2，前者大于后者說明研究區內植物正處于生長旺季，較容易的水分來源與多樣化的水分利用方式決定了其更高的LE過程，并主要以蒸散發的形式消耗能量。隨著綠洲感熱通量變小，潛熱通量變大，濕度增強，對綠洲的保護范圍加大，所以增加地面植被的覆蓋程度，對綠洲應對災害天氣有一定的減輕作用。

4)研究區由于降水量較少，生態系統簡單而脆弱，實際蒸散發波動一定會引起區域水熱資源的多層次變化。一方面，實際蒸散發可減小輻射向感熱的轉化，增加空氣濕度，提高最低氣溫及降低最高氣溫，起到調節氣候的作用；但另一方面，則會增加流域農業需水量，促使農墾區增加灌溉，進而造成精河徑流減少，加速湖泊面積萎縮，加劇流域的干旱和荒漠化程度，使原本脆弱的生態環境更加惡化。因此，有必要開展保護生態環境，加強生態管理，促進生態系統良性發展。