沂蒙山區(qū)土壤侵蝕時(shí)空變化及其驅(qū)動(dòng)因素分析

黃 鑫,陳 紅,徐小任,唐小倩,肖 新,陳世勇,梅大偉,王 梁① (.安徽科技學(xué)院資源與環(huán)境學(xué)院,安徽 滁州 3300;.山東省水土保持與環(huán)境保育重點(diǎn)實(shí)驗(yàn)室/臨沂大學(xué)資源環(huán)境學(xué)院,山東 臨沂 76005)

中國(guó)是世界上土壤侵蝕最嚴(yán)重的國(guó)家之一,2020年全國(guó)水蝕和風(fēng)蝕土壤面積為269.27萬(wàn)km2,約占國(guó)土面積的28.0%[1]。北方土石山區(qū)是我國(guó)水力侵蝕最嚴(yán)重的區(qū)域之一,土壤侵蝕已成為制約區(qū)域發(fā)展的首要問(wèn)題。沂蒙山區(qū)作為北方土石山區(qū)的典型代表,其特有的土壤地理環(huán)境、極不均勻且強(qiáng)度大的降水分布、復(fù)雜多樣的地質(zhì)構(gòu)造以及不合理的人類活動(dòng)成為土壤侵蝕的內(nèi)外驅(qū)動(dòng)力,淺薄的土層流失殆盡造成嚴(yán)重的水土流失問(wèn)題,進(jìn)而導(dǎo)致土地退化、自然災(zāi)害頻發(fā),生態(tài)系統(tǒng)瀕臨崩潰。因此,定量評(píng)價(jià)該區(qū)域的土壤侵蝕現(xiàn)狀,科學(xué)認(rèn)識(shí)和定量評(píng)估自然因素、人類活動(dòng)及其耦合作用對(duì)土壤侵蝕的影響,對(duì)遏制該地區(qū)的水土流失及推進(jìn)生態(tài)安全建設(shè)至關(guān)重要。

對(duì)于土壤侵蝕這一影響全球生態(tài)安全的重要問(wèn)題,多年來(lái)國(guó)內(nèi)外學(xué)者在土壤侵蝕機(jī)理、模型模擬和驅(qū)動(dòng)因素等方面進(jìn)行多領(lǐng)域研究,建立了具有不同特色的侵蝕預(yù)報(bào)模型,主要分為物理過(guò)程模型、分布式模型與經(jīng)驗(yàn)統(tǒng)計(jì)模型3類。典型的物理過(guò)程模型主要有水蝕預(yù)報(bào)模型(WEEP)[2]、歐洲土壤侵蝕模型(EUROSEM)[3]等,分布式模型以分布式水文模型(SHE)[4]最為經(jīng)典。但由于物理過(guò)程模型與分布式模型所需參數(shù)眾多,難以推廣應(yīng)用,因此以通用土壤流失方程(USLE)/修正版通用土壤流失方程(RUSLE)為代表的經(jīng)驗(yàn)統(tǒng)計(jì)模型應(yīng)用最為廣泛。RUSLE在大量實(shí)測(cè)數(shù)據(jù)基礎(chǔ)上對(duì)USLE的局限進(jìn)行修正,形式簡(jiǎn)單靈活,參數(shù)獲取方便,便于與地理信息系統(tǒng)(GIS)和遙感(RS)技術(shù)結(jié)合,是目前發(fā)展較為成熟的土壤侵蝕模型[5-6]。在土壤侵蝕驅(qū)動(dòng)因素研究方面,目前常用的主要方法,如線性與非線性擬合、主成分分析、偏最小二乘回歸、結(jié)構(gòu)方程等,對(duì)于主控因子識(shí)別、因子共線性處理以及剖析因子的直接或間接作用等方面起到了關(guān)鍵作用。但是由于驅(qū)動(dòng)因素的時(shí)空變異性及其之間的交互作用致使土壤侵蝕具有時(shí)空異質(zhì)性。地理探測(cè)器是王勁峰等[7]基于空間方差理論提出的一種分析地理現(xiàn)象空間異質(zhì)性的新統(tǒng)計(jì)學(xué)方法,不僅可以揭示地理現(xiàn)象背后的驅(qū)動(dòng)力,而且能夠量化因子間的交互作用、有效識(shí)別高風(fēng)險(xiǎn)區(qū),成為當(dāng)前運(yùn)用較為廣泛的前沿方法之一[8]。

