朱曉倩，金曉媚，張緒財，張 京

(中國地質大學(北京)水資源與環(huán)境學院，北京 100083)

蒸散(Evapotranspiration，ET)是土壤蒸發(fā)和植物蒸騰的總稱。蒸散量是自然水循環(huán)系統(tǒng)中十分重要的組成部分，是干旱內陸地區(qū)地表水和地下水的主要排泄方式，干旱區(qū)地表降雨量的90%以上是通過蒸散的形式返回到大氣[1]。格爾木河是柴達木盆地第二大河，不僅為格爾木地區(qū)社會經濟發(fā)展提供水資源保障，同時也是山前平原區(qū)地下水主要補給源和涵養(yǎng)植被生態(tài)的主要水源。河流自昆侖山流出后，經山前戈壁礫石區(qū)時大量下滲對地下水進行補給，后穿過細土平原區(qū)至洪積扇前緣，地下徑流受巖性影響溢出形成泉集河，最終流經鹽沼平原匯入終端達布遜湖排泄。由于格爾木河沖洪積扇下部地下水埋深淺，則會受強烈蒸發(fā)作用，使地下水位變化受到影響[2]，加之近年來對地下水的開采，徑流匯入的終端達布遜湖面積逐漸縮小。因此對格爾木河流域長時間序列的蒸散量時空分布及其影響因素的分析就顯得尤為必要。雖然前人對格爾木河流域蒸散發(fā)有過研究，但缺乏長時間序列蒸散量的時間與空間變化，無論時間尺度還是空間尺度均難準確反映格爾木河流域的年際與年內變化趨勢。

就蒸散量估算的研究方法而言，傳統(tǒng)的蒸散量估算方法多為基于Penman公式[3]及Penman-Monteith公式[4]估算潛在蒸發(fā)及非飽和下墊面蒸發(fā)。后發(fā)展的傳統(tǒng)實測方法如波文比法、渦度相關法、蒸滲儀法等。波文比法[5-6]簡單且精度較高，但其適用范圍卻因其要求下墊面均勻且無平流影響而受到很大的限制;蒸滲儀法精度最高可達0.01～0.02 mm，但該法所需費用高，無法進行大面積布點測量，還需專業(yè)技術人員進行監(jiān)控。渦度相關法是直接測量亂流通量的方法，較其它實測方法更為精密和可靠，但也有必須設立在平坦、植被單一的下墊面上的缺點，且其維護成本較高[7-8]。而在研究區(qū)以往的水資源評價中也不可避免地涉及到潛水蒸發(fā)問題，多位學者[9-10]基于潛水蒸發(fā)強度與潛水埋深的關系利用阿維里揚諾夫經驗公式對格爾木河流域蒸發(fā)量進行計算，而受包氣帶土質、巖性及氣候等的影響，計算結果誤差較大，同時計算所得結果僅是某時段內的蒸發(fā)量?？傊?，傳統(tǒng)估算蒸散量的方法受下墊面、氣象因素、植被等因素的影響，其結果存在的誤差難以估計，且傳統(tǒng)方法存在難以獲取大尺度時空范圍內參數(shù)的問題，這也使得長時間序列下大尺度范圍的蒸散量時空分布的研究無法準確地進行。

而隨著遙感技術的不斷推廣，基于遙感手段推出了許多地表蒸散發(fā)的估算模型，擴展了傳統(tǒng)方法在時空尺度上的應用范圍，同時減小了傳統(tǒng)方法以點代面引發(fā)的誤差，使得蒸散量的估算結果更為精準。Su Z[11-12]提出地表能量平衡系統(tǒng)(Surface Energy Balance System, SEBS)，其目的在于將衛(wèi)星對地觀測的可見光、近紅外和熱紅外波段資料與大氣模式輸出數(shù)據結合，從而更連貫地估算大氣湍流通量，該模型目前應用較為廣泛。之后多位學者[13-18]將該模型應用于不同的空間及時間尺度，并進行了不同程度的驗證，都得到了較好的蒸散量估算結果，證明了SEBS模型計算蒸散量的合理可行性。同時，利用全球陸面同化系統(tǒng)(GLDAS)生成的氣象數(shù)據較氣象站觀測數(shù)據具有更好的時空連續(xù)性，由GLDAS氣象數(shù)據估算得出的蒸散量結果準確性也相對較高。

