楊亮彥，黎雅楠*，范鴻建，王雅婷

（1.陜西省土地工程建設(shè)集團有限責(zé)任公司，西安 710075；2.陜西地建土地工程技術(shù)研究院有限責(zé)任公司，西安 710021；3.山東省城鄉(xiāng)規(guī)劃設(shè)計研究院有限公司，濟南 250013）

蒸散發(fā)（evapotranspiration，ET）主要包括水面和土壤蒸發(fā)、植被蒸騰[1]，在地表水循環(huán)和能量平衡中具有重要的作用。全球大約60%的地表降水通過蒸散發(fā)的方式返回大氣，在干旱區(qū)高達90%[2]。同時蒸散發(fā)直接影響土壤濕度、地表溫度[3]、地下水深度和植被生長狀況[4]，是氣候變化[5]、水文過程[6]及農(nóng)作物估產(chǎn)[7]等的重要參數(shù)。隨著生態(tài)環(huán)境退化、水資源短缺等問題的出現(xiàn)，明晰ET 的時空變化和影響因素具有重要的研究意義。

由于地表異質(zhì)性，傳統(tǒng)的地表蒸散發(fā)量獲取方法以氣象站提供的水面蒸發(fā)量為主，但該數(shù)據(jù)不能有效地反映區(qū)域尺度的地表蒸散發(fā)量，無法滿足區(qū)域水資源合理開發(fā)利用的需要[8]。遙感技術(shù)的發(fā)展為地表蒸散發(fā)量的獲取提供了新的方法，利用遙感數(shù)據(jù)可估算出凈輻射通量、土壤熱通量和顯熱通量，結(jié)合能量平衡公式得到潛熱通量。國內(nèi)外學(xué)者在此基礎(chǔ)上相繼提出了SEBAL（surface energy balance algorithm for land）、SEBS（surface energy balance system）、SEBI（surface energy balance index）等[9]ET 模型，以此獲取地表實際蒸散發(fā)量。基于區(qū)域ET反演模型，眾多學(xué)者及團隊研發(fā)了地表蒸散發(fā)遙感數(shù)據(jù)產(chǎn)品，Miralles等[10]基于Priestley-Taylor 方程研發(fā)GLEAM 產(chǎn)品（0.25°，Daily）；NASA 基于MODIS 衛(wèi)星遙感數(shù)據(jù)，利用彭曼（Penman Monteith，P-M）公式及多種阻力算法制作全球陸地蒸散發(fā)量產(chǎn)品MOD16（500 m，8 d）[11]；美國地 質(zhì)調(diào)查 局（United States Geological Survey，USGS）利用陸面能量平衡模型獲取實際蒸散量的SSEBop產(chǎn)品（1 km，month）[12]；Chen 等[13]基 于SEBS 模型獲得了ET-EB產(chǎn)品；NASA 基于陸面模型和數(shù)據(jù)同化技術(shù)生成全球陸面通量數(shù)據(jù)GLDAS（0.25°，month）[14]。其 中，MOD16 產(chǎn)品因時空分辨率高、監(jiān)測精度高、數(shù)據(jù)容易獲取等優(yōu)勢，被廣泛應(yīng)用于典型地表ET時空分布研究。張?zhí)氐萚15]利用P-M 公式在澴河流域驗證了MOD16 數(shù)據(jù)的精度，發(fā)現(xiàn)均相對誤差（mean relative error，MRE）為7.5%，可滿足區(qū)域ET 時空變化研究需求。Ma 等[16]基于MODIS 數(shù)據(jù)探究了2000—2015 年中國黃土高原ET 變化，結(jié)果表明，ET 在2000—2010 年有明顯增加，但2011—2015年沒有明顯變化趨勢。郭曉彤等[17]基于MOD16數(shù)據(jù)采用復(fù)相關(guān)關(guān)系探究了淮河流域ET 時空變化的驅(qū)動力，發(fā)現(xiàn)在淮河流域，人類活動是ET 時空變化的主要驅(qū)動力。