土壤侵蝕在年際之間存在一定變化,目前關(guān)于沂蒙山區(qū)土壤侵蝕的研究主要集中于單一年份土壤侵蝕,對(duì)于長(zhǎng)時(shí)間序列土壤侵蝕與其驅(qū)動(dòng)因素的交互作用研究關(guān)注較少。深入這方面的研究有助于進(jìn)一步認(rèn)識(shí)土壤侵蝕的機(jī)理與變化規(guī)律,特別是自然因素與人類活動(dòng)及其交互作用下的土壤侵蝕定量歸因。因此,筆者以北方土石山區(qū)典型代表沂蒙山區(qū)為對(duì)象,基于GIS和遙感RS技術(shù),利用RUSLE探究沂蒙山區(qū)2000—2018年土壤侵蝕的時(shí)空變化特征,進(jìn)而運(yùn)用地理探測(cè)器模型量化不同驅(qū)動(dòng)因素之間的交互作用及其對(duì)土壤侵蝕的貢獻(xiàn)率,揭示不同時(shí)期影響土壤侵蝕的主導(dǎo)因素及其交互因素,結(jié)果可為沂蒙山區(qū)的水土保持規(guī)劃及國(guó)土空間規(guī)劃提供參考,也可以為該區(qū)域乃至北方土石山區(qū)生態(tài)環(huán)境建設(shè)提供科學(xué)依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

沂蒙山區(qū)位于山東省中南部(34°22～36°23′ N,116°34′～119°39′ E),范圍涉及臨沂市、淄博市、日照市、濟(jì)寧市、棗莊市等27個(gè)縣(市、區(qū))(圖1),總面積約3.41萬(wàn)km2。研究區(qū)的氣候類型屬于暖溫帶季風(fēng)型大陸性氣候,年平均氣溫12～14 ℃,年平均降水量700～900 mm,夏季(6—8月)降水占全年降水量的60%以上。土壤以棕壤和褐土為主,地勢(shì)西北高東南低,地貌類型復(fù)雜,以山地丘陵為主。自然植被以落葉闊葉林為主,但受長(zhǎng)期人類活動(dòng)影響,自然植被只見于海拔較高的山區(qū),以人工植被或次生林為主。

圖1 沂蒙山區(qū)行政區(qū)域及數(shù)字高程模型Fig.1 Administrative division and digital elevation model of Yimeng Mountainous Area

沂蒙山區(qū)是我國(guó)水土流失治理的重點(diǎn)區(qū)域,在全國(guó)水土保持區(qū)劃中被定為魯中南低山丘陵土壤保持功能區(qū),在全國(guó)水土流失重點(diǎn)防治區(qū)劃分中被列為國(guó)家級(jí)重點(diǎn)治理區(qū),在國(guó)家水土保持重點(diǎn)工程建設(shè)中被納入第4期(2008—2012年)和第5期(2013—2017年)規(guī)劃項(xiàng)目。截至2018年底,國(guó)家和地方在沂蒙山區(qū)開展了以修建梯田、淤地壩為代表的水土保持綜合治理工程和以退耕還林還草為代表的一系列生態(tài)修復(fù)工程,小流域水土流失治理度達(dá)70%以上,林草覆蓋度達(dá)30%以上,水土流失情況得到了一定改善。

1.2 數(shù)據(jù)來(lái)源

研究所使用的主要數(shù)據(jù)有:沂蒙山區(qū)矢量范圍,來(lái)源于全國(guó)地理信息資源目錄服務(wù)系統(tǒng)(http:∥www.webmap.cn);1∶100萬(wàn)中國(guó)土壤數(shù)據(jù)集,來(lái)源于中國(guó)科學(xué)院南京土壤研究所(http:∥vdb3.soil.csdb.cn/),包含土壤類型、質(zhì)地、有機(jī)質(zhì)等數(shù)據(jù);數(shù)字高程數(shù)據(jù)ASTER GDEM,來(lái)源于地理空間數(shù)據(jù)云平臺(tái)(http:∥www.gscloud.cn);沂蒙山區(qū)及周邊29個(gè)氣象站點(diǎn)2000—2018年的逐日降雨數(shù)據(jù),來(lái)源于中國(guó)氣象局氣象數(shù)據(jù)中心(http:∥data.cma.cn);土地利用數(shù)據(jù)、地貌形態(tài)數(shù)據(jù),來(lái)源于中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn);植被覆蓋指數(shù)NDVI數(shù)據(jù),采用NASA發(fā)布的MOD13Q1產(chǎn)品(http:∥lpdaac.usgs.gov)。

2 研究方法

2.1 RUSLE模型

采用RUSLE模型評(píng)估沂蒙山區(qū)土壤侵蝕,計(jì)算公式為

A=R×K×L×S×C×P。

(1)

式(1)中,A為年土壤侵蝕模數(shù),t·hm-2·a-1;R為降雨侵蝕力因子,MJ·mm·hm-2·h-1·a-1;K為土壤可蝕性因子,t·hm2·h·MJ-1·mm-1·hm-2;L、S為坡長(zhǎng)、坡度因子;C為植被覆蓋與管理因子;P為水土保持措施因子。

2.1.1降雨侵蝕力因子(R)

選用WISCHMEIER提出并由ARNOLDUS等[9]修正的基于月平均降雨量和年平均降雨量的經(jīng)驗(yàn)公式來(lái)估算降雨侵蝕力。

式(2)中,pi為月平均降雨量,mm;p為年平均降雨量,mm。

