韓宇平，徐 丹，黃會平，賈冬冬，裴 鐸，張俊峰

(1.華北水利水電大學水資源學院，河南 鄭州 450046；2.河南省黃河流域水資源節(jié)約集約利用重點實驗室，河南 鄭州 450046；3.華北水利水電大學建筑學院，河南 鄭州 450046)

近年來，氣候及植被覆被變化引起的水文效應(yīng)顯著增強，進而影響到區(qū)域內(nèi)蒸散發(fā)，并引起區(qū)域儲水量變化[1-2]。蒸散發(fā)(Evapotranspiration，ET)包括土壤、水面蒸發(fā)和植被蒸騰，是地表水文循環(huán)過程的重要組成部分，對保持區(qū)域能量及水量平衡具有重要意義[3-4]。黃土高原是中國典型的生態(tài)脆弱帶，絕大部分屬于干旱半干旱地區(qū)，1999年以來開始實施退耕還林(草)工程，目前植被在很大程度上得到恢復，植被覆蓋度大幅增加[5-6]。植被的生長會改變地表粗糙度和地表反射率，影響匯流產(chǎn)生過程及蒸騰過程，使區(qū)域ET發(fā)生變化，進而對區(qū)域生態(tài)水文過程產(chǎn)生影響[7-9]。因此，研究變化環(huán)境下ET的時空分布規(guī)律，分析并探討其與各影響因子的關(guān)系，對于了解區(qū)域水文條件、提高水資源利用效率有一定的指導意義[10]。

遙感數(shù)據(jù)產(chǎn)品具有多時間尺度、精度高及監(jiān)測范圍廣的特點，可提供多時間、空間尺度、長時間序列的精準數(shù)據(jù)[11]。目前遙感已成為研究區(qū)域ET、植被覆蓋度、區(qū)域水量存儲等的重要手段[12]?；赑enman-Monteith公式研究發(fā)布的全球陸地蒸散發(fā)數(shù)據(jù)集MOD16已通過全球通量塔數(shù)據(jù)驗證，模擬精度可達86%[13]，被廣泛應(yīng)用于各區(qū)域的蒸散發(fā)研究中。郭曉彤等[14]基于MOD16/ET數(shù)據(jù)產(chǎn)品，分析了淮河流域2000—2014年蒸散發(fā)的時空分布特征及對土地覆被類型變化、氣候變化的響應(yīng)。何韶陽等[15]選擇了華北平原PML_V2、SSEBop_V4和MOD16A2 3種高分辨率ET產(chǎn)品，對不同產(chǎn)品進行精度驗證和時空對比，以期選擇出適合研究區(qū)域的數(shù)據(jù)產(chǎn)品。黃土高原因地表覆被變化劇烈成為研究的熱點地區(qū)，韓盟偉等[16]基于Penman-Monteith模型計算潛在蒸發(fā)量，研究黃土高原1959—2015年ET時空分布規(guī)律及與氣象因子的相互關(guān)系。Ma等[17]利用黃土高原2001—2015年地表反射率(MOD09GA)、地表溫度(MOD11A1)數(shù)據(jù)，基于能量平衡方程反演ET，并進一步探討了ET年際時空變化及與氣象因子、植被變化的關(guān)系。周志鵬等[18]基于MOD16A2數(shù)據(jù)產(chǎn)品，分析了黃土高原全年及四季ET時空分布特征。退耕還林還草工程使黃土高原地區(qū)大量坡耕地轉(zhuǎn)換為耕地和林地，地表植被度提高25%左右，隨著植被覆蓋度提高，區(qū)域ET顯著增加[19-20]。大量研究已表明植被生長發(fā)育對ET有明顯影響，但已有研究較少從植被生長的不同階段來探討ET的時空變化特征及影響因子，歸因分析中也缺乏綜合考慮各因素在空間上對ET的影響。因此本文采用MOD16A2 ET、MOD13Q1 NDVI和氣象數(shù)據(jù)，分析黃土高原地區(qū)蒸散發(fā)在全年、植被生長季和非生長季3個時間尺度的時空變化特征，并在此基礎(chǔ)上研究NDVI、氣象要素對蒸散發(fā)的綜合影響。研究成果有助于掌握黃土高原3個不同時間尺度蒸散發(fā)的時空規(guī)律及影響因子，以期更好地探討植被生長發(fā)育與ET的關(guān)系，為黃土高原生態(tài)恢復和水資源綜合治理、高質(zhì)量發(fā)展提供理論依據(jù)。