毛烏素沙地是中國西北典型的農(nóng)牧交錯區(qū)，生態(tài)環(huán)境十分脆弱。由于該地區(qū)降水稀少、氣候干燥、太陽輻射強、蒸發(fā)強烈等自然因素，毛烏素沙地長期存在水資源短缺問題，土地荒漠化進程迅速，對西部和京津冀地區(qū)的生態(tài)環(huán)境造成嚴重危害，因此毛烏素沙地一直是中國干旱的熱點研究區(qū)域。但前人的研究多集中在毛烏素沙地植被覆蓋變化、荒漠化、土地利用變化等方面，在蒸散發(fā)方面以ET 反演模型估算局部地區(qū)日蒸散發(fā)量為主，鮮有針對毛烏素沙地降水量和實際ET時空變化的宏觀研究。本研究基于MOD16 ET 數(shù)據(jù)集，利用趨勢分析法明晰毛烏素沙地2001—2019年ET 時空變化特征，結(jié)合氣象數(shù)據(jù)探究ET 與水熱條件的驅(qū)動力關(guān)系，并分析土地利用∕覆被（land use∕land cover，LULC）和植被指數(shù)對ET時空變化的影響，為落實毛烏素沙地水資源管理及開發(fā)利用制度和生態(tài)環(huán)境保護政策提供參考。

1 材料與方法

1.1 研究區(qū)概況

毛烏素沙地地處中國西北部、黃土高原腹地，橫亙榆林市北部、鄂爾多斯市南部及鹽池縣和靈武 市（37.45° —39.37°N，107.67° —110.67°E）（圖1）。毛烏素沙地處于干旱與半干旱過渡區(qū)，是我國典型的農(nóng)牧交錯帶，其西北部以沙地和草地為主，東南部以耕地和草地為主。研究區(qū)以溫帶大陸性氣候為主，干旱少雨且分布不均，降水多集中在6—8月，年降水量250～440 mm，由西向東南方向遞增。毛烏素沙地屬于西北典型的生態(tài)脆弱區(qū)，是我國沙漠化最嚴重的地區(qū)之一，也是我國防風(fēng)固沙重要的生態(tài)功能區(qū)之一[18-19]。

圖1 毛烏素沙地區(qū)位Fig.1 Distribution of land use types in Mu Us Sandy Land

1.2 數(shù)據(jù)來源及處理

本研究采用的數(shù)據(jù)主要包括ET、氣象、LULC和植被指數(shù)，其中ET 數(shù)據(jù)為MOD16 產(chǎn)品，由NASA 基于P-M 模型計算的數(shù)據(jù)集，空間分辨率為500 m，時間分辨率為8 d。選用2000—2019 年MOD16A3GF 數(shù)據(jù)集（https:∕∕search.earthdata.nasa.gov∕），由MOD16 產(chǎn)品全年統(tǒng)計合成的L4 級產(chǎn)品，時間分辨率為年。本研究利用MODIS Reprojection Tool（MRT）工具對HDF 類型數(shù)據(jù)進行格式轉(zhuǎn)換、投影轉(zhuǎn)換、裁剪、鑲嵌等預(yù)處理。歸一化植被指數(shù)（normalized difference vegetation index，NDVI）為MOD13數(shù)據(jù)（https:∕∕ladsweb.modaps.eosdis.nasa.gov），其空間分辨率為500 m，時間分辨率為16 d，每月可獲取2 期影像，取二者平均值作為當月植被指數(shù)，并利用最大值合成法獲取每年NDVI數(shù)據(jù)。

氣象數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)（https:∕∕data.cma.cn∕），研究區(qū)內(nèi)共有8個國家氣象站點，分別為橫山站（編號53740）、靖邊站（編號53735）、定邊站（編號53725）、鹽池站（編號53723）、神木站（編號53651）、榆林站（編號53646）、鄂托克旗站（編號53529）和烏審旗站（編號53644）（圖1）。氣象數(shù)據(jù)包含地表溫度、氣溫、降水量等8種日值數(shù)據(jù)，本研究通過統(tǒng)計分析獲取每個站點的年平均氣溫和降水量數(shù)據(jù)，并計算毛烏素沙地7個氣象站點的年均氣溫和降水量。LULC 柵格數(shù)據(jù)為30 m 空間分辨率 的GlobeLand30（http:∕∕www.globallandcover.com∕），由中國研制，該數(shù)據(jù)采用WGS-84 坐標系，包括耕地、林地、草地、灌木地、建設(shè)用地、未利用地等10 個一級LULC。目前已發(fā)布2000、2010 和2020版。

1.3 評價指標

1.3.1 ET 變化趨勢 趨勢分析可模擬每個像元多年ET 的變化趨勢，分析ET 在研究時段的時空變化規(guī)律。計算整個研究區(qū)ET 長時間序列的像元回歸斜率K，若斜率為正則表示ET 隨時間變化呈上升趨勢，為負則表示ET 隨時間呈下降趨勢，且斜率的絕對值越大，說明ET 變化趨勢越明顯。其計算公式如下。

式中，K為ET 變化趨勢斜率；n為監(jiān)測時間段的年數(shù)；ETi表示第i年的ET。