2.1.2土壤可蝕性因子(K)

選用WILLIAMS等[10]在侵蝕/生產(chǎn)力影響模型(EPIC)中土壤可蝕性因子的計(jì)算方法進(jìn)行估算。

(3)

(4)

式(3)～(4)中,SAN、SIL、CLA、Co分別為砂粒、粉粒、黏粒和有機(jī)質(zhì)的含量,%;SNI為計(jì)算參數(shù)。

2.1.3坡長(zhǎng)、坡度因子(L、S)

坡長(zhǎng)因子的計(jì)算選用WISCHMERIER等[11]提出的計(jì)算公式。當(dāng)坡度≤18%,選用MCCOOL等[12]提出的計(jì)算公式;當(dāng)坡度>18%,選用LIU等[13]改進(jìn)后的計(jì)算公式進(jìn)行估算。

(5)

(6)

式(5)～(6)中,λ為坡長(zhǎng),m;α為坡長(zhǎng)指數(shù);β為坡度修正值;θ為DEM提取的坡度值,(°);標(biāo)準(zhǔn)小區(qū)坡長(zhǎng)為22.13 m。

2.1.4植被覆蓋和管理因子(C)

選用蔡崇法等[14]建立的c值和植被覆蓋度之間的回歸方程計(jì)算C值,核心算法如下

c=(INDV-INDV,s)/(INDV,v-INDV,s),

(7)

(8)

式(7)～(8)中,INDV為歸一化植被指數(shù);INDV,s為像元在裸地覆蓋狀態(tài)下的NDVI值;INDV,v為在植被全覆蓋狀態(tài)下的NDVI值。

2.1.5水土保持措施因子(P)

水土保持措施因子P是指采用一定的水土保持措施后的土壤流失量與不采取任何措施的土壤流失量之比。P值范圍在0～1之間,0代表不發(fā)生土壤侵蝕,1代表未采取水土保持措施或措施完全失效。對(duì)于P值,目前通常采用結(jié)合實(shí)地考察資料給不同土地利用類型賦值的方法確定。結(jié)合文獻(xiàn)[15-17]和研究區(qū)土地利用情況,將水田賦值為0.15,旱地賦值為0.35,林地、疏林地、草地、未利用地賦值為1,水域、居民點(diǎn)、建設(shè)用地、裸巖賦值為0。

2.2 地理探測(cè)器

研究主要利用地理探測(cè)器中的因子探測(cè)器、交互探測(cè)器和風(fēng)險(xiǎn)探測(cè)器,定量分析沂蒙山區(qū)土壤侵蝕的驅(qū)動(dòng)因子以及因子間的相互影響,并識(shí)別土壤侵蝕高風(fēng)險(xiǎn)區(qū)域。

(9)

(10)

ST=Nσ2。

(11)

式(9)～(11)中,q為土壤侵蝕因子的影響強(qiáng)度;L為因變量Y或自變量X的分層;Nh和N分別為層h和全區(qū)單元數(shù);σh2和σ2分別為層h和全區(qū)的Y值方差;SW為層內(nèi)方差和;ST為全區(qū)總方差。q的范圍在0～1之間,q值越大說(shuō)明土壤侵蝕Y的空間分異性越明顯。如果分類是由因子X造成的,q值越大表示因子X對(duì)土壤侵蝕Y的解釋力越強(qiáng),反之則越弱。

土壤侵蝕的發(fā)生與發(fā)展受眾多因素的影響和控制,是自然因素和人類活動(dòng)綜合作用的結(jié)果。土壤侵蝕的驅(qū)動(dòng)因素可以分為內(nèi)在因素和外在因素2類,內(nèi)在因素在較長(zhǎng)時(shí)間尺度上具有高度的穩(wěn)定性,其空間異質(zhì)性顯著影響土壤侵蝕的空間分布,而外在因素在較小尺度的時(shí)空變化是驅(qū)動(dòng)土壤侵蝕過(guò)程的直接誘因。基于此,選取地形地貌、坡度和海拔等內(nèi)在因素,降雨、土地利用類型和植被覆蓋度等外在因素作為自變量X,2000、2005、2010、2015和2018年土壤侵蝕強(qiáng)度作為因變量Y代入地理探測(cè)器中進(jìn)行運(yùn)算。由于地理探測(cè)器的輸入變量要求為類別數(shù)據(jù),首先需對(duì)連續(xù)型變量做離散化處理,結(jié)合數(shù)據(jù)離散化方法及先驗(yàn)知識(shí)[9,18-19],地形地貌使用《中華人民共和國(guó)地貌圖集1∶100萬(wàn)》中的類別編號(hào),分為26類,土地利用類型分為耕地、林地、草地、水域、建設(shè)用地、未利用地6類,將植被覆蓋度分為8類(≤0.3、>0.3～0.4、>0.4～0.5、>0.5～0.6、>0.6～0.7、>0.7～0.8、>0.8～0.9、>0.9～1.0),海拔分為3類(≤500、>500～1 000、>1 000～1 500 m),年平均降雨量等間距分為9類,坡度分為6類(≤5°、>5°～10°、>10°～15°、>15°～20°、>20°～25°、>25°),并將土壤侵蝕強(qiáng)度及驅(qū)動(dòng)因子分類值賦予1 km×1 km格網(wǎng)點(diǎn)上,作為地理探測(cè)器的運(yùn)行數(shù)據(jù)。

