沂河流域植被覆蓋時空演變及其與徑流的關系研究

石智宇, 趙 清, 王雅婷, 張連蓬

(江蘇師范大學 地理測繪與城鄉規劃學院, 江蘇 徐州 221116)

植被是聯結地表與大氣之間物質、能量交換的關鍵環節，在陸地表面能量交換、水分循環和生物地球化學循環過程中起著至關重要的作用[1]。植被是地理環境的重要組成部分，具有涵養水源、保持水土、調節氣候等作用，植被變化與氣候、土壤、地形條件、人類活動等要素密切相關。徑流是流域水循環的重要環節，徑流量的變化會影響水文系統的演化和水資源的開發利用，進而影響生態環境及社會經濟的發展[2-3]。植被變化可以通過改變冠層截留、植被蒸騰作用和土壤入滲特性等水文過程，顯著地改變地表徑流，進而影響流域水資源可利用量[4-7]。因此，基于植被覆蓋時空演變過程，逐像元分析流域植被覆蓋變化對徑流變化的影響，對揭示流域生態環境演變過程和改善流域水資源開發利用有重要意義。

近些年來，國內外學者圍繞植被覆蓋變化及其對徑流的影響方面開展了大量研究。畢早瑩等[5]基于遙感植被指數和水文、氣象數據，通過基于Budyko假設的彈性系數法，分析了1982—2015年窟野河流域植被變化及其對徑流變化的影響，發現窟野河流域內植被得到了極大的恢復，植被變化成為決定徑流變化的主導因素，其次為氣候變化和人類活動；Kim等[8]選取4個不同冠層的徑流區研究植被冠層對地表徑流的影響，結果表明枯落物層的蓄水能力和葉片形狀影響地表徑流的產生，除此之外，土壤水分條件也是影響地表徑流的重要因素；支童等[9]以陜北禿尾河流域為研究對象，分析流域植被變化特征，探討了降水與植被覆蓋變化對徑流的影響并指出影響徑流縮減的物理機制；李天生等[10]基于Budyko理論具體分析了珠江流域中上游地區氣候和植被變化對該區域徑流的影響，結果發現流域降雨減少是徑流減少的主導因素，而植被在研究時段沒有顯著的變化趨勢，其對徑流減少的貢獻不顯著。

遙感技術因能提供反映不同時空尺度上的植被覆蓋信息及其動態變化，方便快捷，逐漸成為監測植被覆蓋時空演變的主要手段，其中歸一化植被指數(Normalized Difference Vegetation Index，NDVI)對綠色植被信息敏感，可以很好地反映地表植被生長狀況，其變化在一定程度上能代表地表植被變化[11-12]，是目前應用最為廣泛的植被指數之一[13-14]。SCS(Soil Conservation Service)徑流模型是美國農業部水土保持局開發的流域水文模型，由于該模型簡單易行，所需參數較少，對觀測數據要求不高且有效考慮到流域下墊面的特點，被廣泛采用[15-17]。使用SCS模型并結合地理信息技術可以獲取地表徑流的空間分布，從而實現在像元尺度上研究植被覆蓋變化與徑流的關系[15,18]。

沂河流域地處魯中南山地丘陵區，是沂蒙山區核心組成部分和重要的水源涵養區，流域內地形起伏大，生態脆弱，水土流失較嚴重[19]。近年來沂河流域內經濟社會持續發展，生態環境發生較大變化，人水關系趨于緊張，但針對沂河流域的研究多集中在水文變化、水土流失、土地利用變化等方面[19-22]，對于植被覆蓋變化及其水文效應研究相對較少。在人類活動日益加劇和經濟社會快速發展的背景下，沂河流域植被覆蓋時空格局是否發生變化？如果沂河流域植被覆蓋呈現變化趨勢，那么其與地表徑流變化的相關關系如何？這些問題都有待解答。因此本文選取沂河流域作為研究區，基于2000—2020年MOD13Q1 NDVI遙感數據，使用趨勢分析法，分析沂河流域20年間植被覆蓋時空演變特征，在此基礎上，結合SCS模型，使用相關性分析法，在像元尺度上探討植被覆蓋變化與徑流的相關關系，以期為流域生態環境保護和水資源規劃利用提供理論依據，同時為中小型流域植被覆蓋變化和徑流關系研究提供一定參考。

1 研究區概況

