胡夢甜,張 慧,2①,喬亞軍,劉 坤,王 智②,徐網(wǎng)谷

(1.生態(tài)環(huán)境部南京環(huán)境科學研究所,江蘇 南京 210042;2.南京信息工程大學江蘇省大氣環(huán)境與裝備技術協(xié)同創(chuàng)新中心/ 地理科學學院,江蘇 南京 210044)

風蝕對地表土壤的大量搬運和堆積,是導致干旱半干旱區(qū)土地沙化和沙地化進程最重要、最直接的作用過程之一[1-2]。土壤風蝕研究方法主要有野外調查觀測、風洞模擬、元素示蹤法、遙感和GIS等[3],但這些方法無法滿足大范圍區(qū)域風蝕量快速動態(tài)估算的需求[4]。為了定量模擬大范圍區(qū)域風蝕量,自20世紀60年代起,國外先后開發(fā)了大量的風蝕模型,主要包括風蝕方程(wind erosion equation,WEQ)[5]、修正風蝕方程(revised wind erosion equation,RWEQ)[6]、德克薩斯模型(taxas erosion analysis model,TEAM)、風蝕評價模型(wind erosion assessment model,WEAM)[7]、風蝕預報系統(tǒng)(wind erosion prediction system,WEPS)[8]等。其中,WEPS系統(tǒng)較為先進,但該系統(tǒng)建模過程復雜,數(shù)據(jù)要求繁雜,且缺乏相關參數(shù)的本地化工作,故該模型在國內應用較少[9]。相比先進的WEPS系統(tǒng),RWEQ模型因其參數(shù)較易獲取、操作更為簡單的特點,在我國不少地區(qū)的土壤風蝕和防風固沙功能計算與評估中得到了廣泛應用,如在內蒙古自治區(qū)錫林郭勒盟、青海省渾善達克沙地、黑河流域[10-11]等地都開展了相關工作。該模型也被作為《生態(tài)保護紅線劃定指南》和《全國生態(tài)狀況評估技術規(guī)范》中防風固沙功能重要性的評估方法。

對于風蝕量變化的驅動因素,學者大多采用相關性分析、主成分或聚類分析等傳統(tǒng)方法進行研究[12],缺少從地理分異的角度對風蝕量變化的定量歸因和驅動因素的空間差異性研究。地理探測器是一種強大的、能直接量化驅動因素及其交互作用影響的方法,它不必遵循傳統(tǒng)統(tǒng)計方法的假設,且不涉及復雜的參數(shù)設置過程[13]。近年來,地理探測器在土地利用[14]、生態(tài)服務功能[15]等地理現(xiàn)象的空間驅動力分析領域應用廣泛。

呼倫貝爾森林-草原生態(tài)交錯帶是我國生態(tài)系統(tǒng)結構保存完整、健康狀況良好的林草交錯帶,是我國北方地區(qū)重要的生態(tài)屏障[16]。研究區(qū)冬春季寒冷多風,該區(qū)的草原是以栗鈣土、風沙土為主的干草原,易形成風蝕[17],因此探究該地區(qū)土壤風蝕量的變化及其驅動因素對維護區(qū)域生態(tài)安全具有重要意義。近年來,雖然有學者研究了呼倫貝爾森林-草原生態(tài)交錯帶土壤風蝕,但主要集中于區(qū)域的生態(tài)效益評估[18-19],鮮有對區(qū)域風蝕量變化的驅動力開展相關研究。筆者利用RWEQ風蝕修正模型研究呼倫貝爾森林-草原生態(tài)交錯帶2000、2010和2018年的土壤風蝕量時空分布變化,并識別其驅動因素,以期為該區(qū)域科學防治土壤風蝕、遏制區(qū)域沙化趨勢提供科學支撐。

1 研究區(qū)概況