3 結(jié)果與分析

3.1 土壤侵蝕強(qiáng)度的時(shí)間變化

沂蒙山區(qū)2000、2005、2010、2015與2018年的平均土壤侵蝕模數(shù)分別為5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1,2005年侵蝕最為嚴(yán)重,2015年侵蝕程度最低。將研究區(qū)5期土壤侵蝕強(qiáng)度按照SL 190—2007《土壤侵蝕強(qiáng)度分類分級(jí)標(biāo)準(zhǔn)》劃分為微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈、劇烈6個(gè)級(jí)別(圖2),各侵蝕級(jí)別面積占比見表1。由表1可知,2000—2018年各侵蝕級(jí)別整體變化不大,微度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕面積分別增加2.95%、0.81%、0.72%和0.71%,輕度和中度侵蝕面積分別下降4.18%和1.01%。

表1 2000—2018年沂蒙山區(qū)不同土壤侵蝕等級(jí)面積比例Table 1 Percentage change in the area of different soil erosion levels in the Yimeng Mountainous Area from 2000 to 2018

圖2 2000—2018年沂蒙山區(qū)土壤侵蝕強(qiáng)度分級(jí)Fig.2 Soil erosion intensity level in the Yimeng Mountainous Area from 2000 to 2018

為進(jìn)一步探討不同土壤侵蝕強(qiáng)度之間的時(shí)間變化特征,分別對(duì)5個(gè)時(shí)期內(nèi)不同侵蝕強(qiáng)度面積數(shù)據(jù)進(jìn)行統(tǒng)計(jì),制成土壤侵蝕強(qiáng)度轉(zhuǎn)移矩陣(表2)。

表2 2000—2018年沂蒙山區(qū)土壤侵蝕強(qiáng)度轉(zhuǎn)移矩陣Table 2 Soil erosion intensity transfer matrix for Yimeng Mountainous Area from 2000 to 2018

結(jié)果表明,2000—2005年微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕的穩(wěn)定率(2005與2000年土壤侵蝕強(qiáng)度等級(jí)相同的面積之比)分別為99.82%、0.38%、0.43%、0.00%、59.37%和0.00%。除微度侵蝕外,其他級(jí)別侵蝕強(qiáng)度均發(fā)生極大變化,整體趨勢(shì)為低等級(jí)向高等級(jí)轉(zhuǎn)移。2005—2010年微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕的穩(wěn)定率分別為67.46%、3.80%、9.14%、15.20%、11.39%和27.66%。除微度侵蝕外,其他級(jí)別侵蝕強(qiáng)度均發(fā)生較大變化,整體趨勢(shì)為高等級(jí)向低等級(jí)轉(zhuǎn)移。2010—2015年土壤侵蝕強(qiáng)度整體變化趨勢(shì)仍為高等級(jí)向低等級(jí)轉(zhuǎn)移,輕度、中度、強(qiáng)烈和極強(qiáng)烈侵蝕轉(zhuǎn)化為微度侵蝕的面積比例分別為77.14%、47.27%、45.45%和38.46%,劇烈侵蝕全部轉(zhuǎn)化為極強(qiáng)烈侵蝕。2015—2018年微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕的穩(wěn)定率分別為70.17%、18.55%、10.37%、12.15%、10.18%和30.94%,整體趨勢(shì)依舊表現(xiàn)為高等級(jí)向低等級(jí)轉(zhuǎn)移,但也存在低等級(jí)向高等級(jí)轉(zhuǎn)移的現(xiàn)象:輕度侵蝕中27.62%轉(zhuǎn)化為極強(qiáng)烈侵蝕,中度侵蝕中15.49%轉(zhuǎn)化為極強(qiáng)烈侵蝕,強(qiáng)烈侵蝕中12.55%轉(zhuǎn)化為極強(qiáng)烈侵蝕,極強(qiáng)烈侵蝕中21.25%轉(zhuǎn)化為劇烈侵蝕。2000—2018年微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕的穩(wěn)定率分別為98.74%、28.27%、12.50%、0.00%、21.43%和50.00%,輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕向低等級(jí)轉(zhuǎn)移的面積占比分別為67.70%、60.00%、80.00%、78.57%和50.00%。可以看出,2000—2018年土壤侵蝕強(qiáng)度向低等級(jí)轉(zhuǎn)化趨勢(shì)明顯,說(shuō)明研究區(qū)土壤侵蝕整體上呈現(xiàn)改善趨勢(shì)。