1 材料與研究方法

1.1 研究區(qū)概況

黃土高原(33°43′～41°16′N，100°54′～114°33′E)位于中國西北部、黃河中上游地區(qū)，高程為81～4 945 m，總面積約64.9萬km2，主要包括山西、陜西、甘肅、青海、寧夏、河南和內(nèi)蒙古自治區(qū)7個省(區(qū))的部分地區(qū)。該地區(qū)為溫帶大陸性季風氣候，全年降水主要集中在7—9月，年降水量在150～750 mm之間[21]，降水量由東南向西北減少。1999年起，國家開始實行大規(guī)模的退耕還林還草工程，黃土高原為首批試點區(qū)域之一。

1.2 數(shù)據(jù)來源

ET(MOD16A2)和NDVI (MOD13Q1)數(shù)據(jù)來源于美國國家航空航天局NASA網(wǎng)站(https://ladsweb.modaps.eosdis.nasa.gov/ search/)，時間序列為2001—2020年，時間分辨率分別為8、16 d，空間分辨率分別為500、250 m。氣象數(shù)據(jù)(降水、溫度、風速)來源于中國氣象科學數(shù)據(jù)共享網(wǎng)(http://cdc.nmic.cn/home.do)，時間序列為2001—2020年。Landsat-8影像數(shù)據(jù)來源于美國地質(zhì)勘探局(USGS)(https://earthexplorer.usgs.gov/)，時間序列為2013—2020年，空間分辨率為30 m。

首先利用MRT(MODIS Reprojection Tool)軟件對下載的數(shù)據(jù)進行拼接轉(zhuǎn)換等預處理工作，再利用ArcGIS軟件裁剪出黃土高原圖，ET圖像中的無效值賦值為0。黃土高原55個氣象站點的降水、溫度和風速等數(shù)據(jù)，利用創(chuàng)建的泰森多邊形(Voronoi diagram)求平均值，采用克里金(Kriging)插值得到空間分布。研究區(qū)位置及氣象站點的空間分布見圖1。

1.3 研究方法

1.3.1空間回歸分析

黃土高原植被生長發(fā)育特點為每年4月返青發(fā)芽，10月份枝葉枯黃凋萎。因此，按照植被生長時段規(guī)律劃分為生長季(4—10月)、非生長季(1—3、11—12月)。利用slope趨勢分析[10]，對每個柵格的ET數(shù)值進行線性回歸分析，計算公式如下：

(1)

式中Eslope——統(tǒng)計時間尺度內(nèi)的變化速率，Eslope>0表示20年間ET為增加趨勢，Eslope<0表示ET為減少趨勢，Eslope=0表示沒有變化；N——總年數(shù)，本文為20；Ei——ET第i年的值。

1.3.2Morlet小波分析

Morlet小波分析[22]被廣泛應(yīng)用于水文過程中，是一種揭示周期性變化規(guī)律的研究方法。數(shù)學表達見式(2)：

(2)

1.3.3主成分分析

主成分分析(PCA)法[23]是一種把各變量通過正交變換重新組合成一組線性不相關(guān)變量的多元統(tǒng)計方法。通過降維技術(shù)對多個變量進行處理，在保留原信息的基礎(chǔ)上更加客觀地篩選出主導的影響因子，根據(jù)其特征值及貢獻率確定主成分。本文基于PCA分析方法提取各主成分以反映不同因子的綜合效應(yīng)對ET時空變化的影響。在p=0.05 顯著性水平下，通過方差最大化正交旋轉(zhuǎn)，提取前2個主成分。并根據(jù)主成分初始因子載荷量，計算各線性組合中系數(shù)，確定各因子在綜合得分模型中的系數(shù)并進行歸一化處理[24]，得到各因子指標權(quán)重。

1.4 MODIS數(shù)據(jù)驗證

選取2013—2020年行列號為127/35、130/35及126/33的Landsat-8影像數(shù)據(jù)，每年選擇2幅圖像，每幅影像時間間隔約為半年，基于SEBAL模型，結(jié)合洛川、榆中、五寨3個地面氣象站點的風速、氣溫等氣象數(shù)據(jù)，反演各個氣象站點相應(yīng)時間的ET值并擴展至8 d，利用Pearson相關(guān)系數(shù)法對MODIS ET數(shù)據(jù)進行精度驗證(圖2)。結(jié)果表明，在95%置信區(qū)間相關(guān)系數(shù)為0.88，證明MODIS ET數(shù)據(jù)在黃土高原是適用的，可以進行實際蒸散發(fā)的相關(guān)研究。

圖2 SEBAL模型反演ET與MODIS ET關(guān)系

2 結(jié)果與分析

2.1 ET的時空變化特征

2.1.1ET時間序列特征

a)年際線性變化。2001—2020年黃土高原地區(qū)ET年際變化趨勢見圖3。黃土高原全年、生長季和非生長季ET總體均呈增長趨勢，增速分別為6.780 mm/a(p<0.01)、6.256 mm/a(p<0.01)和0.524 mm/a(p<0.05)。3個時段ET數(shù)值差異較大，但變化趨勢基本一致，其中生長季與全年的波動變化具有高度一致性。ET年均值為279.5 mm，最大值和最小值分別出現(xiàn)在2018年(347.8 mm)和2001年(191.1 mm)；ET生長季均值為225.6 mm，占全年的80.71%，最大值和最小值分別出現(xiàn)在2018年(286.2 mm)和2001年(135.2 mm)；非生長季ET均值為57.8 mm，占全年的19.29%，最大值和最小值分別出現(xiàn)在2017年(71.9 mm)和2010年(45.4 mm)。由此可知，黃土高原ET增速由快到慢依次為全年>生長季>非生長季，且非生長季與全年變化波動的一致性小于生長季與全年變化波動的一致性，其原因可能是全年的蒸散發(fā)主要集中在生長季，因此全年ET主要是受生長季的影響較大，非生長季影響相對較小。