呼倫貝爾森林-草原生態(tài)交錯帶處于大興安嶺西麓山地向亞洲中部蒙古高原東北部過渡的區(qū)域,東北區(qū)域為林區(qū),海拔700～1 700 m,由東北向西南依次為農田、森林草原、草甸草原和干旱草原,海拔550～1 000 m。研究區(qū)位于內蒙古自治區(qū)呼倫貝爾市中部,地理位置處于北緯46°10′～53°26′,東經117°33′～122°55′ 之間,行政區(qū)域涉及呼倫貝爾市的額爾古納市、根河市、牙克石市、陳巴爾虎旗、海拉爾區(qū)、鄂溫克自治旗和新巴爾虎左旗共7個旗市(圖1)。該地區(qū)的經濟活動主要包括種植業(yè)和畜牧業(yè),其中種植業(yè)以冬小麥、油菜和苜蓿等種植為主;畜牧業(yè)以養(yǎng)羊、牛為主。研究區(qū)處在溫帶-寒溫帶氣候區(qū),氣候較干燥,多大風,研究區(qū)年平均氣溫在 -2.2～2.4 ℃之間,年降水量為290～450 mm[16],年均蒸發(fā)量590～960 mm,年日照時數(shù)為2 600～2 800 h,年均風速1.8～2.45 m·s-1。研究區(qū)土地利用類型以森林、草地和濕地為主,總占比超過90%。其中,森林集中分布在研究區(qū)的東北部山區(qū);草地主要分布在西南部;農田集中分布在林草交錯帶,主要位于額爾古納市、呼倫貝爾市、牙克石市區(qū)域內;濕地主要分布在額爾古納河、海拉爾河、根河、輝河區(qū)域。

該圖基于審圖號為蒙S(2020)027號的標準地圖制作。

2 數(shù)據(jù)來源與研究方法

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

研究區(qū)采用的遙感數(shù)據(jù)源為2000和2010年的Landsat 5 TM衛(wèi)星影像(空間分辨率30 m)、2018年的Landsat 8 TM 衛(wèi)星影像(空間分辨率30 m)(https:∥glovis.usgs.gov/)。依據(jù)《全國30米分辨率土地利用分類系統(tǒng)》并結合研究區(qū)實際情況,考慮到大興安嶺地區(qū)在1987年發(fā)生過特大火災,2000年過火林區(qū)的植被還未恢復,因此在土地利用分類中增加火燒跡地類型,將土地利用類型分為森林、草地、農田、濕地、城鎮(zhèn)用地、沙地和火燒跡地7類。采用人工目視解譯方法分類,精度達95%以上,土地利用變化檢測總體精度達85%以上。氣象數(shù)據(jù)包括月平均風速,數(shù)據(jù)來源于國家青藏高原科學數(shù)據(jù)中心(http:∥data.tpdc.ac.cn/zh-hans/data),空間分辨率為1 km。雪蓋因子采用中國雪深長時間序列數(shù)據(jù)集(http:∥data.tpdc.ac.cn/zh-hans/data/df40346a-0202-4ed2-bb07-b65dfcda9368/);土壤特性因子數(shù)據(jù)來自世界土壤數(shù)據(jù)庫(HWSD,http:∥webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database),包括土壤粗砂、細砂、黏粒及有機質含量等理化性質;NDVI數(shù)據(jù)采用資源環(huán)境科學與數(shù)據(jù)中心發(fā)布的數(shù)據(jù)產品(http:∥www.resdc.cn/),空間分辨率為1 km;土壤濕度數(shù)據(jù)采用干旱植被指數(shù)法計算得到,其中地表溫度采用NASA發(fā)布的MOD11A2(https:∥ladsweb.modaps.eosdis.nasa.gov)數(shù)據(jù);地表糙度因子采用GDEM DEM 30 m數(shù)字高程數(shù)據(jù)(http:∥www.gscloud.cn/)計算得到。將所有因子統(tǒng)一為250 m分辨率的柵格數(shù)據(jù),保證柵格計算過程準確。

2.2 風蝕量計算模型

風蝕量計算采用RWEQ模型, RWEQ作為美國農田風蝕模型,其計算參數(shù)均來源于美國本土[20]。遲文峰等[21]在內蒙古高原利用137Cs示蹤技術檢驗RWEQ模型的模擬效果,發(fā)現(xiàn)該模型的風蝕量計算結果擬合性較好(R2=0.83,P<0.01)?！渡鷳B(tài)保護紅線劃定指南》[22]將該模型作為風蝕量和防風固沙功能的推薦評估方法,RWEQ方法計算公式為