3.2 土壤侵蝕強(qiáng)度的空間變化

疊加分析沂蒙山區(qū)不同時(shí)期土壤侵蝕強(qiáng)度分布圖與沂蒙山區(qū)行政區(qū)劃圖,獲得沂蒙山區(qū)2000—2018年各縣區(qū)土壤侵蝕強(qiáng)度的面積比例(圖3),提取各縣域平均土壤侵蝕模數(shù),根據(jù)SL 190—2007劃分土壤侵蝕等級(jí)(圖4)。微度、輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕在不同年份各縣區(qū)均有分布,其中微度和輕度侵蝕分布面積廣泛,其他4個(gè)等級(jí)在不同年份的分布特征差異明顯。2000年各縣區(qū)強(qiáng)烈侵蝕、極強(qiáng)烈侵蝕和劇烈侵蝕面積占比分別為1.05%～8.33%、1.12%～6.75%和1.05%～7.49%,其中山亭區(qū)占比最大,微山縣最小;中度侵蝕在沂源縣和山亭區(qū)面積占比較大,分別為18.36%和21.41%。2005年各縣區(qū)強(qiáng)烈侵蝕和極強(qiáng)烈侵蝕面積占比分別為2.18%～10.86%和1.58%～10.23%;中度侵蝕在沂源縣和山亭區(qū)的面積占比最大,分別為15.80%和15.40%。2010年各縣區(qū)強(qiáng)烈侵蝕和極強(qiáng)烈侵蝕面積占比分別為2.33%～14.81%和1.23%～11.45%;中度侵蝕在蒙陰縣、費(fèi)縣、沂源縣和嵐山區(qū)的面積占比較大,分別為16.02%、14.59%、14.02%和14.29%。2005和2010年劇烈侵蝕面積占比在山亭區(qū)(13.08%和13.97%)和鄒城市(10.28%和11.00%)最大。2015年各縣區(qū)中度侵蝕和強(qiáng)烈侵蝕面積占比分別為1.81%～13.45%和0.78%～6.06%;極強(qiáng)烈侵蝕在莒南縣面積占比最大(4.86%);劇烈侵蝕除莒縣(5.72%)外其他26個(gè)縣區(qū)面積占比均較小。2018年各縣區(qū)極強(qiáng)烈侵蝕和劇烈侵蝕的面積占比分別為0.76%～10.02%和1.88%～11.92%;強(qiáng)烈侵蝕在沂源縣(10.69%)和山亭區(qū)(10.55%)的面積占比較大;中度侵蝕在蒙陰縣、平邑縣、沂源縣和山亭區(qū)的面積占比較大,為14.98%～15.34%。

圖3 2000—2018年沂蒙山區(qū)各縣區(qū)土壤侵蝕強(qiáng)度的面積比例Fig.3 Percentage change in the area of different soil erosion intensity in each county in the Yimeng Mountainous Area from 2000 to 2018

圖4 2000—2018年沂蒙山區(qū)各縣區(qū)土壤侵蝕分級(jí)Fig.4 Soil erosion classification in each country in the Yimeng Mountainous Area from 2000 to 2018

綜上可知,2000—2018年強(qiáng)烈及以上等級(jí)的侵蝕主要分布在西部和北部,集中于山亭區(qū)和沂南縣一帶,2個(gè)縣區(qū)的平均土壤侵蝕模數(shù)分別為9 217.55 和8 005.60 t·km-2·a-1,屬于極強(qiáng)烈侵蝕,強(qiáng)烈以下等級(jí)的侵蝕主要分布在西南部和中東部。隨時(shí)間變化,土壤侵蝕呈2000—2005年向北和西北方向擴(kuò)散、2005—2010年向西北方向擴(kuò)散、2010—2015年向中部擴(kuò)散、2015—2018年向西和東北方向擴(kuò)散的趨勢(shì)。

3.3 土壤侵蝕驅(qū)動(dòng)因素變化分析

運(yùn)用地理探測(cè)器的因子探測(cè)器獲取2000、2005、2010、2015和2018年地理探測(cè)器的q值(表3)。2000年土壤侵蝕驅(qū)動(dòng)因子的貢獻(xiàn)力由大到小依次為q(X5)、q(X1)、q(X2)、q(X6)、q(X4)、q(X3),其中土地利用類型起到主導(dǎo)作用,其次具有較強(qiáng)解釋力的是地形地貌。2005年土壤侵蝕驅(qū)動(dòng)因子的貢獻(xiàn)力由大到小依次為q(X5)、q(X1)、q(X4)、q(X2)、q(X6)、q(X3),其中起主導(dǎo)作用的仍然是土地利用類型,地形地貌、植被覆蓋度都具有較強(qiáng)的解釋力。2010年土壤侵蝕驅(qū)動(dòng)因子的貢獻(xiàn)力由大到小依次為q(X5)、q(X1)、q(X6)、q(X4)、q(X2)、q(X3),其中起主導(dǎo)作用的因子與上一時(shí)期相同,其他有較強(qiáng)解釋力的因子為地形地貌、年平均降雨量與植被覆蓋度。2015年土壤侵蝕驅(qū)動(dòng)因子的貢獻(xiàn)力由大到小依次為q(X5)、q(X1)、q(X4)、q(X2)、q(X6)、q(X3),該時(shí)期主導(dǎo)因子是土地利用類型,另外具有較高解釋力的是地形地貌。2018年土壤侵蝕驅(qū)動(dòng)因子的貢獻(xiàn)力由大到小依次為q(X5)、q(X1)、q(X6)、q(X2)、q(X4)、q(X3),該時(shí)期主導(dǎo)因子依然是土地利用類型,地形地貌、年平均降雨量都具有較強(qiáng)的解釋力,這一時(shí)期驅(qū)動(dòng)因子的解釋力較上一時(shí)期均有明顯提高,土壤侵蝕強(qiáng)度的變化也最劇烈。