沂河是淮河流域沂沭泗水系中主要河流之一，發源于山東省沂源縣西部，流經沂源、沂水、沂南、臨沂、蒙陰、平邑、郯城等縣市，在江蘇省新沂市匯入駱馬湖，河道全長333 km，流域總面積為11 820 km2，是魯南地區和蘇北地區重要的山洪河道，有東汶河、蒙河、祊河等諸多支流匯入。

由于沂河下游開辟新沂河、“分沂入沭”等河道整治工程，流域界限被破壞，因此本文選取沂河臨沂市以上相對完整的流域作為研究區域。沂河臨沂以上河道長223 km，流域面積約9 950 km2，屬暖溫帶半濕潤季風氣候，年均氣溫13～14.3℃，流域地勢西北高、東南低，每年6—9月為汛期，多年平均年降水量830 mm[21]，土壤類型較多，主要有棕壤、褐土、潮土、砂姜黑土等，流域處于暖溫帶落葉闊葉林帶，主要樹種有刺槐、油松、赤松、加拿大楊、旱柳等[19]。降雨是沂河流域徑流補給的主要來源且徑流年內變化十分劇烈，汛期多山洪暴雨，枯水期徑流量很小甚至斷流，導致沂河流域成為為易旱易澇地區，對當地經濟發展帶來不利影響[21,23]。

2 數據與研究方法

2.1 數據來源與處理

本文遙感影像數據采用美國國家航空航天局(NASA)MODIS數據中的MOD13Q1 NDVI數據(http:∥ladsweb.modaps.eosdis.nasa.com)，文件格式為HDF，空間分辨率為250 m，時間分辨率為16 d，數據時間范圍為2000年2月至2020年12月，共480景影像。利用MRT(MODIS Reprojection Tool)軟件對獲取的數據進行投影轉換、裁剪和格式轉換等處理，獲取質量較為可靠的NDVI數據集。采用最大值合成法(Maximum Value Composites，MVC)獲取年最大NDVI作為年NDVI值，能較好地代表當年植被生長最好的狀況[24]，并可進一步消除云、大氣、太陽高度角等的部分干擾[25]。

研究區土地利用數據使用來自于中國科學院資源環境科學與數據中心(http:∥www.resdc.cn)提供的1 km柵格數據，將土地利用類型分為耕地、林地、草地、水域、建設和未利用土地，包括2000年、2010年、2020年三期土地利用數據。

土壤數據來源于聯合國糧農組織(Food and Agriculture Organization of the United Nations，FAO)和維也納國際應用系統研究所(International Institute for Applied System Analysis，IIASA)所構建的世界土壤數據庫(Harmonized World Soil Database，HWSD)，采用FAO 90土壤分類系統，包括上層土壤(0—30 cm土層)和下層土壤(30—100 cm土層)的砂粒、黏粒、碎石等百分比含量及USDA土壤質地分類數據。

DEM數據使用地理空間數據云(http:∥www.gscloud.cn/)提供的SRTMSLOPE 90 m分辨率坡度數據產品，經過ArcGIS裁剪得到研究區坡度數據。降水數據來源于中國氣象數據網(http:∥data.cma.cn)，包含研究區內沂源、沂水、沂南、蒙陰、平邑、費縣和臨沂7個雨量站點2000—2020年降水量數據(表1)，經反距離加權插值法(IDW)獲取與NDVI數據具有相同投影和空間分辨率的柵格數據。

表1 沂河流域雨量站位置信息及數據時間

2.2 研究方法

2.2.1 趨勢分析 本文采用一元線性回歸分析逐像元模擬NDVI年際空間變化趨勢。一元線性回歸分析可以模擬每個柵格的變化趨勢，以單個像元時間變化特征反映整個區域時空格局演變[26]，計算公式為：

(1)

式中:slope為像元NDVI回歸方程的斜率；i為年序號；n代表時間序列長度；NDVIi代表第i年的NDVI值。若slope>0,則表示NDVI隨時間變化呈增加趨勢,若slpoe<0,則表示NDVI隨時間變化呈減少趨勢,擬合斜率絕對值越大,表示變化越明顯。

2.2.2 徑流模擬 集總式水文模型具有結構簡單、參數較少、復雜程度較低等優點，SCS(Soil Conservation Service)水文模型是美國農業部水土保持局于1954年開發的集總式水文模型，實用性強且改進空間大，近年來得到了比較廣泛的應用[27]。該模型綜合考慮了流域降雨、地形、土壤類型、土地利用方式、前期土壤濕潤狀況與徑流的關系，SCS水文模型基本產流方程[28]為：

(2)

(3)