圖3 2001—2020年黃土高原不同時間段ET時間序列變化

b)年際周期性變化。圖4為2001—2020年黃土高原全年、生長季和非生長季3個時段的小波系數(shù)實部和小波方差。小波系數(shù)正值代表ET數(shù)值增大，負值代表ET數(shù)值減小。由圖4a、4c、4e可知，3個時段的系數(shù)實部等值線在4～8、13～16 a時間尺度上比較密集。由圖4b、4d、4f可知3個時段均出現(xiàn)2個峰值，全年及生長季峰值分別出現(xiàn)在6、14 a的時間尺度，其中14 a出現(xiàn)的峰值大，說明全年及生長季時段14 a左右的周期震蕩最強，為ET變化的第一主周期，6 a為ET變化的第二主周期；非生長季2個峰值分別出現(xiàn)為5、14 a，最大峰值為14 a，即非生長季14 a為ET變化的第一主周期，6 a為ET變化的第二主周期。不同時段小波方差大小為全年>生長季>非生長季，即在年尺度上黃土高原ET周期變化特征最明顯。根據(jù)14 a強主周期特征，2020年處于ET值減小階段，預測2020年之后會進入增加階段。

圖4 黃土高原不同時段小波系數(shù)實部圖和小波方差

圖4 黃土高原不同時段小波系數(shù)實部圖和小波方差