運(yùn)用地理探測(cè)器的交互探測(cè)器獲取2000、2005、2010、2015和2018年各個(gè)驅(qū)動(dòng)因子對(duì)土壤侵蝕的交互作用(表4～5)。2000年土壤侵蝕的交互作用均為雙因子增強(qiáng)與非線性增強(qiáng),沒有單獨(dú)起作用的因子,土地利用類型、地形地貌可以增強(qiáng)各個(gè)驅(qū)動(dòng)因子對(duì)土壤侵蝕的解釋。2005年各個(gè)因子對(duì)土壤侵蝕的交互作用大多為非線性增強(qiáng),雙因子增強(qiáng)較少,在這一時(shí)期解釋力較強(qiáng)的交互因子與2000年趨同。2010年各個(gè)因子對(duì)土壤侵蝕的交互作用大多為雙因子增強(qiáng),非線性增強(qiáng)較少,其中與土地利用類型交互的因子解釋力增強(qiáng),土地利用類型、地形地貌、年平均降雨量和植被覆蓋度在這一時(shí)期對(duì)土壤侵蝕有顯著影響。2015年各驅(qū)動(dòng)因子對(duì)土壤侵蝕的交互作用與前3個(gè)時(shí)期相同,仍然還是非線性增強(qiáng)和雙因子增強(qiáng),土地利用類型和地形地貌在這一時(shí)期還是對(duì)土壤侵蝕有顯著影響。2018年各因子對(duì)土壤侵蝕的交互作用仍為非線性增強(qiáng)和雙因子增強(qiáng),在這一時(shí)期土地利用類型、地形地貌、年平均降雨量可以增強(qiáng)各驅(qū)動(dòng)因子的解釋力,坡度對(duì)土壤侵蝕有一定的引導(dǎo)作用。

表4 2000—2018年驅(qū)動(dòng)因子交互探測(cè)結(jié)果Table 4 Results of the interaction of driving factors from 2000 to 2018

表5 2000—2018年沂蒙山區(qū)因子交互貢獻(xiàn)率Table 5 Interaction contribution of factors in the Yimeng Mountainous Area from 2000 to 2018

運(yùn)用地理探測(cè)器的風(fēng)險(xiǎn)探測(cè)器識(shí)別2000、2005、2010、2015和2018年土壤侵蝕高風(fēng)險(xiǎn)區(qū),結(jié)果顯示,不同時(shí)期土壤侵蝕高風(fēng)險(xiǎn)區(qū)分布較為一致。地形地貌為中起伏的區(qū)域土壤侵蝕風(fēng)險(xiǎn)較高,而海拔分類中易發(fā)生土壤侵蝕的級(jí)別為500～1 000 m。疊加這2個(gè)因子分布圖發(fā)現(xiàn),2個(gè)高風(fēng)險(xiǎn)區(qū)有較大重疊,包括沂源縣以及蒙山地區(qū)(蒙陰、平邑、費(fèi)縣交界處),這部分區(qū)域地表起伏大,陡坡較多,侵蝕較為嚴(yán)重。通過(guò)對(duì)比不同坡度分區(qū)對(duì)應(yīng)的土壤侵蝕量可以發(fā)現(xiàn),隨著坡度的增加,土壤侵蝕量呈先增加后減少的趨勢(shì),在15°～ 20°區(qū)間出現(xiàn)拐點(diǎn),即為坡度的高風(fēng)險(xiǎn)區(qū)域。植被覆蓋度處于0.4～0.6之間的區(qū)域土壤侵蝕風(fēng)險(xiǎn)較高。比較發(fā)現(xiàn),植被覆蓋度介于0.4～0.6之間時(shí)為臨界值,低于這個(gè)臨界值土壤侵蝕隨植被覆蓋度的增加而增加,高于這個(gè)臨界值土壤侵蝕隨植被覆蓋度的增加而降低,說(shuō)明植被覆蓋度對(duì)土壤侵蝕有一定的抑制作用。各土地利用類型中,草地發(fā)生土壤侵蝕的風(fēng)險(xiǎn)最高,這是因?yàn)檠芯繀^(qū)草地主要分布在坡度較大的山區(qū),環(huán)境條件較差,同時(shí)科學(xué)有效的水土管理措施實(shí)施并不到位,造成土壤侵蝕現(xiàn)象較為嚴(yán)峻。疊加因子分布圖發(fā)現(xiàn),年平均降雨量的高風(fēng)險(xiǎn)區(qū)在空間上也對(duì)應(yīng)了土地利用類型的高風(fēng)險(xiǎn)區(qū),主要分布在沂源縣、沂水縣、蒙陰縣、平邑縣以及山亭區(qū)等地,這部分區(qū)域山多坡陡,降雨強(qiáng)度大,加之植被覆蓋度低,造成侵蝕作用強(qiáng)烈,水土流失嚴(yán)重。