(1)

S=150.71×(FW×FE×FSC×K′×C)-0.371 1,

(2)

Qmax=109.8×FW×FE×FSC×K′×C。

(3)

式(1)～(3)中,SL為實際土壤侵蝕量,t·km-2·a-1;S為區(qū)域防風固沙系數(shù);Qmax為風沙滯留量,kg·m-1;Z為最大風蝕出現(xiàn)距離,m;FW為氣候因子,kg·m-1;FE為土壤可蝕性因子;FSC為土壤結皮因子;K′為地表粗糙度因子;C為植被因子。

氣候因子(FW)計算公式為

(4)

fW=u2×(u2-u1)2×Nd,

(5)

(6)

ρ=348×(1.013-0.118 3×LE+0.004 8×LE2)/T。

(7)

式(4)～(7)中,fW為風力因子,m3·s-3;ρ為空氣密度,kg·m-3;g為重力加速度,取9.8 m·s-2;WS為各月土壤濕度因子,表征土壤濕度抑制土壤風蝕效果的大小,土壤濕度越大,土壤濕度因子越小,越不易起沙,該因子以往大多采用氣象站點數(shù)據(jù)插值計算得到,由于氣象站點分布不均,插值結果往往有“牛眼”存在,因此采用干旱植被指數(shù)法表征土壤濕度因子[23-24];DS為雪蓋因子;u1為起沙風速,參照江凌[25]的計算方法,取5 m·s-1;j為1個月內日平均風速≥u1的天數(shù),j=1,2,…,m;uj為第j天的日平均風速,m·s-1;TLS,i為i評估區(qū)域的地表溫度,℃;TLS,dry為評估區(qū)域NDVI對應的最高地表溫度,℃,即干邊;TLS,wet為評估區(qū)域NDVI對應的最低地表溫度,℃,即濕邊;LE為海拔高度,km;T為絕對溫度,K,即在各月平均氣溫數(shù)據(jù)t的基礎上加常數(shù)273.15。

土壤可蝕性因子(FE)和土壤結皮因子(FSC)計算公式為

(8)

(9)

式(8)～(9)中,wsa為土壤粗砂含量,%;wsi為土壤粉砂含量,%;wcl為土壤黏粒含量,%;wOM為土壤有機質含量,%;wCa為碳酸鈣含量,%,此次計算未予考慮,其值取0。

植被覆蓋因子(C)計算公式為

C=e-CS×ai。

(10)

式(10)中,CS為植被覆蓋度,%,由12個月的NDVI指數(shù)計算得到年均植被覆蓋度;ai為不同植被類型的系數(shù),因不同植被類型的土壤風蝕效果不同,森林、草地、農田、沙地(包括裸地和沙地)分別取值0.153 5、0.115 1、0.043 8、0.071 3[25]。

地表糙度因子(K′)計算公式為

K′=e1.86Kr-0.127Crr-2.41Kr0.934,

(11)

(12)

式(11)～(12)中,Kr為地形粗糙度長度,m;Crr為隨機糙度,在區(qū)域尺度的計算中可以忽略不計;L為地勢起伏參數(shù),m;△H為距離L范圍內的海拔高程差,m。

2.3 驅動因素選擇