c)突變檢測。區(qū)域年均ET的M-K檢驗見圖5a，在2001—2020年，年均ET持續(xù)增加(UF>0)；ET突變點發(fā)生在2010年附近(99%的置信區(qū)間，u0.01=±2.32)，表明2010年后，ET增速加快。圖5b為年均NDVI突變檢驗，研究時段內(nèi)，NDVI持續(xù)增加(UF>0)。NDVI的持續(xù)增長對ET產(chǎn)生正相關(guān)作用，ET在2010年附近發(fā)生突變，原因可能是退耕還林還林工程實施以后，樹木生長需要一定時間，約10年左右樹木成林，因此2010年以后成年樹木對黃土高原ET影響更大，加快ET的增長速度。

圖5 突變檢驗

2.1.2空間變化特征

圖6為2001—2020年黃土高原3個不同時間段的ET空間分布。全年及生長季均由東南向西北遞減，高值位于黃土高原東南部的陜西省部分區(qū)域；東部和南部的山區(qū)、灌溉區(qū)等地次之；低值主要位于黃土高原西北部。非生長季高值主要分布于南部秦嶺山地和西部賀蘭山部分區(qū)域。

圖6 2001—2020年黃土高原不同時段ET空間分布

由表1可知，年尺度上，ET值在200～300 mm之間的面積占黃土高原總面積的31.48%，100～200、300～400 mm的比例分別為23.29%、23.83%，極小區(qū)間(0～100 mm)和極大區(qū)間(600～800 mm)所占比例僅為4.16%、0.21%。生長季100～200 mm之間的面積占研究區(qū)總面積的比例最大為33.59%，200～300、300～400 mm的區(qū)間占比分別為30.90%、19.87%，極小區(qū)間(0～100 mm)和極大區(qū)間(600～700 mm)所占比例僅為9.99%、0.01%。全年和生長季ET在空間上有明顯的分異特征，但其極值占比很小，90%左右的區(qū)域ET為100～500 mm。非生長季ET數(shù)值整體偏小，94.63%的區(qū)域小于100 mm，100～200 mm的區(qū)域僅占5.37%。

表1 黃土高原不同時段ET多年均值占比

為進一步研究ET時間變化趨勢的空間分布特征，對每個像元進行線性回歸分析，得到2001—2020年3個時間段ET變化斜率(圖7)。由圖可知，全年及生長季ET變化空間分布情況較為一致，東南部及東部增長幅度較大；非生長季西部地區(qū)增長幅度較大。3個時段呈增加趨勢的面積不同，全年及生長季增加趨勢面積占比分別為93.34%和93.48%，非生長季為75.33%，非生長季變化幅度遠小于生長季和全年。

圖7 2001—2020年黃土高原不同時間段ET變化斜率

2.2 NDVI及氣象要素時空變化

2.2.1NDVI時空變化特征

由圖8知，全年、生長季和非生長季的NDVI整體呈增長趨勢，增速分別為0.005 /a(p<0.01)、0.007 /a(p<0.01)和0.003 /a(p<0.01)，可見退耕還林(草)政策實施以后，黃土高原地區(qū)生態(tài)環(huán)境改善，植被覆蓋率逐年上升。20年間，不同時段NDVI均值為0.305(全年)、0.387(生長季)和0.189(非生長季)，NDVI最大值均在2020年，其值分別為0.463(生長季)、0.219(非生長季)和0.362(全年)，NDVI最小值均在2001年，其值分別為0.305(生長季)、0.153(非生長季)和0.242(全年)。黃土高原不同時段NDVI的平均值為生長季>全年>非生長季，且3個時段NDVI時間序列變化基本一致。

圖8 2001—2020年黃土高原不同時間段NDVI時間序列變化