土地利用類型是該研究區(qū)土壤侵蝕的主導(dǎo)因子,為進(jìn)一步分析不同土地利用類型的土壤侵蝕狀況,對(duì)沂蒙山區(qū)不同土地利用類型的土壤侵蝕強(qiáng)度平均值和面積占比進(jìn)行統(tǒng)計(jì)(表6)。土地利用類型侵蝕強(qiáng)度排序?yàn)椴莸?、林地、耕地、水域、建設(shè)用地和未利用地,不同土地利用類型侵蝕強(qiáng)度差異較大,2005、2010和2018年草地的平均侵蝕強(qiáng)度達(dá)到極強(qiáng)烈侵蝕,2010、2015與2018年未利用地與2015年建設(shè)用地的平均土壤侵蝕強(qiáng)度處于輕度侵蝕。

表6 2000—2018年不同土地利用類型土壤侵蝕平均模數(shù)及面積占比Table 6 Average modulus of soil erosion and area share of different land use types in the Yimeng Mountainous Area from 2000 to 2018

在不同土地利用類型中,侵蝕面積最大的為耕地,占侵蝕總面積的60%,其次是草地、建設(shè)用地和林地,未利用地的侵蝕面積較小,侵蝕面積占比不足1%。結(jié)合地理探測(cè)器分析結(jié)果和研究區(qū)土地利用類型面積占比和分布情況,草地是研究區(qū)亟需治理的土地利用類型,這與鐘旭珍等[20]、王猛等[21]進(jìn)行的土壤侵蝕定量歸因研究相符合,耕地不僅侵蝕面積占比最大,而且侵蝕強(qiáng)度突出,應(yīng)成為研究區(qū)應(yīng)加強(qiáng)預(yù)防和治理的重點(diǎn)區(qū)域,而對(duì)于林地、水域、建設(shè)用地和未利用地則需要重點(diǎn)關(guān)注并預(yù)防進(jìn)一步擴(kuò)展。

4 討論

4.1 土壤侵蝕時(shí)空變化

RUSLE模型模擬結(jié)果表明,沂蒙山區(qū)2000、2005、2010、2015與2018年的平均土壤侵蝕模數(shù)分別為5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1,2005年侵蝕最為嚴(yán)重,2015年侵蝕程度最低。2005年土壤侵蝕強(qiáng)度最高,主要原因在于2005年是研究區(qū)工業(yè)化、城鎮(zhèn)化進(jìn)程發(fā)展最快的時(shí)期,農(nóng)業(yè)及城鄉(xiāng)建設(shè)等人類活動(dòng)不斷擴(kuò)大,導(dǎo)致人為水土流失加重。2015年土壤侵蝕強(qiáng)度最低,主要是由于2015年沂蒙山區(qū)降雨出現(xiàn)異常情況,侵蝕性降雨量大幅減少且空間分布較為均勻。2015年平均降雨侵蝕力為1 323.93 MJ·mm·hm-2·h-1·a-1,較2000、2005、2010和2018年平均降雨侵蝕力(1 805.02、2 474.22、2 519.69和2 466.41 MJ·mm·hm-2·h-1·a-1)分別下降26.7%、46.5%、47.5%和46.3%。采用模擬降雨法[22]計(jì)算2000、2005、2010與2018年在2015年降雨下的土壤侵蝕情況,得出沂蒙山區(qū)2000、2005、2010、2015與2018年的平均土壤侵蝕模數(shù)為3 435.29、3 837.27、3 541.73、3 493.23和3 434.30 t·km-2·a-1,計(jì)算2000、2005、2010與2015年在2018年降雨下的土壤侵蝕情況,得出沂蒙山區(qū)2000、2005、2010、2015與2018年的平均土壤侵蝕模數(shù)為6 306.29、6 344.00、6 189.78、6 115.75和5 669.03 t·km-2·a-1,可以看出降雨之間的年際差別是2015年土壤侵蝕強(qiáng)度大幅度下降的原因。近20 a來(lái)強(qiáng)烈及以上等級(jí)的侵蝕主要分布在西部和北部,集中在山亭區(qū)、沂南縣一帶。這除了與不可避免的人類活動(dòng)不斷擴(kuò)大有關(guān)外,這部分地區(qū)長(zhǎng)期以來(lái)開采礦產(chǎn)資源,形成大量廢棄礦山,也在一定程度上加劇了該區(qū)域的土壤侵蝕。