土壤風蝕的產生受到土壤、地形等自然本底因素的制約,同時氣溫、降水、風速等氣候因素以及人類活動也會對其產生影響[10]。根據(jù)研究區(qū)的實際情況,從氣候因素和人類活動方面選取5個影響因子,包括降水變化(X1)、氣溫變化(X2)、植被覆蓋度變化(X3)、風速變化(X4)和土地利用類型(X5),作為探測研究區(qū)土壤風蝕變化量的驅動因素。以研究區(qū)范圍為基礎創(chuàng)建漁網(wǎng),參照各影響因子空間分辨率將漁網(wǎng)大小設置為1 000 m×1 000 m。采用ArcGIS軟件Spatial Analysis Tool工具包中的Zonal Statistic工具,將各影響因子按照平均值統(tǒng)計到各個漁網(wǎng)。因為地理探測器中自變量為類型量,對于順序量需要進行離散化[13]。在R語言環(huán)境下,采用等間隔法、分位數(shù)法、自然斷點法和標準差法比較降水變化、氣溫變化、植被覆蓋度變化、風速變化4個因子的離散化效果,采用自然斷點法將植被覆蓋度變化、降水變化、氣溫變化、風速變化分別分為8、6、9、9類。土地利用之間轉化類型達38類,因此將土地利用變化因子分為38類。為深入探究呼倫貝爾森林-草原生態(tài)交錯帶土壤風蝕量時空變化驅動機制,針對2000—2018年土壤風蝕量變化的主要影響因素開展研究。

2.4 地理探測器模型

采用地理探測器模型進行土壤風蝕量侵蝕變化的驅動力分析。地理探測器模型主要包括因子探測器和交互探測器[26]。因子探測器用于分析自變量對因變量的解釋程度,解釋力的強弱通過比較q值的大小反映[13],q取值在0～1之間,q值越大,表示該影響因子對土壤風蝕變化量的影響越大。交互作用探測器可探索2個自變量的聯(lián)合效應是否會增加、減少對因變量的解釋力[27]。通過比較因子單獨作用時的q值〔q(x1)和q(x2)〕與交互作用時的q值〔q(x1∩x2)〕,對2個因子之間的關系進行界定。q(x1∩x2)>q(x1)+q(x2)、q(x1∩x2)同時大于q(x1)和q(x2)、q(x1∩x2)處于q(x1)和q(x2)之間、q(x1∩x2)同時小于q(x1)和q(x2)、q(x1∩x2)=q(x1)+q(x2)分別表征非線性增強、雙因子增強、單因子非線性減弱、非線性減弱、獨立。

3 結果與分析

3.1 呼倫貝爾森林-草原生態(tài)交錯帶土壤風蝕量及強度的時空分布變化

從風蝕量的年際變化來看,2000、2010和2018年研究區(qū)土壤風蝕總量分別為9.74×107、1.33×108、8.51×107t,總體呈現(xiàn)先上升后下降的趨勢,2000—2018年土壤風蝕總量減少1.23×107t,平均年降速為68.33萬t,年均下降率為0.70%(圖1)。

評估期內,研究區(qū)單位面積風蝕量在7.03～10.88 t·hm-2之間波動,單位面積風蝕量的年均下降率為0.54%,研究區(qū)單位面積風蝕量的下降幅度低于土壤風蝕總量的下降幅度。

從研究區(qū)土壤風蝕量及強度的空間分布來看,土壤風蝕量整體呈現(xiàn)由東北部林區(qū)向西南部草原區(qū)逐漸增加的特征(圖1)。依據(jù)《土壤侵蝕分類分級標準》將研究區(qū)土壤風蝕強度分為5級(劇烈、極強烈、強烈、中度、輕度和微度)(表1)。結果表明,2018年劇烈風蝕區(qū)域面積占研究區(qū)總面積的7.61%,該區(qū)域年均單位面積風蝕量達378.91 t·hm-2,主要包括新巴爾虎左旗西部、陳巴爾虎旗南部和鄂溫克族自治旗西部區(qū)域;土壤極強烈風蝕區(qū)域面積占研究區(qū)總面積的5.86%,該區(qū)域年均單位面積風蝕量達105.91 t·hm-2,集中分布在新巴爾虎左旗中部,零散分布在鄂溫克族自治旗;土壤強烈和中度風蝕區(qū)域面積占比相近,分別為4.13%和3.97%,分布于新巴爾虎左旗、陳巴爾虎旗和鄂溫克族自治旗交界處;研究區(qū)有13.86%的土壤屬輕度風蝕區(qū),主要位于林草過渡帶;土壤微度風蝕區(qū)面積占研究區(qū)總面積的64.57%,主要分布在林區(qū)。