受水熱條件的影響，3個時段多年平均NDVI空間分布范圍基本一致，由東南向西北遞減，東南部值最高，其次為東部，西北部值最低(圖9)。西南地區(qū)植被類型主要是針闊葉林、寒溫性灌林和草叢，而西北地區(qū)受水熱條件限制，植被類型多為草原、草甸和高山稀疏植物，因此西北部NDVI值最低[25]。全年、生長季和非生長季NDVI空間分布規(guī)律與其對應(yīng)時段ET空間分布基本一致。

圖9 2001—2020年黃土高原不同時間段NDVI空間分布

2.2.2氣象要素時空變化特征

20年間黃土高原年平均降水量為444.99 mm，生長季平均降水量為406.30 mm，非生長季僅為39.12 mm。由圖10a可知，3個時段的降水量均呈增加趨勢，全年、生長季、非生長季降水量增速分別為3.820 mm/a(p<0.05)、3.295 mm/a(p>0.05)、0.410 mm/a(p>0.05)。圖10b表明全年、生長季和非生長季的平均氣溫分別為9.14 ℃、16.50 ℃、 -1.24 ℃，增長速率分別為0.015、0.010、0.014 ℃/a，但均不顯著(p>0.05)。圖10c表明全年、生長季和非生長季平均風速增速分別為2.14、2.16、2.13 m/s。20年間黃土高原全年、生長季和非生長季平均風速增速分別為0.011、0.010、0.013 m/(s·a)，增長速率較低，但平均風速年際間波動較大。

圖10 黃土高原不同時段降水、溫度、風速變化

由圖11可知，3個統(tǒng)計時段內(nèi)不同氣象要素空間分布規(guī)律一致，降水和溫度由東南向西北遞減，與ET分布規(guī)律基本一致；風速由東向西遞減，東北部風速最高。

圖11 黃土高原不同時段降水(mm)、溫度(℃)、風速(m/s)空間分布

2.3 ET影響因子分析

2.3.1ET影響因子的主成分分析

選擇NDVI、氣溫、降水、風速作為影響因子，運用SPSS 23.0 提取主成分特征值大于等于1.0的主成分(表2)。生長季、非生長季和全年的前2個主成分累計貢獻率分別為79.11%、76.19%和83.60%，這表明前2個主成分基本能夠反映這3個不同時間尺度NDVI與氣象因子對ET的貢獻度。

表2 特征值、貢獻率、累計貢獻率

由表3可知生長季第1主成分特征值為1.958，貢獻率為48.95%，NDVI、降水和風速荷載較大，均呈正相關(guān)；第2主成分特征值為1.207，貢獻率為30.16%，溫度在第2主成分中有較大荷載，且呈正相關(guān)。非生長季中，風速、NDVI和溫度在第1成分中占比較大，特征值為1.874，貢獻率為46.86%，第2成分特征值為1.173，貢獻率為29.33%，降水有較高荷載。在年時間尺度上，第1主成分特征值和貢獻率分別為2.161和54.04%，溫度和風速有較大荷載，第2主成分降水和NDVI為主要荷載，其特征值和貢獻率分別為1.183和29.56%。

表3 旋轉(zhuǎn)成分矩陣

各因子指標權(quán)重計算結(jié)果見表4。生長季時段，風速、氣溫和NDVI所占權(quán)重較大，降水權(quán)重比例較低，其指標權(quán)重大小為風速>氣溫>NDVI>>降水，說明在生長季時段對ET影響有較大貢獻的是風速、氣溫和NDVI。非生長季中，NDVI、降水和氣溫權(quán)重占比較大，風速權(quán)重比例相對較低，其指標權(quán)重大小為NDVI>降水>氣溫>風速。在年尺度上，NDVI、降水和風速權(quán)重比例較高，氣溫權(quán)重比例較低，其權(quán)重大小為NDVI>降水>風速>>氣溫，即在全年中，黃土高原NDVI和降水對ET的貢獻較大，其次為風速，氣溫對黃土高原ET影響并不明顯。