4.2 驅(qū)動(dòng)因素時(shí)空變化

不同時(shí)期各驅(qū)動(dòng)因子對(duì)土壤侵蝕的影響程度顯著不同。土地利用類型在近20 a對(duì)土壤侵蝕始終具有最高的貢獻(xiàn)率,表明土地利用類型是影響沂蒙山區(qū)土壤侵蝕的最主要因素。研究區(qū)地貌復(fù)雜,以山地、丘陵為主,這樣的地形地貌從宏觀上控制地表過(guò)程的發(fā)生和發(fā)展,對(duì)土壤侵蝕也具有主導(dǎo)性作用。雖然坡度在近20 a對(duì)土壤侵蝕的貢獻(xiàn)率不足20%,但坡度與其他因子交互作用可大大增強(qiáng)對(duì)土壤侵蝕的控制作用,由此也表明坡度是影響土壤侵蝕的主要因素之一。植被覆蓋度受人類活動(dòng)的影響強(qiáng)烈,1995—2000年間通過(guò)實(shí)施退耕還林政策,整體的植被覆蓋度得到了一定程度的增加,導(dǎo)致了植被覆蓋度空間分布趨于均勻的變化。受此影響,植被覆蓋度對(duì)土壤侵蝕的影響程度在2000年減弱。年均降雨量對(duì)土壤侵蝕的影響在2000、2005和2015年較小,這可能與這3個(gè)時(shí)期降雨量的空間分布較為均勻有關(guān)。隨著2010和2018年降雨量分布空間差異的增大,年均降雨量成為影響土壤侵蝕的重要因素之一。海拔對(duì)土壤侵蝕的影響程度在5個(gè)時(shí)期均未達(dá)到顯著水平,其原因可能在于海拔差異小且空間離散化程度較低。

由于影響土壤侵蝕的因素較多,同時(shí)受數(shù)據(jù)收集的限制,研究仍存在一些不足:RUSLE模型因子計(jì)算過(guò)程中各個(gè)因子的數(shù)據(jù)精度不同,在柵格計(jì)算過(guò)程中會(huì)產(chǎn)生一定的誤差。在地理探測(cè)器的應(yīng)用中主要選擇了部分典型的驅(qū)動(dòng)因子,今后的研究中若加入地層巖性、人口等更多因子進(jìn)行探討,將更有利于沂蒙山區(qū)土壤侵蝕重點(diǎn)治理區(qū)的判別及因地制宜的進(jìn)行治理。

5 結(jié)論

(1)沂蒙山區(qū)2000、2005、2010、2015與2018年的平均土壤侵蝕模數(shù)分別為5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1。2000—2018年輕度、中度、強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕向低等級(jí)轉(zhuǎn)移的面積占比分別為67.70%、60.00%、80.00%、78.57%和50.00%,土壤侵蝕強(qiáng)度向低等級(jí)轉(zhuǎn)化趨勢(shì)明顯。從研究區(qū)的總體特征來(lái)看,土壤侵蝕先升高后降低,土壤侵蝕呈改善趨勢(shì)。

(2)近20 a來(lái)沂蒙山區(qū)土壤侵蝕較為嚴(yán)重的區(qū)域主要分布在西部和北部,主要發(fā)生強(qiáng)烈及以上等級(jí)的侵蝕(≥5 000 t·km-2·a-1),強(qiáng)烈以下等級(jí)的侵蝕(<5 000 t·km-2·a-1)主要分布在西南部和中東部。隨時(shí)間變化,土壤侵蝕呈現(xiàn)2000—2005年向北和西北方向擴(kuò)散,2005—2010年向西北方向擴(kuò)散,2010—2015年向中部擴(kuò)散,2015—2018年向西和東北方向擴(kuò)散的趨勢(shì)。

(3)驅(qū)動(dòng)因子在不同時(shí)期的解釋能力不同,土地利用類型是影響沂蒙山區(qū)土壤侵蝕的主導(dǎo)因子,對(duì)土壤侵蝕的解釋力在42.52%～66.84%之間。不同時(shí)期驅(qū)動(dòng)因子的交互作用均會(huì)增強(qiáng)其對(duì)土壤侵蝕的解釋能力,土地利用類型和地形地貌與其他因子交互后對(duì)侵蝕強(qiáng)度的貢獻(xiàn)率明顯提高。土地利用類型為草地、地形地貌為中起伏中山、坡度位于15°～20°區(qū)間以及植被覆蓋度處于0.4～0.6區(qū)間的地區(qū)是水土防護(hù)治理的重點(diǎn)區(qū)域。