尹作堂, 常 軍
(山東師范大學(xué) 地理與環(huán)境學(xué)院, 山東 濟(jì)南 250358)
土壤侵蝕是指土壤及土壤母質(zhì)等成分在水、風(fēng)力和凍融等外力作用下被破壞、運(yùn)輸和沉積的過程,這是近年來糧食安全和生態(tài)安全的主要威脅[1,2]。因此,探索土壤侵蝕和驅(qū)動(dòng)因素的時(shí)空特征,并采取有效措施保護(hù)土壤資源和修復(fù)生態(tài)環(huán)境至關(guān)重要[3,4]。當(dāng)前,越來越多的研究開始關(guān)注土壤侵蝕的空間異質(zhì)性和尺度效應(yīng),如考慮行政邊界變化對土壤侵蝕影響因素造成的不確定性,以便將研究結(jié)果應(yīng)用于水土保持規(guī)劃與治理[1,3,5,6]。如Guo等[3]分析了城市尺度、縣級(jí)尺度與鄉(xiāng)鎮(zhèn)尺度上的京津冀地區(qū)土壤侵蝕空間分布特征與主要控制因素的變化情況;Liu等[6]分析了縣級(jí)尺度下,珠江三角洲水土保持與城市化水平、平均高程和歸一化植被指數(shù)(NDVI)地理加權(quán)回歸系數(shù)的空間分布。這些研究提出了更有針對性的水土保持治理與土地利用政策。
基于修正通用土壤流失方程(RUSLE)[7]進(jìn)行土壤侵蝕分析是當(dāng)前常用的研究方法之一。在中國,RUSLE模型已在青藏高原[8,9]、黃土高原[10-13]、喀斯特地貌[14,15]和其他地區(qū)[16-18]得到了廣泛應(yīng)用,并取得了良好的評(píng)估效果。地理探測器是探測地理現(xiàn)象空間分異性,揭示地理現(xiàn)象與因變量關(guān)系的一種空間統(tǒng)計(jì)學(xué)方法[19],因其能夠量化不同驅(qū)動(dòng)因素及其交互作用對土壤侵蝕的實(shí)質(zhì)性貢獻(xiàn),在最近幾年被廣泛應(yīng)用于土壤侵蝕驅(qū)動(dòng)因素的研究[3-5,15,20]。但相關(guān)學(xué)者在使用地理探測器進(jìn)行土壤侵蝕驅(qū)動(dòng)因素分析時(shí),主要依靠王勁峰等[19]提出的方法及先驗(yàn)經(jīng)驗(yàn)進(jìn)行人為設(shè)定,忽視了數(shù)據(jù)離散化方法與分類數(shù)可能會(huì)對探測結(jié)果造成的影響[21,22],故本研究選用參數(shù)最優(yōu)地理探測器(optimal parameters-based geographical detector, OPGD)[22]對數(shù)據(jù)離散方法與間斷數(shù)量進(jìn)行最優(yōu)化設(shè)置,以提高地理探測精度。
作為土壤侵蝕最嚴(yán)重的區(qū)域之一,1990—2015年黃河年均土壤流失量2.0×109t[23],造成土地資源退化、土壤養(yǎng)分流失、河道淤積,給黃河流域的生態(tài)保護(hù)和高質(zhì)量發(fā)展造成巨大阻力。近年來,大量學(xué)者對黃河流域及相關(guān)地區(qū)開展了土壤侵蝕及其驅(qū)動(dòng)因素分析[4,23-25],但關(guān)于黃河流域土壤侵蝕驅(qū)動(dòng)因素的空間分布還有待深入研究。鑒于此,本文基于RUSLE模型與熱點(diǎn)分析,探討2000年、2010年和2020年黃河流域及縣域尺度下的土壤侵蝕狀況,并基于參數(shù)最優(yōu)地理探測器評(píng)估黃河流域與縣域尺度下的土壤侵蝕驅(qū)動(dòng)因素的解釋力,分析縣域尺度下土壤侵蝕驅(qū)動(dòng)因素的時(shí)空分布,以期為黃河流域水土保持及高質(zhì)量發(fā)展提供有效建議。
黃河流域地處95°02′ E ~ 119°43′ E,31°28′N~ 41°33′ N之間,整體地勢西高東低,流域內(nèi)海拔高差達(dá)6 200 m,地形起伏差異明顯,自西向東形成橫跨青藏高原、河套平原、鄂爾多斯高原、黃土高原和黃淮海平原的三級(jí)階梯(圖1)。區(qū)域內(nèi)氣候差異明顯,季節(jié)差異大,且降水時(shí)空分布不均,夏季降水量約占總降水量的70%,氣溫日較差和年較差較大;土壤包括高山土、干旱土、半淋溶土等多種類型,自然植被包括高寒草甸、草原、落葉林等。流域內(nèi)地理特征復(fù)雜,土壤侵蝕以水力侵蝕營力為主,同時(shí)存在風(fēng)力、凍融和重力等侵蝕營力[26]。近年來,黃河流域因水資源短缺、生態(tài)脆弱和水沙關(guān)系不協(xié)調(diào)等問題受到廣泛關(guān)注。
圖1 研究區(qū)概況圖Fig.1 Overview map of the study area 注:此圖基于國家自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站審圖號(hào) 為GS2020(4619)的標(biāo)準(zhǔn)地圖制作,底圖無修改。
1.2.1修正通用土壤流失方程
修正通用土壤流失方程(RUSLE)為:
A=R·K·LS·C·P
(1)
式中:A為年土壤侵蝕模數(shù)(t/(hm2·a));R為降雨侵蝕力((MJ·mm)/(hm2·h·a));K為土壤可蝕性因子((t·h)/(MJ·mm));LS為坡度因子S與坡長因子L(無量綱);C為地表植被覆蓋管理因子(無量綱);P為水土保持措施因子(無量綱)。
其中,基于章文波等[27]提出的日雨量估算侵蝕力方法,計(jì)算流域各氣象站點(diǎn)年降雨侵蝕力,通過Kriging空間插值處理,得到2000年、2010年和2020年黃河流域降雨侵蝕力因子;基于泛第三極(20國)土壤可蝕性因子(K)數(shù)據(jù)集[28],并在陳同德[26]、李天宏[29]、張科利[30,31]等人研究基礎(chǔ)上做出修正,K因子取值范圍為0.017 2~0.086 9 (t·h)/(MJ·mm);基于泛第三極(20國)坡度坡長因子(LS)數(shù)據(jù)集[32]提取黃河流域LS因子數(shù)據(jù)集,LS因子取值范圍為0~75.59;C因子以蔡崇法等[33]提出的基于植被覆蓋度的C因子估算方法進(jìn)行估算;基于RS與GIS提取法,參考文獻(xiàn)[4]、[23]、[34]并結(jié)合研究區(qū)實(shí)際情況,根據(jù)土地利用類型(LUCC)數(shù)據(jù)與坡度數(shù)據(jù)對P因子進(jìn)行賦值(表1)。
表1 不同土地利用類型的P值Tab.1 P values of different land use types
1.2.2熱點(diǎn)分析
熱點(diǎn)分析是一種空間自相關(guān)方法,可用于識(shí)別土壤侵蝕在空間上的高值(熱點(diǎn))與低值(冷點(diǎn))的聚類情況[35]。基于ArcGIS Pro熱點(diǎn)分析(Getis-OrdGi*)工具,統(tǒng)計(jì)識(shí)別具有統(tǒng)計(jì)顯著性的土壤侵蝕熱點(diǎn)與冷點(diǎn),計(jì)算公式為:
(2)
(3)
式中:Gi*為空間關(guān)聯(lián)指數(shù);E(Gi*)為Gi*(d)的期望值;Var(Gi*)為變異系數(shù);Wij(d)為空間權(quán)重矩陣;Xj為第j級(jí)土壤侵蝕強(qiáng)度等級(jí)的面積;n為土壤侵蝕強(qiáng)度等級(jí)數(shù);Z(Gi*)為Gi*的標(biāo)準(zhǔn)化Z值。標(biāo)準(zhǔn)化Z值為正且值大,表示土壤侵蝕高值的空間聚類(熱點(diǎn)),Z值為負(fù)且值小,則表示土壤侵蝕低值的空間聚類(冷點(diǎn))。
1.2.3參數(shù)最優(yōu)地理探測器(OPGD)
地理探測器是探測驅(qū)動(dòng)因素的一種統(tǒng)計(jì)學(xué)方法,包括分異及因子探測、交互作用探測、風(fēng)險(xiǎn)區(qū)探測與生態(tài)探測[19]。本研究使用參數(shù)最優(yōu)地理探測器R語言程序包GD[22]來分析單個(gè)因子與因子交互作用對土壤侵蝕的解釋力,其中探測結(jié)果q值表征解釋力的大小。具體設(shè)置如下:以土壤侵蝕模數(shù)作為因變量,以LUCC(X1)、地貌類型(X2)、土壤類型(X3)、年平均降雨量(X4)、坡度(X5)、海拔(X6)、植被覆蓋度(X7)、人口密度(X8)為驅(qū)動(dòng)因素,選用五種連續(xù)變量離散化算法(自然斷點(diǎn)法、標(biāo)準(zhǔn)差法、分位數(shù)法、等間隔法和幾何間隔法),間斷數(shù)量設(shè)置為3~8類,進(jìn)行土壤侵蝕驅(qū)動(dòng)因素分析。
主要數(shù)據(jù)來源:①黃河流域矢量面數(shù)據(jù)、地貌類型數(shù)據(jù)(1 km×1 km)、土壤類型數(shù)據(jù)(1 km×1 km)和土地利用類型數(shù)據(jù)(1 km×1 km)來源于中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心(https://www.resdc.cn/);②ASTER GDEM 數(shù)據(jù)(30 m×30 m)來源于中國科學(xué)院地理空間數(shù)據(jù)云平臺(tái)(http://www.gscloud.cn/);③日降雨資料(站點(diǎn)數(shù):113)來源于中國氣象數(shù)據(jù)網(wǎng)(https://data.cma.cn/);④NDVI數(shù)據(jù)(250 m×250 m)是在美國國家航空航天局提供的MOD13Q1數(shù)據(jù)產(chǎn)品基礎(chǔ)上,采用最大值合成法生成;⑤人口密度數(shù)據(jù)(1 km×1 km)來源于WorldPop project官方網(wǎng)站(https://www.worldpop.org/)。其中柵格數(shù)據(jù)在經(jīng)過精度分析、糾偏等處理后,重采樣為1 km分辨率。
依據(jù)水利部頒布的《土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)》(SL 190—2007)[36],將土壤侵蝕強(qiáng)度分為微度(0 ~1 000 t/(hm2·a))、輕度(1 000 ~ 2 500 t/(hm2·a))、中度(2 500 ~ 5 000 t/(hm2·a))、強(qiáng)烈(5 000 ~ 8 000 t/(hm2·a))、極強(qiáng)烈(8 000 ~ 15 000 t/(hm2·a))和劇烈(> 15 000 t/(hm2·a))六個(gè)級(jí)別。圖2中(a)~(c)分別2000年、2010年和2020年黃河流域以柵格像元為單元?jiǎng)澐值耐寥狼治g強(qiáng)度圖,圖2中(d)~(f)分別為2000年、2010年和2020年黃河流域以縣域?yàn)閱卧獎(jiǎng)澐值耐寥狼治g強(qiáng)度圖。
圖2 土壤侵蝕強(qiáng)度圖Fig.2 Map of soil erosion intensity 注:此圖基于國家自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站審圖號(hào)為GS2020(4619)的標(biāo)準(zhǔn)地圖制作,底圖無修改。
2000—2020年,黃河流域土壤侵蝕狀況總體呈現(xiàn)好轉(zhuǎn)趨勢,土壤侵蝕強(qiáng)度減弱區(qū)域約1.9×105km2,占流域總面積的24%,土壤侵蝕強(qiáng)度增強(qiáng)區(qū)域約6.5×104km2,占流域總面積的8%。2000年、2010年和2020年,黃河流域平均土壤侵蝕模數(shù)分別為1 994、1 860和1 384 t/(km2·a),水土流失面積分別為2.80×105、2.62×105和2.14×105km2。侵蝕強(qiáng)度較高的區(qū)域集中于黃土高原地區(qū),并且呈現(xiàn)東北—西南的條帶狀走向。但在此期間,河套平原北部的陰山山脈地區(qū),土壤侵蝕強(qiáng)度具有顯著的增強(qiáng)趨勢。
由圖2(d)~(f)可知,縣域尺度下不同區(qū)縣土壤侵蝕的空間差異性明顯。2000年,土壤侵蝕強(qiáng)度呈現(xiàn)極強(qiáng)烈的縣分別為陜西省子洲縣、子長縣、延川縣、清澗縣、延長縣,山西省石樓縣和甘肅省會(huì)寧縣,各縣平均土壤侵蝕模數(shù)分別為8 572、10 160、10 812、11 449、9 502、9 056和8 651 t/(hm2·a);2010年,土壤侵蝕強(qiáng)度呈現(xiàn)極強(qiáng)烈的縣為甘肅省白銀區(qū)、平川區(qū)、甘谷縣、環(huán)縣和慶城縣,各縣平均土壤侵蝕模數(shù)分別為8 907、8 622、9 104、9 873和8 045 t/(hm2·a);2020年,土壤侵蝕強(qiáng)度呈現(xiàn)極強(qiáng)烈的縣為寧夏回族自治區(qū)大武口區(qū)、惠農(nóng)區(qū)和內(nèi)蒙古自治區(qū)烏拉特后旗(流域內(nèi)部分區(qū)域),平均土壤侵蝕模數(shù)分別為14 575、8 329和13 288 t/(hm2·a)??h域尺度下,土壤侵蝕熱點(diǎn)分析結(jié)果(圖3)與土壤侵蝕強(qiáng)度呈現(xiàn)相似的時(shí)空特征,侵蝕聚集于黃土高原地區(qū),寧夏回族自治區(qū)大武口區(qū)、惠農(nóng)區(qū)和內(nèi)蒙古自治區(qū)烏拉特后旗逐漸成為新的侵蝕聚集中心。
土壤侵蝕受自然因素(降水、坡度、土壤類型等)和人類活動(dòng)(LUCC、人口密度、植被覆蓋度)的共同影響,因此確定自然因素與人類活動(dòng)對土壤侵蝕的解釋力對制定后續(xù)的水土保持措施尤為重要。黃河流域土壤侵蝕驅(qū)動(dòng)因素地理探測結(jié)果如圖4所示。
總體來說,各因子對土壤侵蝕的解釋力均有限,其中植被覆蓋度因子對流域土壤侵蝕的解釋力最強(qiáng)(2000年q值最大,為0.160 5),海拔、土壤類型與年平均降雨次之,另外,兩個(gè)人為因素(LUCC和人口密度)的解釋力均很有限。因子交互作用對土壤侵蝕的解釋力略強(qiáng)于單一因子,但最多也僅在2000年坡度與植被覆蓋度交互作用時(shí),解釋了約32%的土壤侵蝕。植被覆蓋度作為主要的控制因素,與其他因素交互作用時(shí)q值相對較大。在時(shí)間上,土壤侵蝕因子探測的結(jié)果總體呈減小趨勢,植被覆蓋度q值從0.160 5減小至0.125 5,人口密度q值從0.096 9減小至0.009 9,LUCC的q值從0.064 5減小至0.047 0。因子交互作用探測的結(jié)果變化趨勢相似,但植被覆蓋度與地貌類型交互時(shí)的q值呈增加趨勢,植被覆蓋度與LUCC交互時(shí)的q值呈先增加后減少的趨勢。驅(qū)動(dòng)因素的變化趨勢表明,黃河流域土壤侵蝕狀況在好轉(zhuǎn)的同時(shí),驅(qū)動(dòng)機(jī)理變得更加復(fù)雜。
圖3 熱點(diǎn)分析結(jié)果圖Fig.3 Map of hot spot analysis result 注:此圖基于國家自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站審圖號(hào) 為GS2020(4619)的標(biāo)準(zhǔn)地圖制作,底圖無修改。
在縣域尺度下,對各縣土壤侵蝕驅(qū)動(dòng)因素因子探測q值求平均,結(jié)果如表2所示。在以縣域?yàn)閱卧獎(jiǎng)澐贮S河流域進(jìn)行土壤侵蝕驅(qū)動(dòng)因素分析后,土壤侵蝕的驅(qū)動(dòng)機(jī)理更加明顯,平均q值大大超過以流域?yàn)閱卧_展因子探測的q值,且LUCC這一人為因素成為第二重要的驅(qū)動(dòng)因素。這表明在一定程度內(nèi),自然因素對土壤侵蝕的影響能力隨尺度減小而逐漸減小,而人為因素對土壤侵蝕的影響隨尺度減小而增大,人為因素在局部區(qū)域發(fā)揮了更大作用。這一發(fā)現(xiàn)與Guo等[3]、Li等[5]的結(jié)論相似。
圖4 黃河流域地理探測結(jié)果圖Fig.4 Map of geographical exploration results in the Yellow River Basin
表2 因子探測結(jié)果平均值
縣域尺度下,各縣土壤侵蝕驅(qū)動(dòng)因素因子探測結(jié)果表明(圖5中,空值主要因地理探測結(jié)果中的sig大于0.05,無法通過顯著性檢驗(yàn)造成),黃河流域各縣土壤侵蝕驅(qū)動(dòng)因素的空間分異明顯,但總體上無明顯分布規(guī)律。在空間上,LUCC、地貌類型、土壤類型、降水、海拔和坡度在寧夏回族自治區(qū)北部部分區(qū)縣(平羅縣、賀蘭縣、西夏縣等)對土壤侵蝕的解釋力較強(qiáng);地貌類型與海拔在內(nèi)蒙古自治區(qū)部分區(qū)縣(土默特右旗、土默特左旗、東河區(qū)和賽罕區(qū))與山西省部分區(qū)縣(離石區(qū)、中陽縣和汾陽市等)對土壤侵蝕的解釋力較強(qiáng);人口密度在絕大多數(shù)區(qū)縣對土壤侵蝕的解釋力有限。此外,植被覆蓋度與其他因素相比,展現(xiàn)出不同的空間分布格局,其在山西省、陜西省和甘肅省對土壤侵蝕的解釋力較強(qiáng),而這些區(qū)縣的土壤侵蝕強(qiáng)度減弱最為明顯,表明植被恢復(fù)能有效減少土壤侵蝕。
除植被覆蓋度外,其他因素的q值在時(shí)間上總體呈現(xiàn)減小趨勢,這與流域土壤侵蝕驅(qū)動(dòng)因素的時(shí)間變化特征相似。為探究植被覆蓋度因子q值增強(qiáng)的原因,將植被覆蓋度因子探測結(jié)果與植被覆蓋度數(shù)據(jù)進(jìn)行疊加分析,發(fā)現(xiàn)2000—2020年,黃河流域植被覆蓋度從49%上升至59%,而植被覆蓋度因子在植被覆蓋度較高的區(qū)縣對土壤侵蝕的解釋力更強(qiáng)。又因在RUSLE模型中,當(dāng)植被覆蓋度大于78.3%時(shí)C因子被賦值為0,說明當(dāng)研究區(qū)植被覆蓋度高于特定值時(shí),植被控制土壤侵蝕的益處趨于穩(wěn)定,這與Teng等[37]的研究結(jié)論相似。
圖5 各縣因子探測結(jié)果圖Fig.5 Map of factor detection results for each county 注:此圖基于國家自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站審圖號(hào)為GS2020(4619)的標(biāo)準(zhǔn)地圖制作,底圖無修改。
為進(jìn)一步探究新的侵蝕聚集中心形成的驅(qū)動(dòng)因素,對寧夏回族自治區(qū)大武口區(qū)、惠農(nóng)區(qū)和內(nèi)蒙古自治區(qū)烏拉特后旗三個(gè)區(qū)域的因子探測結(jié)果進(jìn)行分析。由表3可知,2020年,新的侵蝕熱點(diǎn)的形成主要受坡度、地貌類型和LUCC影響。分析LUCC、地貌與遙感影像數(shù)據(jù)可知,三個(gè)區(qū)域的LUCC類型均以草地為主,大武口區(qū)與惠農(nóng)區(qū)地處賀蘭山丘陵地區(qū),烏拉特后旗地處陰山山脈區(qū)域,均為土壤侵蝕易發(fā)區(qū),加之人類活動(dòng)的影響,極易造成土壤侵蝕的發(fā)生。Li等[4]的研究發(fā)現(xiàn),黃河上游凈土壤侵蝕速率的空間格局主要由坡度控制,這也從一定程度上印證了所得結(jié)論。
表3 大武口區(qū)、惠農(nóng)區(qū)和烏拉特后旗因子探測結(jié)果Tab.3 Factor detection results of Dawukou, Huinong and Urad Rear Banner
由于不同學(xué)者采用的數(shù)據(jù)源不同,而且RUSLE模型的因子計(jì)算方法也不盡相同,因此,各文獻(xiàn)對黃河流域或相關(guān)區(qū)域土壤侵蝕模數(shù)的估算也存在較大差異。此外,考慮到黃河流域土壤侵蝕情況較為復(fù)雜,從水利部《中國河流泥沙公報(bào)》僅能獲知黃河部分水文站點(diǎn)的輸沙模數(shù),而在泥沙輸移比未知時(shí)無法估算侵蝕模數(shù),故對比了當(dāng)前研究與其他研究來驗(yàn)證本文結(jié)果的準(zhǔn)確性(表4)。當(dāng)前研究結(jié)果與其他研究結(jié)果具有較好的線性擬合(R2=0.96),但模型精度仍有待進(jìn)一步提高。
表4 其他研究的土壤侵蝕模數(shù)Tab.4 Soil erosion modulus obtained from previous studies
綜上,黃河流域土壤侵蝕時(shí)空特征明顯,水土流失嚴(yán)重的區(qū)域主要位于黃土高原地區(qū),流域土壤侵蝕狀況有所好轉(zhuǎn),但近年來陰山山脈與賀蘭山區(qū)的水土流失較為嚴(yán)重。各因子對黃河流域土壤侵蝕的解釋力均有限,植被覆蓋度始終是黃河流域土壤侵蝕的重要控制因素,該結(jié)論與周璐紅[44]、賈磊[45]等人的分析結(jié)果一致。黃河流域各縣土壤侵蝕驅(qū)動(dòng)因素的空間分異明顯,但空間分布規(guī)律不顯著。
不同季節(jié)的土壤侵蝕及其主要控制因素不同,Li等[4]研究發(fā)現(xiàn),黃河上游6~8月土壤水蝕最為強(qiáng)烈,且NDVI在夏季對土壤侵蝕的解釋力更強(qiáng),降水在春秋季節(jié)、河源地區(qū)對土壤侵蝕的解釋力更強(qiáng)。在后期研究中,還應(yīng)探討不同季節(jié)、不同地貌特征(梯田等)下的黃河流域土壤侵蝕及驅(qū)動(dòng)因素的時(shí)空變化,以期為黃河流域水土保持精細(xì)化管理提供幫助。
1) 2000年、2010年和2020年黃河流域平均土壤侵蝕模數(shù)分別為1 994、1 860和1 384 t/(km2·a),侵蝕熱點(diǎn)呈東北—西南走向,并呈帶狀集中于黃土高原地區(qū)。在坡度、地貌類型與LUCC的影響下,陰山西部與賀蘭山區(qū)逐漸成為新的侵蝕熱點(diǎn)區(qū)域。
2) 植被覆蓋度為黃河流域土壤侵蝕最重要的控制因素。在縣域尺度下,LUCC成為僅次于植被覆蓋度的影響因素,人為因素在局部區(qū)域發(fā)揮了更大作用。各縣土壤侵蝕驅(qū)動(dòng)因素差異性特征明顯,但除植被覆蓋度在高植被覆蓋區(qū)域解釋力較強(qiáng)外,其他驅(qū)動(dòng)因素?zé)o明顯空間分布特征。
3) 2000—2020年,黃河流域土壤侵蝕狀況好轉(zhuǎn),但侵蝕驅(qū)動(dòng)因素的解釋力呈現(xiàn)減弱趨勢,土壤侵蝕驅(qū)動(dòng)機(jī)理更加復(fù)雜。