式中:P為降雨量(mm);Q為徑流深(mm);S為最大可能滯留量(mm);CN值是反映降雨前流域特征的一個綜合參數,用于描述降雨—徑流關系。CN值通常由土壤性質、土地利用方式和降雨事件前期土壤濕潤情況等數據來確定[28]。

根據研究區各類型土壤上層與下層土壤質地數據，基于SPAW軟件計算研究區上層土壤與下層土壤最小下滲率并選取較小一組作為該土壤類型最小下滲率，參考SCS徑流模型的水文土壤分類標準[18](表2)劃分流域水文土壤組，在ArcGIS軟件下重分類，得到研究區水文土壤組數據(圖1)。

表2 SCS模型水文土壤組定義指標

圖1 沂河流域水文土壤組

SCS徑流模型考慮前期降水對徑流的影響，將前期土壤濕潤程度根據徑流事件發生前5天的降雨總量(即前期降雨指數API)劃分為干旱(AMCⅠ)、平均(AMCⅡ)和濕潤(AMCⅢ)3種狀態，本文研究中假設前期土壤濕潤程度為平均狀態(AMCⅡ)。

SCS徑流模型中不同的下墊面類型對產匯流的影響不同，本文選取沂河流域2000年、2010年、2020年三期土地利用數據，通過GIS軟件分析三期土地利用(圖2)，得到研究區三期土地利用類型變化，其中2000—2005年、2006—2015年、2016—2020年分別基于2000年、2010年、2020年土地利用數據生成的CN值模擬地表徑流深。

結合研究區概況，在SCS模型提供的CN值查算表基礎上，參考國內有關學者的研究成果[29]，確定正常條件下的CN值。考慮到研究區地形坡度對CN值的影響，本文采用Huang坡度修正公式[30]對CN值進行修正，公式如下：

(4)

式中:CN2為AMCⅡ條件下CN值；CNa為坡度修正后的CN2值;α為多邊形坡度值,用百分比(%)表示。

基于GIS軟件，將研究區水文土壤、土地利用數據兩層數據交叉，得到包含土壤、土地利用雙重信息的數據圖層，輸入修正后CN值，生成AMCⅡ條件下CN值圖層，結合研究區降水數據，利用ArcGIS柵格運算功能得到逐年地表徑流深圖層。

圖2 沂河流域土地利用變化

2.2.3 植被覆蓋變化與徑流的相關分析 為了分析研究區植被覆蓋變化與徑流的相關關系，本文采用相關分析方法，基于像元尺度對2000—2020年沂河流域年最大NDVI與徑流深進行逐像元相關分析，使用相關系數來反映植被NDVI和地表徑流深的相關程度，在p<0.05置信水平對相關性結果進行顯著性分析。相關分析計算公式如下：

(5)

3 結果與分析

3.1 NDVI年際變化特征

提取沂河流域年最大NDVI代表當年植被覆蓋最好情況，并將其作為該年的NDVI值(圖3)，可以看出，2000—2020年沂河流域NDVI整體呈現波動增加趨勢，20年間流域年最大NDVI介于0.67～0.79之間，增長速率為1.46%/10 a(R2=0.118)。植被覆蓋時間變化大致可以分為3個階段：(1) 2000—2002年植被覆蓋處于較低水平，其中2002年出現20 a間NDVI最低值；(2) 2003—2013年植被覆蓋較上一階段出現較大幅度提升，NDVI值整體保持小幅波動上升趨勢；(3) 2014—2020年植被NDVI值較上一階段出現較大程度波動，2014年出現該階段NDVI最低值，該年植被覆蓋情況幾乎退化至2000—2001年水平，但總體上該階段植被覆蓋依然呈現增加趨勢。總體而言，2000—2020年沂河流域NDVI呈現增加趨勢，流域植被覆蓋程度得到了改善，植被生長情況總體向較好方向發展。

圖3 2000-2020年沂河流域NDVI年際變化趨勢

3.2 NDVI空間分布特征及其變化趨勢

在空間分布上，2000—2020年沂河流域年均NDVI空間分布見圖4A。可以看出，沂河流域NDVI多年均值范圍在0.16～0.91之間，流域內植被覆蓋空間分異較大。植被NDVI高值區多集中在流域南部與東部，其中蒙陰縣、平邑縣與費縣三縣交界處的蒙山山區NDVI高值最為集中，部分地區NDVI值在0.9以上，該區位于沂蒙山區腹地，生態環境較好，同時建有國家森林公園，植被覆蓋度較高。植被NDVI低值區主要集中在流域北部以及流域內各市(縣)域城鎮建成區，其中流域北部沂山山區和魯山山區NDVI值除少數地區外整體較南部蒙山山區低，表明流域內沂山山區和魯山山區多年來植被覆蓋情況較蒙山山區差。此外，流域內各市(縣)域城鎮建成區NDVI值均較低，其中臨沂市蘭山區低植被覆蓋區面積最大，表明城鎮用地的增加導致以上區域地表植被覆蓋情況較差。