表4 主成分因子賦權(quán)計算

2.3.2各主成分因子空間分布

利用ArcGIS10.2軟件Multivariate中的Principal Component功能，取2個主成分分析各自空間分布(圖12)。生長季第1主成分由西向東、由北向南逐漸增大，南部和東北部最大，即NDVI、降水和風速對東部、南部影響大，對西部、北部影響較小。第2主成分為溫度，對東南部影響較大，對南部影響較小。非生長季第1主成分從西北、西部向東南部逐漸升高，即風速、NDVI和溫度對東南部地區(qū)影響較大，對西北、西部地區(qū)影響較小。第2主成分為降水，高值呈散點狀分布在東北及南部和西南小部分地區(qū)。全年尺度第1主成分溫度和風速由西北向東南逐漸升高，即溫度和風速對ET影響為東南高、西北低。第2主成分由東南、西北部向外逐漸減小，即降水和NDVI對東南和西北地區(qū)影響較大，對東部、西部及東北部地區(qū)影響較小。

圖12 ET各主成分空間分布

3 討論

黃土高原自1999年實施退耕還林(草)工程以來，植被覆蓋度明顯提高，伴隨著全球氣候變化，區(qū)域的實際蒸散發(fā)量及水文條件都發(fā)生了明顯改變。本文研究結(jié)果表明實際蒸散發(fā)以6.256 mm/a的速度顯著增加，年均ET值自東南向西北遞減，這與韓盟偉等[16]、周志鵬等[18]研究結(jié)果基本一致。不同時段NDVI及氣象要素對黃土高原ET的影響程度不同，但均呈正作用。在干旱半干旱的黃土高原地區(qū)，區(qū)域蒸散發(fā)受植被密度和降水影響較顯著[2,26]，區(qū)域NDVI、降水的空間分布規(guī)律與ET高度一致，且在年尺度上黃土高原NDVI和降水對ET的貢獻較大，表明NDVI、降水對ET的影響較大。植被大規(guī)?；謴图皻夂蜃兓觿×苏羯l(fā)，增加了黃土高原水資源需求量，進而影響到水資源承載力[27]。影響ET的因子除NDVI、氣象要素外，還包括海拔[28]等自然因素以及土地利用變化[8]、修建水工建筑物[29]等人類活動的影響。后續(xù)研究將考慮其他因素對實際蒸散發(fā)的影響，并進一步分析實際蒸散發(fā)量變化對區(qū)域水文環(huán)境及陸地水儲量的影響。

4 結(jié)論

基于MOD16A2 ET、MOD13Q1 NDVI和氣象數(shù)據(jù)，采用線性回歸分析、Morlet小波分析和主成分分析等方法，探討了黃土高原2001—2020年ET時空演變特征及NDVI和氣象因子的影響，結(jié)論如下。

a)黃土高原全年、生長季和非生長季ET均增加，增速為全年(6.256 mm/a)>生長季(0.524 mm/a)>非生長季(6.780 mm/a)；3個時段均存在2個主周期，第一主周期為14 a；ET值全年及生長季由東南向西北遞減，90%以上區(qū)域為100～500 mm，非生長季ET小于200 mm。

b)全年、生長季和非生長季NDVI均顯著增加，增速為生長季>全年>非生長季；不同時段年均NDVI由東南向西北遞減，與ET空間分布基本一致。降水量在全年和生長季顯著增加，非生長季增加不明顯，溫度和風速在3個時段內(nèi)均無明顯上升趨勢。

c)生長季時段對ET影響有較大貢獻的是風速、氣溫和NDVI，降水對其影響最??；非生長季中，NDVI、降水和氣溫對ET貢獻較大，風速權(quán)重比例即貢獻率相對較低；全年中黃土高原NDVI和降水對ET的貢獻較大，其次為風速，氣溫對黃土高原ET影響并不明顯。