因此本文以格爾木河流域山前平原區(qū)為研究區(qū)，應用連續(xù)序列的MODIS數(shù)據及GLDAS氣象數(shù)據基于SEBS模型估算研究區(qū)2001-2016年間16年的蒸散量，并就蒸散量時空分布特征及其影響因素進行分析。通過準確估算研究區(qū)長時間序列的蒸散發(fā)時空分布，客觀評價格爾木河流域可持續(xù)利用量，為研究區(qū)水資源合理開發(fā)利用及生態(tài)環(huán)境保護提供參考。

1 研究區(qū)概況

研究區(qū)位于格爾木河流域的山前平原區(qū)，北臨達布遜湖，南靠昆侖山，地勢南高北低，平均海拔為2 780 m。格爾木流域山前平原區(qū)具有典型的高原大陸性氣候特征，干旱，多風少雨，蒸發(fā)強烈。多年平均降水量為42.7 mm，多年平均蒸發(fā)量約為2 600 mm，年均氣溫4.7 ℃，相對濕度32%。研究區(qū)地形由西向東傾斜且南高北低，地貌復雜，地貌類型由南向北分別為山前洪積平原、沖積平原、沖湖積平原以及湖積平原，地下水的賦存及運動依賴于巨厚的第四系松散沉積物，從而形成了研究區(qū)孔隙地下水系統(tǒng)，由北向南主要可分為鹽沼地區(qū)、細土平原區(qū)及山前戈壁礫石區(qū)。研究區(qū)工業(yè)以鹽化工為主，植被類型主要為枸杞、藜麥、蘆葦、紅柳等。

格爾木河為青海境內河流，發(fā)源于昆侖山北麓，經格爾木市，最終流入東達布遜湖，全長468 km，流域面積1.786×104km2。河水補給來源以高山融雪為主，次為河谷潛水，平均流量21.3 m3/s，年徑流量6.41×108m3。

圖1 研究區(qū)地理位置示意圖Fig.1 Geographical location map of the study area

2 研究方法及數(shù)據

2.1 研究方法

(1)SEBS模型

SEBS(Surface Energy Balance System)模型主要通過地表能量平衡方程估算地表凈輻射通量和土壤熱通量，之后計算出感熱通量，由地表凈輻射通量、土壤熱通量和感熱通量得出潛熱通量及蒸發(fā)比，最終估算得出日蒸散量。

任一時刻，地表的凈輻射通量與該區(qū)域感熱通量、潛熱通量及土壤熱通量處于平衡，由此得任一時刻的表面能量平衡方程為：

Rn=G0+H+λE

(1)

式中：Rn——凈輻射通量/(W·m-2)；

G0——土壤熱通量/(W·m-2)；

H——感熱通量/(W·m-2)；

λE——潛熱通量/(W·m-2)；

λ——水的汽化潛熱/(J·kg-1)；

E——實際水分蒸散總量/(mm·d-1)。

潛熱通量指沿垂直方向，單位時間、單位面積上地表和植被以及大氣之間因為水的相變而交換的熱量，而在地表和大氣的能量交換中，潛熱通量所占的比例即為蒸發(fā)比。其計算公式為：

(2)

式中：Λ——蒸發(fā)比；

Λr——相對蒸發(fā)比；

H——實際蒸散量與最大蒸散量之比。

日蒸散量計算公式為：

(3)

式中：Edaily——地表實際日蒸散量/(mm·d-1)；

λ——汽化潛熱/(J·kg-1)；

ρw——水的密度/(1 000 kg·m-3)。

(2)一元線性回歸分析

回歸分析是確定兩種或兩種以上變量間相互依賴的定量關系的一種統(tǒng)計分析方法。在回歸分析中，只包括一個自變量和一個因變量，且二者的關系可用一條直線近似表示，這種回歸分析稱為一元線性回歸分析。通過一元線性回歸分析可模擬研究區(qū)每個柵格的變化趨勢，本文通過該方法分析研究區(qū)2001—2016年每年7月份蒸散量的變化趨勢。其計算公式為：