1.3.2 相對變化率 相對變化率是相對值對時間的導(dǎo)數(shù)，主要用于探究研究時間段內(nèi)因素的波動變化狀況。

式中，r為相對變化率為研究時段ET 多年平均值。

1.4 數(shù)據(jù)分析

用分析軟件EXCEL 和SPSS 進行數(shù)據(jù)統(tǒng)計和相關(guān)性分析。相關(guān)系數(shù)r取值范圍-1～1，r為正表示參數(shù)為正相關(guān)關(guān)系，r為負表示參數(shù)呈負相關(guān)關(guān)系，絕對值越大相關(guān)程度越高。

2 結(jié)果與分析

2.1 毛烏素沙地ET時空分布

2.1.1 毛烏素沙地ET 空間變化特征 圖2 為毛烏素沙地2000—2019年ET逐年空間變化趨勢，可以看出，毛烏素沙地ET分布具有較強的空間差異性，整體呈西北低、東南高的空間規(guī)律，并且不同行政區(qū)域內(nèi)ET分層明顯，陜西省明顯高于內(nèi)蒙古自治區(qū)，這一特征可能與LULC和地表植被覆蓋度有關(guān)。毛烏素沙地2000—2019年ET變化明顯，除研究區(qū)西部的靈武市農(nóng)業(yè)區(qū)ET較高外，20年內(nèi)研究區(qū)ET在不斷上升，其中在2000—2009年，ET在300 mm·a-1以下，2010年后高于300 mm·a-1的區(qū)域明顯增多，主要分布在陜西省境內(nèi)的定邊縣、靖邊縣、神木市、榆陽區(qū)和橫山區(qū)。

2.1.2 毛烏素沙地ET 時間變化特征 為進一步分析不同時期ET的分布規(guī)律，統(tǒng)計了研究區(qū)2000、2005、2010、2015 和2019 年ET 分布區(qū) 間及占 比（表1），可以看出，2000和2005年，研究區(qū)ET主要集中在250 mm·a-1以下，分別占研究區(qū)的99.0%和89.5%，比例略微下降；到了2010 逐漸增加至250～300 和300～350 mm·a-1，分別占研究區(qū)的55.2%和23.4%；2015 年，研究區(qū)ET 出現(xiàn)回落，主要集中在200～300 mm·a-1之間，占研究區(qū)的74.6%；2019 年毛烏素沙地ET主要區(qū)間恢復(fù)到250～350 mm·a-1，其中300～350 mm·a-1首次成為ET 分布最多的區(qū)間。綜上所述，毛烏素沙地年ET有上升趨勢，且存在地域性差異，毛烏素沙地ET變化受行政區(qū)劃生態(tài)環(huán)境政策和自然因素共同影響。

表1 毛烏素沙地ET分布區(qū)間及占比Table 1 Distribution range and proportion of ET in the Mu Us Sandy Land

從圖3 可以看出，2000—2019 年研究區(qū)年均ET 在179.7～321.4 mm·a-1，多年平 均ET 為258.8 mm·a-1。研究時段內(nèi)，高于和低于平均水平的年份分別有9 和11 個，其中高于平均水平的年份集中分布在2010—2019 年，低于平均水平的年份分布在2000—2009 年，距平相對變化率最大的年份為2000 和2016 年，ET 分別為179.7 和321.4 mm·a-1。以上數(shù)據(jù)表明，2000—2019 年毛烏素沙地ET 在波動中呈明顯增加的趨勢，且2010年之后波動范圍變小。

圖3 2000—2019年毛烏素沙地ET年際變化及距平相對變化率Fig.3 Interannual change and relative change rate of ET in the Mu Us Sandy Land from 2000 to 2019

2.2 毛烏素沙地ET時空變化趨勢

根據(jù)毛烏素沙地ET 年際變化將研究時段分為2 個階段：2000—2009 和2010—2019 年，并將ET 的趨勢變化率分為7 個等級（圖4，表2）?？梢钥闯?，毛烏素沙地各階段ET 存在一定的差異性。在2000—2009年，毛烏素沙地ET平均變化速率為4.60 mm·a-1，呈快速上升趨勢，ET增長速率集中在2～10 mm·a-1，占研究區(qū)的96.86%，其 中K＞5.00 mm·a-1的區(qū)域主要分布在研究區(qū)的周邊，下降的區(qū)域僅占0.04%，主要位于靈武市的農(nóng)業(yè)區(qū)；2010—2018年，毛烏素沙地ET的平均變化速率為5.71 mm·a-1，ET 繼續(xù)保持增加趨勢的同時，速率有所加快，與2000—2009年相似，ET的增長速率集中在2～10 mm·a-1，但K＞5 mm·a-1區(qū)域成為毛烏素沙地最主要區(qū)域，主要分布在研究區(qū)的中部和東部。整體分析，2000—2019 年毛烏素沙地ET 以上升為主，平均變化速率為6.87 mm·a-1，其中2～5 mm·a-1的區(qū)域占研究區(qū)的14.89%，分布在西北部；5～10 mm·a-1的區(qū)域占研究區(qū)的79.75%，主要分布在中部、南部和東部區(qū)域；K＞10 mm·a-1的區(qū)域占研究區(qū)的5.25%，散落在研究區(qū)的東部和南部；基本不變和下降的區(qū)域僅占0.11%。