在空間變化趨勢上，基于一元線性回歸分析在像元尺度上計算2000—2020年沂河流域NDVI空間變化趨勢，參考史曉亮等[31]對淮河流域植被NDVI劃分的變化趨勢等級，將沂河流域NDVI變化趨勢劃分為3個等級，并統計各個等級的面積及所占比例(圖4B，表3)。可以看出，植被NDVI基本不變的區域約占研究區總面積的78.35%，表明流域大部分地區植被覆蓋保持相對穩定，不存在明顯的增減趨勢。植被NDVI呈增加趨勢的區域所占比例為15.05%，主要分布在流域北部，以沂源縣和蒙陰縣境內最為集中，分析認為近些年來對沂蒙山區進行的一系列生態保護與修復措施，使得該區域植被NDVI增加面積持續增加。植被NDVI出現下降的區域面積相對較小，約占研究區總面積的6.60%，主要分布在流域東部和南部，其中流域南部平邑縣至蘭山區境內出現植被NDVI呈連續片狀減少區域，該區域主要集中在祊河流域，結合圖2可以看出，該區域20 a間草地面積不斷減少，同時耕地和建設用地面積不斷增加，分析認為該區域人類活動的增加導致植被NDVI的減少。此外，研究還發現在流域各城鎮建成區及其周邊區域內NDVI值呈減少趨勢，表明城鎮化的發展導致NDVI下降，植被覆蓋水平降低。

3.3 沂河流域地表徑流變化特征

基于GIS軟件結合SCS模型，在AMCⅡ條件下生成沂河流域20年逐年地表徑流深圖層，并依此生成2000—2005年，2006—2010年，2011—2015年，2016—2020年4個時期地表徑流深變化圖(圖5)。結果表明：(1) 流域不同時期徑流深均呈現出不同程度的由西北向東南逐漸增加趨勢，徑流深低值主要集中在流域北部及南部蒙山地區，高值主要集中在流域東部及南部。(2) 不同時期，沂河流域徑流深出現較大幅度波動，徑流存在明顯的豐枯變化，徑流深高值區和低值區空間范圍和面積均發生了不同程度變化。

和繼軍等[32]提出，在臨界值下，坡面產流量隨坡度的增加而增加。隨著坡度的增大，水的勢能變大，在地表中的流速增大，縮短了入滲時間，徑流深減小，從而使地表流量增大，沂河流域北部魯山山區、沂山山區及南部蒙山山區地形坡度較大，坡面產流能力較強，坡面產流主要向東南部地形較平坦的平原地區匯集。同時結合沂河流域土地利用圖發現，20年間流域南部耕地及建設用地面積不斷擴大，而北部土地利用結構相對穩定，林地與草地面積比例較南部高。建設用地所造成的不透水面導致地表水滲透比例減小，而耕地由于農業活動頻繁，促使土壤壓實和結皮，入滲速率和土壤蓄水含量降低[33]，從而導致地表徑流深增加。因此結合研究區高程及土地利用結構，分析認為沂河流域地表徑流深空間分布趨勢主要是地形坡度和土地利用結構共同作用的結果。此外，沂河流域徑流變化與降水變化具有顯著相關性，流域內季風氣候顯著，降水年內分配不均且年際變化大，導致沂河流域徑流出現較大幅度波動[23]。

圖4 2000-2020年沂河流域植被NDVI空間分布及變化趨勢分布

表3 2000-2020年沂河流域NDVI變化趨勢統計

3.4 植被NDVI與地表徑流的相關性分析

基于SCS模型模擬沂河流域20年間逐年地表徑流深，逐像元計算20年間年最大NDVI和徑流深之間的相關性，并對相關性進行顯著性檢驗(圖6)。結果顯示，沂河流域植被NDVI和地表徑流之間正相關和負相關共存，統計表明，流域內88.79%的區域植被NDVI與徑流呈正相關，11.21%的區域植被NDVI與徑流呈負相關。由圖6A可知，流域北部呈正相關區域面積遠大于呈負相關區域，呈負相關的區域則主要集中在流域南部，尤其以蘭山區為最。此外，流域內植被NDVI與徑流變化相關系數有較大分異，流域正相關系數較大，部分地區可達0.8以上，而在流域內部分呈負相關的區域，相關系數則可接近-0.7。