(4)

式中：n——總年數(shù)；

MET·i——第i年的ET值，i為年序號(1、2、3……)。

計算得出的Slope可反映研究區(qū)各柵格點ET變化趨勢。若Slope>0，則該柵格點ET呈增長趨勢；若Slope<0，則該柵格點ET呈減少趨勢。

(3)Mann-Kendall顯著性檢測法[19-21]

Mann-Kendall是一種已經廣泛應用于水文數(shù)據分析非參數(shù)檢驗方法，該方法能很好地揭示時間序列的變化趨勢，但是不需要樣本數(shù)據服從一定的分布，同時也不受少數(shù)異常值的干擾。其檢測過程如下：

對于時間序列X(x1、x2、x3……)：

(5)

其中，n為數(shù)據樣本長度，sgn為符號函數(shù)，其定義如下：

(6)

當n<10時，直接使用統(tǒng)計量S進行雙邊趨勢檢驗。

當n≥10時，統(tǒng)計量大致服從正態(tài)分布，使用檢驗統(tǒng)計量Z進行趨勢檢驗，Z值計算公式為：

(7)

其中，

(8)

式中：m——序列中結(重復出現(xiàn)的數(shù)據組)的個數(shù)；

ti——結的寬度(第i組重復數(shù)據組中的重復數(shù)據個數(shù))。

2.2 研究數(shù)據

(1)MODIS數(shù)據

本文選取2001—2016年中等分辨率的MODIS數(shù)據，共計1472景，基于SEBS模型估算地表蒸散量。其中包括空間分辨率分別為500 m、1 000 m，時間分辨率均為8 d的MOD09A1數(shù)據(地表反射率產品)和MOD11A2數(shù)據(地表比輻射率與地表溫度產品)。