表2 2000—2019年毛烏素沙地ET的趨勢變化率不同等級的比例Table 2 Proportion of different grades of ET trend change rate in the Mu Us Sandy Land from 2000 to 2019

圖4 2000—2019年毛烏素沙地ET變化速率空間分布Fig.4 Spatial distribution of ET change rate in the Mu Us Sandy Land from 2000 to 2019

2.3 ET影響因素分析

2.3.1 氣候變化對毛烏素沙地ET 時空變化的影響 氣候變化是影響區(qū)域水熱條件空間分布的重要因素，ET 的變化與區(qū)域水熱變化規(guī)律密切相關(guān)[20]。在區(qū)域研究中，年均氣溫表征區(qū)域熱力條件，降水量表征水分條件，研究氣溫和降水量的年際變化有助于探究區(qū)域ET 的影響因素。圖5 為2000—2019 年毛烏素沙地ET 與氣溫、降水量年際變化。研究區(qū)年均降水為337.8 mm，降水年際波動較大，上升趨勢不明顯，其中烏審旗站點2001—2010 年年降水量為383.1 mm，2011—2019年年降水量為430.1 mm，具有明顯上升趨勢；毛烏素沙地氣溫年際變化相對較小，年均氣溫在9.3 ℃左右。研究區(qū)ET 與降水量的年際變化趨勢基本一致，呈波動上升趨勢，除2019 年外，二者的相關(guān)系數(shù)為0.74；與氣溫變化趨勢無明顯相關(guān)性。以上數(shù)據(jù)表明，毛烏素沙地ET 受水熱條件的影響，其中降水量為主要驅(qū)動因素。

圖5 2000—2019年毛烏素沙地ET與氣溫、降水量年際變化Fig.5 Interannual changes of ET，temperature and precipitation in the Mu Us Sandy Land from 2000 to 2019

2.3.2 LULC 對毛烏素沙地ET 時空變化的影響毛烏素沙地LULC 以草地和耕地為主，其次為裸地、灌木地、林地、建設(shè)用地和水體。以2010 年為例（圖6），草地、耕地、裸地、灌木地、建設(shè)用地、林地和水體及濕地占比依次為73.05%、14.63%、9.59%、1.73%、0.50%、0.26%和0.26%，草地遍布整個研究區(qū)，耕地分布在東部和南部，未利用地主要分布在研究區(qū)中部，呈條帶狀分布，其他LULC 散落分布在研究區(qū)內(nèi)。由于MOD16 數(shù)據(jù)缺乏水體ET數(shù)據(jù)，本研究不再分析水體ET變化。

圖6 毛烏素沙地土地利用∕覆被分布Fig.6 Distribution of LULC in the Mu Us Sandy Land

圖7 所示為2000—2019 年毛烏素沙地不同LULC 面積及ET 變化。2000—2019 年間，建設(shè)用地和耕地面積增加，草地、裸地面積減少，未利用地面積變化不大。6 種LULC 的ET 均呈顯著增加的趨勢，與2000年相比，2019年耕地、林地、草地、灌木地、建設(shè)用地和未利用地ET 分別增加70.33%、81.43%、64.52%、57.50%、43.51% 和58.28%。不同LULC的平均ET從大到小依次為耕地＞林地＞灌木地＞建設(shè)用地＞草地＞裸地，耕地和建設(shè)用地面積的增加，裸地、草地面積的減少顯著增加了毛烏素沙地ET，與研究區(qū)ET 增加的趨勢相符。綜上所述，LULC 是影響毛烏素沙地ET 變化的重要因素。但同時發(fā)現(xiàn)2019 年各LULC 的ET普遍高于2000年，說明其他因素對毛烏素沙地ET變化也具有重要影響。

圖7 2000—2019年毛烏素沙地不同土地利用∕覆被面積及ET變化Fig.7 Changes of different LULC areas and ET in the Mu Us Sandy Land from 2000 to 2019