表1 2000—2018年呼倫貝爾森林-草原生態(tài)交錯帶土壤風蝕變化

從研究區(qū)風蝕量及強度空間變化來看,2000—2018年,研究區(qū)有90.86%的區(qū)域土壤風蝕強度未發(fā)生變化;有9.04%的區(qū)域土壤風蝕強度減輕(表1),主要分布在西南部的新巴爾虎左旗、陳巴爾虎旗和鄂溫克族自治旗,單位面積風蝕量減少了27 t·hm-2;有0.10%的區(qū)域土壤風蝕強度加重,主要分布在新巴爾虎器西北角和鄂溫克族自治旗正北角,單位面積風蝕量增加了20 t·hm-2。可見,2000—2018年研究區(qū)土壤風蝕程度整體較為穩(wěn)定,部分地區(qū)加重(圖2),其中西南部地區(qū)土壤風蝕程度減輕,而新巴爾虎器西北角和鄂溫克族自治旗正北角部分地區(qū)土壤風蝕強度惡化。

該圖基于審圖號為蒙S(2020)027號的標準地圖制作。

3.2 土壤風蝕量變化的驅動力分析

依據(jù)地理探測器分析得出,單個影響因子變化量對土壤風蝕量變化量的解釋力排序為土地利用變化>降水變化>風速變化>植被覆蓋度變化>氣溫變化。整體來看,各因子變化量對土壤侵蝕量的解釋力q值均很小,最大值亦不超過0.1,表明單一因子變化量對土壤風蝕量變化驅動作用有限。與風速變化和氣溫變化因子相比,降水變化對土壤風蝕量變化的影響力較高,說明降水量增加對區(qū)域風蝕量減少發(fā)揮著重要作用[11]。進一步探究各因子變化量的交互作用,發(fā)現(xiàn)任意2種因子交互作用的解釋力高于單個因子(表2),且降水變化與其他因子變化均呈雙因子非線性增強作用。

表2 各驅動因子對土壤風蝕量變化量的解釋力(q值)

其中,降水變化和土地利用變化的協(xié)同作用對土壤風蝕量變化的解釋力最大,解釋力q值達0.22。同時,植被覆蓋度變化與其他因子變化也均呈雙因子非線性增強作用,且植被覆蓋度變化協(xié)同降水變化解釋力最強,q值為0.14;植被覆蓋度變化協(xié)同土地利用變化次之,q值為0.13。

3.3 土壤風蝕量變化的人為影響因素分析

利用研究區(qū)2000—2018年土地利用現(xiàn)狀圖和土地利用轉換圖(圖3～4),定量分析土地利用變化對風蝕量的影響。2000—2018年大興安嶺林草交錯帶各土地利用類型發(fā)生了復雜的相互轉換,農田轉出面積變大,草地面積增加,濕地面積萎縮,退耕還草和城市化是農田的主要轉出方向。2000—2010年,退耕還草是農田面積減少、草地面積增加的主要轉換方式,退耕還草面積達到201.71 km2,以鄂溫克族自治旗林草交錯帶內的農田轉換為主。2010年以后,城鎮(zhèn)用地增加成為農田主要的轉出方向,期間城鎮(zhèn)用地增加占用農田40.35 km2,城鎮(zhèn)用地占用農田擴張趨勢明顯。2000—2018年間,濕地面積凈減少117.66 km2,尤其是2000—2010年,濕地大面積退化為草地,多發(fā)生在新巴爾虎旗境內的呼倫湖、海拉爾流域(表3)。2000—2018年,土地利用類型轉換分別造成了24.32 t風蝕量的增加和12.31×106t風蝕量的減少。

表3 研究區(qū)2000—2018年土地利用轉移面積表

該圖基于審圖號為蒙S(2020)027號的標準地圖制作。

地理探測器研究表明,土地利用變化協(xié)同植被覆蓋度變化將增大對土壤風蝕量變化的影響,因此將土地利用變化與植被覆蓋度變化協(xié)同分析土壤風蝕量變化的人為影響因素。將植被覆蓋度按照馬志勇等[28]提出的植被覆蓋度分級標準,按0～0.45、>0.45～0.75和>0.75～1.00分別將研究區(qū)植被覆蓋度分為低、中、高3類,研究不同土地利用類型、不同植被蓋度之間的轉化對土壤風蝕量的影響(植被覆蓋度分類中濕地為沼澤濕地,不包含河流和湖泊)。