MODIS數(shù)據均由美國NASA網站獲取(https://reverb.echo.nasa.gov/reverb)。為保證研究中數(shù)據的一致性，將數(shù)據統(tǒng)一重采樣為500 m分辨率。

(2)GLDAS氣象數(shù)據

由于研究區(qū)氣象站較少，本文采用GLDAS(Global Land Data Assimilation System)數(shù)據中提取的比濕/(kg·kg-1)、風速/(m·s-1)、氣溫/K和氣壓/Pa四個波段參與SEBS模型蒸散發(fā)的計算。選取2001—2016年時間分辨率為1月、空間分辨率為0.25°的GLDAS氣象數(shù)據，共計192景。GLDAS數(shù)據來源于https://disc.gsfc.nasa.gov/datasets。為保證研究中數(shù)據的一致性，數(shù)據需重采樣為500 m分辨率。

(3)DEM數(shù)據

為定量地反映研究區(qū)地形的空間變化特征，本文選取了空間分辨率90 m的SRTM-DEM數(shù)據，并于ENVI軟件中進行空間投影、重采樣等預處理，使其與MODIS、GLDAS數(shù)據空間分辨率一致。數(shù)據來源于中國科學院計算機網絡信息中心國際科學數(shù)據鏡像網站(http://www.gscloud.cn)。

(4)氣象數(shù)據

為研究氣象因素對蒸散量的影響，本文選取了研究區(qū)內格爾木氣象站氣溫、降水、相對濕度、蒸發(fā)等氣象數(shù)據。數(shù)據來源于 “中國氣象數(shù)據網(http://data.cma.cn)”的中國地面氣候資料日值數(shù)據集。

3 結果與分析

3.1 蒸散量時空分布特征

本文基于SEBS模型結合GLDAS氣象數(shù)據得到格爾木河流域山前平原區(qū)2001-2016年的日蒸散量，乘以相應月份的天數(shù)得到月蒸散量，將各月蒸散量相加即為年蒸散量。通過研究區(qū)2016年蒸散量空間分布圖(圖2)可知：研究區(qū)北部為大片鹽池，蒸散量最高，最高可達780 mm；研究區(qū)內蒸散量與河流分布特征具有較好的一致性，自南向北的河流流經區(qū)域蒸散量整體較大；中南部格爾木市東西兩側為研究區(qū)植被主要發(fā)育帶，市西側為格爾木西農場，以耕地為主，農作物類型主要為枸杞、藜麥、小麥等，市東側除農田主要發(fā)育蘆葦、紅柳、沙柳、芨芨草等，蒸散量較高，約200～300 mm；中北部裸土地區(qū)和南部荒漠地區(qū)蒸散量較小，約100 mm。

圖2 研究區(qū)2016年蒸散量空間分布圖Fig.2 Annual ET distribution of the study area in 2016

對研究區(qū)2001—2016年蒸散量年際變化及多年平均蒸散量年內變化進行統(tǒng)計(圖3)，結果表明：研究區(qū)2001—2004年蒸散量略有升高，到2005年有大幅度增長，增長到82.94 mm，總蒸散量約3.79×108m3，2005—2008年蒸散量基本保持不變，浮動于90 mm左右，2008—2016年蒸散量有所波動，但2001—2016年蒸散量總體呈上升趨勢。其中2015年蒸散量最大，為185 mm，總蒸散量約8.45×108m3；2001年最小，為25 mm，總蒸散量約1.16×108m3。研究區(qū)多年平均年內夏季即6、7、8三個月蒸散量最大，約20 mm；冬季蒸散量小，11、12、1、2月蒸散量值在0.5～3 mm之間。通過與前人對格爾木河流域蒸散量估算結果[22～24]對比表明，本文計算所得蒸散量空間分布規(guī)律與前人計算結果吻合，這也進一步證明了SEBS模型估算蒸散量結果的準確可靠性；但由于研究區(qū)范圍不同，蒸散量絕對值上有差異；同時，前人計算區(qū)域蒸散發(fā)時只采用格爾木一個氣象站的氣象值代表全區(qū)，計算結果有一定誤差。本文采用GLDAS全區(qū)氣象數(shù)據計算區(qū)域蒸散發(fā)，所得結果相對更為可靠。

圖3 研究區(qū)2001—2016蒸散量年際變化規(guī)律圖Fig.3 Annual variation of actual evapotranspiration in the study area from 2001 to 2016

3.2 基于線性回歸和Mann-Kendall的蒸散量變化趨勢分析

圖4 研究區(qū)連續(xù)時間序列蒸散量空間變化趨勢(a)及顯著性變化分類(b)圖Fig.4 Spatial variation trend(a) and classification of significant changes(b) of continuous time series evapotranspiration in the study area

為提高分類結果的精確性，將兩種方法結合得到研究區(qū)16年間蒸散量空間變化趨勢圖(圖4)，即當Slope<0且|Z|>1.96,序列呈顯著下降趨勢；當Slope<0且|Z|≤1.96,序列呈下降但不顯著趨勢；當Slope≥0且Slope≤10，序列呈基本穩(wěn)定趨勢；當Slope>10且|Z|>1.96,序列呈顯著增長趨勢；當Slope>10且|Z|≤1.96,序列呈增長但不顯著趨勢。最終將研究區(qū)連續(xù)時間序列蒸散量空間變化趨勢分為五類：顯著增長、輕微增長、基本穩(wěn)定、輕微降低、顯著降低。結果表明：線性回歸分析法和Mann-Kendall顯著性檢測法得出的結果具有良好的一致性。研究區(qū)大部分地區(qū)蒸散量呈穩(wěn)定不變狀態(tài)，其面積占研究區(qū)總面積的68.58%；蒸散量顯著增長的地區(qū)主要集中在研究區(qū)中南部植被發(fā)育帶及北部鹽池區(qū)，占研究區(qū)總面積的29.70%；僅極少地區(qū)有輕微或顯著的減少趨勢，蒸散量減少地區(qū)約占研究區(qū)總面積的1%；其余輕微增長地區(qū)占研究區(qū)總面積的0.72%。同時為分析格爾木河流域山前平原區(qū)各季節(jié)蒸散量變化趨勢，將一年12個月按照春(3—5月)、夏(6—8月)、秋(9—11月)、冬(12—2月)劃分為四季。分別估計四季蒸散量穩(wěn)定時的值，統(tǒng)計了研究區(qū)2001-2016年四季占全年蒸散量的比重，并求其平均值，得到16年間四個季節(jié)在年蒸散量中所占比例值(表1)，結合實際，最終確定四季的Slopesta值：春季Slopesta=1.5，夏季Slopesta=6，秋季Slopesta=3，冬季Slopesta=0.5。

表1 研究區(qū)四季蒸散量占年蒸散量比例表Table 1 Proportion of four seasons evapotranspiration to annual evapotranspiration in the study area

按照確定的slope值得到研究區(qū)16年間四季蒸散量空間變化趨勢圖及顯著性變化趨勢圖(圖5、6)，結合顯著性變化趨勢分類的面積統(tǒng)計結果進行分析。

結果表明：格爾木河流域山前平原區(qū)夏季蒸散量增減幅度遠高于其他三季，秋季和春季次之，冬季增減幅度最小。夏季和秋季均有大面積顯著增長地區(qū)，主要集中在研究區(qū)北部鹽池區(qū)及中南部植被發(fā)育帶，研究區(qū)中北部及南部保持基本穩(wěn)定趨勢，小部分區(qū)域有輕微降低趨勢。冬季研究區(qū)蒸散量除鹽池部分有顯著和輕微的增長外，其余大部分地區(qū)蒸散量都有輕微降低趨勢。春季蒸散量也主要集中在北部鹽池區(qū)，中南部植被發(fā)育帶較夏季、秋季增幅減少，但也有少量顯著增漲趨勢。蒸散量變化趨勢面積統(tǒng)計(圖7)看，呈現(xiàn)顯著增長趨勢：夏季(26.46%)>秋季(25.49%)>春季(20.99%)>冬季(16.66%)；輕微增長趨勢：冬季(4.87%)>春季(4.38%)>夏季(0.49%)>秋季(0.36%)；基本穩(wěn)定趨勢：秋季(73.66%)>夏季(71.54%)>春季(67.79%)>冬季(52.25%)；輕微降低趨勢：冬季(22.49%)>春季(5.19%)>夏季(1.06%)>秋季(0.45%)；顯著降低趨勢：冬季(3.72%)>春季(1.64%)>夏季(0.45%)>秋季(0.04%)。總體上格爾木河流域山前平原區(qū)大部分地區(qū)蒸散量處于基本穩(wěn)定狀態(tài)，16年來北部鹽池地區(qū)不斷擴大導致鹽池區(qū)域蒸散量顯著增加。

圖5 研究區(qū)蒸散量季節(jié)空間變化趨勢圖Fig.5 Seasonal spatial variation trend of evapotranspiration in the study area

圖6 研究區(qū)蒸散量季節(jié)顯著性變化趨勢圖Fig.6 Seasonal significance variation trend of evapotranspiration in the study area

圖7 研究區(qū)全年及四季蒸散量變化趨勢柱狀圖Fig.7 Histogram of annual and seasonal evapotranspiration variation trend in the study area

為進一步分析研究區(qū)蒸散量的時空分布，將研究區(qū)按山前戈壁礫質平原、綠洲平原、鹽沼平原及鹽湖區(qū)劃分為四個地貌及水文地質單元，分別對不同地貌單元蒸散量進行統(tǒng)計分析(圖8)。結果表明：研究區(qū)蒸散量空間分布與地貌類型空間格局具有較好的一致性?？傮w來看，研究區(qū)蒸散量鹽湖區(qū)>鹽沼平原>綠洲平原>山前戈壁礫質平原。四個地貌單元蒸散量從2001—2016年總體呈增長趨勢，其中鹽湖區(qū)增長速率最快。2010年后，由于洪水及耕地面積增加的影響，綠洲平原蒸散量增長速率加快，鹽沼平原和山前戈壁礫質平原蒸散量16年間增長速率較慢。通過對不同地貌四季蒸散量的對比發(fā)現(xiàn)(圖9)，任一地貌類型對應蒸散量隨季節(jié)的變化與研究區(qū)16年間四季蒸散量的空間變化趨勢(圖5)一致，即除鹽湖區(qū)外，其余三個地貌類型蒸散量均呈夏季>秋季>春季>冬季的變化規(guī)律。

圖8 不同地貌單元對應蒸散量變化趨勢Fig.8 Variation trend of evapotranspiration of different geomorphic units

圖9 不同地貌單元對應蒸散量四季變化趨勢Fig.9 Variation trend of evapotranspiration of different geomorphic units in four seasons

4 蒸散量影響因素分析

4.1 氣象因素

本文利用研究區(qū)內格爾木氣象站2001—2016年每月的氣象數(shù)據與對應月的蒸散量進行相關性分析(圖10)。結果表明：研究區(qū)蒸散量與氣溫相關性較好，確定系數(shù)達0.604。氣溫對研究區(qū)內鹽湖區(qū)的影響較大。蒸散量與降水量相關性一般，確定系數(shù)為0.405。通過氣溫、降水與蒸散量的年內與年際變化趨勢對比分析可以看出，蒸散量總體與氣溫、降水量變化趨勢一致，呈正相關關系，即蒸散量隨氣溫和降水量的增大而升高。

圖10 氣象因素對蒸散影響分析圖Fig.10 Effect of meteorological factors on evapotranspiration

4.2 植被指數(shù)

植被蒸騰也是區(qū)域主要蒸發(fā)途徑之一，研究區(qū)內植被主要發(fā)育帶分布于綠洲平原中南部格爾木市東西兩側，因此植被對研究區(qū)蒸散量的影響主要表現(xiàn)在綠洲平原區(qū)。近年來，隨著農田面積的不斷擴大，綠洲平原區(qū)的蒸散量也在不斷增大。為分析植被對研究區(qū)蒸發(fā)的影響，本文將連續(xù)序列的植被指數(shù)MODIS NDVI數(shù)據與計算得出的蒸散量進行相關性分析(圖11)，結果表明：兩者相關性較好，確定系數(shù)為0.612，呈正相關。

圖11 植被指數(shù)與蒸散量相關性分析圖Fig.11 Correlation analysis diagram of NDVI and evapotranspiration

圖12 不同地貌蒸散量與NDVI變化對照圖Fig.12 Comparison of variation of evapotranspiration and NDVI in different physiognomies

為進一步研究植被對研究區(qū)蒸散量的影響，本文利用圖1中AB剖面同時提取2016年7月NDVI及蒸散量空間分布圖中對應點的NDVI及蒸散量值，通過滑動平均法，對比分析不同地貌蒸散量與NDVI值的變化特征(圖12)。結果表明：研究區(qū)由北向南AB剖面的不同地貌下，蒸散量和NDVI值總體變化趨勢具有較好的一致性。除水體月蒸散量較高，最高可達90 mm外，NDVI值較低的裸土地區(qū)和荒漠地區(qū)蒸散量值也相對較低，在0～15 mm之間；NDVI值位于0.3～0.7的灌木地區(qū)，蒸散量也相應增大，波動于15～45 mm間；NDVI值位于0.2～0.4的農田地區(qū)，蒸散量值波動于20～40 mm間。因此，蒸散量與NDVI呈正相關，蒸散量值隨NDVI的增大而增大。

5 結論

(1)研究區(qū)蒸散量2001—2016年總體呈增長趨勢。鹽湖區(qū)蒸散量最大，蒸散量可達780 mm；鹽沼平原和綠洲平原其次，山前戈壁礫質平原蒸散量最小，蒸散量約70 mm。

(2)研究區(qū)內鹽沼平原及山前戈壁礫質平原大部分地區(qū)蒸散量變化在16年間呈穩(wěn)定不變狀態(tài)，鹽湖區(qū)和綠洲平原呈顯著增長趨勢。

(3)研究區(qū)16年間不同地貌蒸散量夏季增減幅度均遠大于其他三季，冬季變化幅度最小。

(4)研究區(qū)蒸散量隨氣溫及降水的增大而增大。植被對蒸散量的影響主要表現(xiàn)在綠洲平原，蒸散量與植被指數(shù)呈較好的正相關。