2.3.3 NDVI 對蒸散發(fā)時空變化的影響 為了更好地探究NDVI 對ET 的影響，本研究將毛烏素沙地多年平均NDVI 和ET 分為6 級（圖8）。毛烏素沙地多年平均NDVI 為0.33，8.56%的區(qū)域NDVI≤0.20，主要分布在西北部和腹部；僅有1.25%的區(qū)域NDVI＞0.60，分布在靈武市農(nóng)業(yè)區(qū)和靖邊、定邊境 內(nèi)；0.20≤NDVI＜0.60 的區(qū)域 占研究區(qū)的90.19%，遍布整個研究區(qū)（表3）。毛烏素沙地NDVI整體呈西北到東南逐漸升高的趨勢，與研究區(qū)ET空間分布較為一致，表明ET與NDVI具有較好的相關(guān)性。為了定量化分析二者的關(guān)系，圖9為NDVI 與ET 年際變化曲線，二者變化曲線具有高度的一致性。圖10 為ET 與NDVI 多年均值的全區(qū)域散點圖（紅色代表點密度大），二者決定系數(shù)為0.628 8，進一步證明了ET 與NDVI 的相關(guān)性。整體來說，ET 的變化情況與植被生長狀態(tài)具有高度的相關(guān)性，植被的生長趨勢是影響ET的重要因素。

圖8 2000—2019年毛烏素沙地多年平均NDVI和ET空間分布Fig.8 Spatial distribution of multi-year average NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

圖9 2000—2019年毛烏素沙地NDVI與ET年際變化Fig.9 Interannual change of NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

圖10 2019年ET與NDVI多年均值的全區(qū)域熱點圖Fig.10 Region-wide heat map of the multi-year average of ET and NDVI in 2019

表3 2000—2019年毛烏素沙地多年平均NDVI和ET不同區(qū)間占比Table 3 Proportion of different intervals of the multi-year average NDVI and ET in the Mu Us Sandy Land from 2000 to 2019

3 討論

ET 作為干旱區(qū)水循環(huán)的重要途徑，明晰其時空變化規(guī)律對毛烏素沙地水資源合理開發(fā)及可持續(xù)發(fā)展具有重要的作用。毛烏素沙地ET 在時空變化上存在一定的規(guī)律性，在空間分布上具有明顯的分異性，呈西北向東南逐漸增加的趨勢，其中陜西行政區(qū)內(nèi)的ET明顯高于內(nèi)蒙古自治區(qū)的，與劉靜等[18]研究結(jié)果一致。ET與降水量年際變化趨勢較為一致，與氣溫?zé)o明顯相關(guān)性，表明在干旱半干旱區(qū)，降水量是ET的主要驅(qū)動因素，這與毛烏素沙地缺水的自然特征相符。在不考慮水體的情況下，耕地和林地是毛烏素沙地ET 最大的2 種LULC，與王卓月等[8]的研究結(jié)果一致。多年平均ET與NDVI空間分布一致，年際變化曲線相似，表明毛烏素沙地植被覆蓋及長勢是ET的主要因素，與谷佳賀在黃河流域的研究結(jié)果相符[1]。

2000—2019 年毛烏素沙地ET 增加明顯，其增加速率略高于降水量的變化速率，考慮其原因是耕地、林地面積的增加導(dǎo)致ET 的增加，過度的植樹造林存在惡化當?shù)厣鷳B(tài)環(huán)境的可能性，因此如何平衡植樹造林工程對研究區(qū)的影響是未來需要考慮的問題。本研究發(fā)現(xiàn)，在不同LULC 中，灌木地和草地的ET明顯低于林地和耕地的，因此增加灌木地和草地面積可有效減少水分流失，調(diào)整林木類型是改善植樹造林政策合理性的重要選擇。

本文針對2000—2019 年毛烏素沙地時空特征及其與水熱條件、LULC 和NDVI 的關(guān)系進行了深入研究，但研究仍存在一定的局限性：MOD16數(shù)據(jù)空間分辨率為500 m，且水體部分的ET 缺失，較粗的空間分辨率和缺失值的問題導(dǎo)致部分研究信息無法獲取，研究缺乏完整性；由于研究區(qū)內(nèi)氣象站點數(shù)量較少，且氣象站周邊缺乏ET 數(shù)據(jù)，因此在時空立方體和新興熱點分析[21]等研究方法上受到一定的局限性；在探究相關(guān)驅(qū)動機制時，不同行政區(qū)域的土地政策不同，間接影響研究區(qū)ET的變化，本文僅從宏觀角度對ET進行分析，并未區(qū)分行政區(qū)域解讀相關(guān)土地利用政策，這將是日后工作重點研究的方向。