表4列出了居前16位的不同植被覆蓋度的土地利用類型下土壤風蝕量變化,研究區(qū)因為草地覆蓋度增加減少了8.53×106t的土壤風蝕量,主要分布在陳巴爾虎旗、新巴爾虎左旗和鄂溫克族自治旗大部分地區(qū);因草地覆蓋度降低增加了24.32 t土壤風蝕量,主要分布在鄂溫克族自治旗西部巴彥烏拉蘇木地區(qū)。沙地轉化為草地、濕地、森林后土壤風蝕量減少了2.38×106t,主要分布在新巴爾虎旗東部、陳巴爾虎旗中部、鄂溫克族自治旗中西部這3條沙帶區(qū)[29];濕地覆蓋度提高后,土壤風蝕量減少了7.34×105t,主要分布在額爾古納河、海拉爾河等沼澤濕地區(qū)域。

表4 研究區(qū)不同植被覆蓋度的土地利用類型下土壤風蝕量變化

該圖基于審圖號為蒙S(2020)027號的標準地圖制作。

草地覆蓋度增加、沙化土地封育(沙地質量改善、沙地轉草地、沙地轉濕地、沙地轉森林)、生態(tài)退耕(農田轉草地、林地和濕地)、天然林保護(森林覆蓋度增加)對風蝕量減少的貢獻率分別為69.32%、19.37%、0.06%和1.81%,占土壤風蝕減少總量的90.56%。其中,草地覆蓋度增加導致的風蝕量降低最為明顯,其次是沙化土地封育,雖然草地覆蓋度增加、沙化土地封育與研究區(qū)降水增加有一定關系,但是這兩者與研究區(qū)的圍封禁牧、沙化土地封育政策關系更加密切,這些保護政策增加了草地生物量、蓋度及高度,減少了草地的風蝕作用,提高了草地的保水能力[30]。因此,今后應繼續(xù)采取圍封禁牧、休牧、輪牧、改良牧草場等措施,增加其覆蓋度,可以有效減少該區(qū)域的風蝕量。

4 結論

(1)研究區(qū)土壤風蝕量整體呈現(xiàn)東北部林區(qū)向西南部草原區(qū)遞增的空間分布特征。2000—2018年土壤風蝕量波動下降,土壤風蝕總量共減少1.23×107t,年均下降率為0.70%;從風蝕強度變化來看,2000—2018年研究區(qū)有90.86%的區(qū)域土壤風蝕強度未發(fā)生變化;有9.04%的區(qū)域土壤風蝕強度減輕;有0.10%的區(qū)域土壤風蝕強度惡化。

(2)2000—2018年間,研究區(qū)土壤風蝕量變化的主要驅動因子解釋力表現(xiàn)為土地利用變化>降水變化>風速變化>植被覆蓋度變化>氣溫變化,整體來看,各因子變化量解釋力q值均很小,表明單一因子變化量對土壤風蝕量變化的驅動作用有限。對各因子變化量的交互作用分析發(fā)現(xiàn),降水量變化與其他因子變化均呈雙因子非線性增強作用,其中,降水變化協(xié)同土地利用變化與降水變化協(xié)同植被覆蓋度變化可顯著增強對土壤風蝕量變化的驅動作用。

(3)研究區(qū)草地覆蓋度增加、沙化土地封育、生態(tài)退耕、天然林保護等轉換方式對土壤風蝕量降低的貢獻率為90.56%,說明該區(qū)域實施的沙化土地封育、天然林保護、退耕還草、退耕還林等一系列生態(tài)保護措施,對該地區(qū)的生態(tài)環(huán)境改善具有重要的推動作用。今后還需進一步采取禁牧、休牧、輪牧等措施,恢復和提升草地覆蓋度,從而更有效地減少區(qū)域風蝕量。