圖5 2000-2020年沂河流域不同時間段地表徑流深

對相關性結果采用p<0.05置信水平進行顯著性檢驗(圖6B，表4)，可以看出研究區內以不顯著正相關為主，研究區內正、負相關達到顯著水平的面積分別占統計總面積的34.07%和4.72%，其中呈顯著正相關的區域主要集中在流域北部，對比圖4B，發現呈顯著正相關的區域與NDVI呈改善趨勢的區域基本一致，表明隨著該區域植被NDVI的增加，地表徑流也增加。呈顯著負相關的區域主要集中在研究區南部，尤其以臨沂市蘭山區最為集中，結合圖2和圖4A，分析認為城鎮建設用地的增加導致城市地表不透水面增加，同時造成植被覆蓋的減少，導致地表徑流的增加。此外，植被NDVI與地表徑流呈顯著正相關與顯著負相關的像元在流域內相互交錯分布，表明流域內植被NDVI與徑流變化關系在小空間尺度下存在明顯差異。

圖6 沂河流域植被NDVI與徑流變化相關系數及顯著性

表4 沂河流域植被NDVI與徑流變化相關系數顯著性統計

4 討論與結論

4.1 討 論

植被覆蓋變化受多種因素耦合作用，在空間上呈現差異性，由其引發的水文效應則更加復雜。由于區域氣候和地理的差異性，研究尺度與方法的不同，以及植被本身的復雜性，特定流域內植被變化對水文過程的影響結果并不一致。近幾十年以來，沂河流域經濟社會快速發展，流域內植被覆蓋情況的改變及其引發的水文效應勢必對現在及未來沂河流域生態環境產生影響。

本文基于遙感影像，在像元尺度分析植被變化及其對徑流的影響，實現在較高精度下分析沂河流域植被變化與徑流的相關關系。但需要指出的是本文在利用SCS模型模擬地表徑流時，只基于前期土壤水分處于平均狀態并結合三期土地利用模擬地表徑流，且沂河流域降雨主要集中在主汛期的7—9月份，多數降雨為集中性暴雨，因此基于SCS模型的沂河流域徑流模擬精度依然尚待改善。此外，本文只是定量描述了植被變化與徑流的關系，而受植被變化時空特征影響冠層截留、蒸散發等水文過程，需要進一步利用分布式水文模型，在小尺度背景下，從產匯流過程上具體分析植被空間格局變化對水文過程的影響機制。

沂河是魯南與蘇北地區重要的山洪河道，同時沂河流域也是魯南山地生態環境的重要組成部分，因此沂河流域植被覆蓋變化及其對徑流的影響對魯南與蘇北地區經濟社會發展意義重大。未來沂河流域經濟社會發展應該更加注重土地利用合理規劃，加強生態環境保護，而未來沂河流域生態環境變化研究也應該將目光放在自然因素與人為因素相互耦合作用研究上。

4.2 結 論

(1) 2000—2020年沂河流域植被NDVI呈波動增加趨勢，增長速率為1.46%/10 a。流域植被NDVI空間分異較大，高值區主要分布在流域南部與東部，低值區主要集中在流域北部以及流域內城鎮建成區。流域植被NDVI空間變化趨勢整體較為穩定，以基本不變為主，NDVI呈增加趨勢的區域主要集中在流域北部，占研究區面積的15.05%，NDVI呈減少趨勢的區域主要分布在流域東部和南部及各城鎮建成區內，占比為6.60%。

(2) 2000—2020年沂河流域不同時期徑流深均呈現出不同程度的由西北向東南逐漸增加趨勢，徑流深出現較大幅度波動，徑流存在明顯的豐枯變化。徑流深低值主要集中在流域北部及南部蒙山地區，高值主要集中在流域東部及南部，徑流深高值區和低值區空間范圍和面積均發生了不同程度變化。

(3) 沂河流域植被NDVI和地表徑流深的相關關系呈正相關與負相關共存的狀態，流域不同空間位置的植被變化與徑流變化的關系存在較大分異。88.79%的區域植被NDVI與徑流呈正相關，以不顯著正相關為主，呈顯著正相關的區域主要集中在流域北部，植被NDVI與徑流呈負相關僅占比11.21%，其中顯著負相關區域主要集中在流域南部蘭山區境內。