黃孟冬,肖 玉,*,秦克玉,甘 爽,謝高地,劉婧雅,王洋洋,牛櫻楠,劉 佳

1 中國科學(xué)院地理科學(xué)與資源研究所,北京 100101 2 中國科學(xué)院大學(xué),北京 100049

風(fēng)蝕是土壤顆粒在風(fēng)力作用下侵蝕、搬運(yùn)和堆積的過程,也是土地退化和荒漠化的主要原因[1]。我國北方干旱半干旱地區(qū)是受風(fēng)蝕現(xiàn)象影響較為嚴(yán)重的地區(qū),是沙漠化防治的重點(diǎn)區(qū)域[2]。防風(fēng)固沙服務(wù)作為風(fēng)蝕地區(qū)生態(tài)系統(tǒng)提供的一項(xiàng)重要防護(hù)型服務(wù),對風(fēng)蝕發(fā)生地與周邊地區(qū)風(fēng)沙災(zāi)害治理、生態(tài)環(huán)境恢復(fù)具有重要意義。渾善達(dá)克地區(qū)位于我國內(nèi)蒙古高原干旱半干旱地區(qū),是我國北方農(nóng)牧交錯帶,生態(tài)環(huán)境脆弱,腹地有渾善達(dá)克沙地,是京津地區(qū)風(fēng)沙主要來源地[2—3]。2000年以來,隨著京津風(fēng)沙源治理工程、三北防護(hù)林工程、退耕還林還草、退牧禁牧政策的實(shí)施,渾善達(dá)克地區(qū)生態(tài)系統(tǒng)狀況呈現(xiàn)改善趨勢,防風(fēng)固沙服務(wù)供給能力也有一定的提升,風(fēng)蝕現(xiàn)象有所好轉(zhuǎn)[4—6]。然而,氣候變化、人類活動與生態(tài)系統(tǒng)之間頻繁的相互影響,使得環(huán)境敏感區(qū)的生態(tài)系統(tǒng)結(jié)構(gòu)并不穩(wěn)定,其提供的生態(tài)系統(tǒng)服務(wù)也相對有限[7—8]。因此,開展長時間序列的渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)時空變化及驅(qū)動因素分析對于區(qū)域生態(tài)系統(tǒng)恢復(fù)以及生態(tài)建設(shè)、土地管理具有重要意義。

防風(fēng)固沙服務(wù)一般通過潛在風(fēng)蝕量和實(shí)際風(fēng)蝕量之差計(jì)算獲得。潛在風(fēng)蝕量和實(shí)際風(fēng)蝕量可通過實(shí)地觀測和模型模擬確定[9—12]。Hu等[9]利用137Cs和210Pb示蹤法模擬渾善達(dá)克南部地區(qū)的風(fēng)蝕狀況。然而,這種方法需要花費(fèi)巨大的人力、物力,且僅能代表風(fēng)蝕測定樣點(diǎn)小范圍的結(jié)果,難以用于區(qū)域風(fēng)蝕量確定?；谀Ｐ偷娘L(fēng)蝕量估算能顯著提高計(jì)算效率,且可用于大面積區(qū)域的風(fēng)蝕量計(jì)算[13—14]。修正風(fēng)蝕方程模型(Revised Wind Erosion Equation,RWEQ)因子參數(shù)相對全面、模型構(gòu)成較為簡單、數(shù)據(jù)更易獲取[13—15],是應(yīng)用最為廣泛的防風(fēng)固沙服務(wù)計(jì)算方法。目前,已經(jīng)有較多基于RWEQ模型的防風(fēng)固沙量的模擬研究。這些區(qū)域防風(fēng)固沙服務(wù)模擬以多年研究為主[16—18],以避免單個年份的氣象數(shù)據(jù)帶來的隨機(jī)誤差。Li等[19]基于RWEQ模型分析了1990—2015年我國西北干旱地區(qū)土地利用變化對防風(fēng)固沙服務(wù)的影響。徐潔等[20]和張彪等[4]利用RWEQ分別分析了防風(fēng)固沙型重點(diǎn)生態(tài)功能區(qū)的防風(fēng)固沙服務(wù)能力和京津風(fēng)沙源工程對防風(fēng)固沙服務(wù)產(chǎn)生的效益。王俊枝等[21]通過RWEQ模型模擬了內(nèi)蒙古自治區(qū)境內(nèi)的渾善達(dá)克生態(tài)功能區(qū)1990—2015年防風(fēng)固沙服務(wù)的時空分布特征。長時間序列的防風(fēng)固沙服務(wù)模擬有助于全面分析一個區(qū)域防風(fēng)固沙服務(wù)的真實(shí)狀況及其時空格局變化,但現(xiàn)有的研究主要集中在1990年以后的長時間序列,更長時間序列的研究還不多見。

由于風(fēng)蝕現(xiàn)象受到自然和人為等多種因素綜合作用影響,開展防風(fēng)固沙影響因素分析可以為提高防風(fēng)固沙服務(wù)提供科學(xué)依據(jù)。目前,大部分研究采用線性相關(guān)性模型[22—24]、結(jié)構(gòu)方程模型[25]、約束效應(yīng)[26—27]等方法分析植被覆蓋度、風(fēng)速、降水、畜牧業(yè)發(fā)展、土地利用等因素與防風(fēng)固沙服務(wù)之間的關(guān)系[28—29]。關(guān)于渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)驅(qū)動因素分析的分析相對較少。申陸等[24]通過主成分分析法探究渾善達(dá)克生態(tài)功能區(qū)單位面積防風(fēng)固沙量變化的驅(qū)動因素。高曉霞等[23]則通過Logistic回歸模型探究渾善達(dá)克沙地動態(tài)變化的驅(qū)動因素。雖然上述研究為探討防風(fēng)固沙服務(wù)時空格局變化機(jī)制提供很好的參考,但多數(shù)研究為時間序列的研究,而且忽略了多因素之間的共同作用對防風(fēng)固沙服務(wù)的影響。從空間分布格局變化的角度探究防風(fēng)固沙服務(wù)變化驅(qū)動機(jī)制的研究也相對較少,類型變量對防風(fēng)固沙服務(wù)的影響也難以衡量。

地理探測器是一種基于空間分異性,揭示某一變量驅(qū)動機(jī)制的空間統(tǒng)計(jì)學(xué)方法[30],它能夠在探測類型與數(shù)值兩種不同屬性的變量對因變量空間分布影響的同時,分析自變量兩兩之間的交互作用及具體作用類型。目前,在資源環(huán)境、社會經(jīng)濟(jì)等多領(lǐng)域的空間格局變化過程分析中得到了廣泛的應(yīng)用[31—34],但用于分析防風(fēng)固沙驅(qū)動因素的研究還較少。應(yīng)用地理探測器分析防風(fēng)固沙服務(wù)空間格局變化的驅(qū)動因素,能夠兼顧類型變量與數(shù)值變量對服務(wù)變化的影響,分析各因素之間的相互作用,彌補(bǔ)防風(fēng)固沙服務(wù)空間分布驅(qū)動因素分析方面的不足。

因此,本研究以渾善達(dá)克重點(diǎn)生態(tài)功能區(qū)為例,通過RWEQ模型模擬1980—2018年防風(fēng)固沙服務(wù)的時空變化趨勢,利用地理探測器開展防風(fēng)固沙服務(wù)空間變化的驅(qū)動因素分析,為渾善達(dá)克地區(qū)生態(tài)建設(shè)與恢復(fù)以及土地管理提供科學(xué)依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

渾善達(dá)克地區(qū)(40.7°N—45.5°N,111.1°E—118.5°E)地處我國內(nèi)蒙古高原中段,總面積16.39萬km2,跨越內(nèi)蒙古自治區(qū)與河北省。該區(qū)空間范圍依據(jù)國家重點(diǎn)功能區(qū)中的渾善達(dá)克沙漠化防治生態(tài)功能區(qū)范圍劃定,包括內(nèi)蒙古8旗1縣(克什克騰旗、阿巴嘎旗、蘇尼特左旗、蘇尼特右旗、太仆寺旗、鑲黃旗、正鑲白旗、正藍(lán)旗、多倫縣)與河北省6縣(豐寧滿族自治縣、圍場滿族蒙古族自治縣、張北縣、康?？h、沽源縣、尚義縣)(圖1)。渾善達(dá)克地區(qū)地勢西南高、東北低,平均海拔1300 m。該地區(qū)屬于中緯度干旱、半干旱大陸性季風(fēng)氣候,年平均氣溫0—3℃,大部分地區(qū)年降水量在250—400 mm,年平均相對濕度為60%。春秋季風(fēng)沙多,年均風(fēng)速度為4—5 m/s,最大風(fēng)速普遍在24—28 m/s,年大風(fēng)日數(shù)(>8級)大約60—80 d。草地為主要的土地覆被類型,自東向西依次為草甸、草原、荒漠草原；土壤類型為栗鈣土、棕鈣土,非地帶性土壤為風(fēng)沙土。渾善達(dá)克地區(qū)整體生態(tài)環(huán)境脆弱,是我國北方重要生態(tài)屏障之一,同時也是京津地區(qū)主要沙源地,沙漠化是該地區(qū)主要的生態(tài)問題。針對土地退化、沙漠化等問題,2000年以來，開始在渾善達(dá)克地區(qū)實(shí)施京津風(fēng)沙源治理工程,2010年渾善達(dá)克地區(qū)被列為國家沙漠化防治生態(tài)功能區(qū)。在各類生態(tài)工程支持下,渾善達(dá)克地區(qū)通過圍欄封育、劃區(qū)輪牧、草地補(bǔ)播等措施實(shí)施草地恢復(fù)與治理,沙漠化趨勢逐漸減緩[35]。

圖1 研究區(qū)地理位置Fig.1 Location of study area

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

本研究使用的數(shù)據(jù)年份為1980年、1990年、1995年、2000年、2005年、2010年、2015年、2018年。20—20時日累計(jì)降水量、日均溫、10 m處日均風(fēng)速等氣象日值數(shù)據(jù)來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://data.cma.cn/)。土壤質(zhì)地?cái)?shù)據(jù)來自第二次全國土地調(diào)查,中國科學(xué)院南京土壤研究所提供的1:100萬土壤數(shù)據(jù)。雪蓋數(shù)據(jù)來源于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心(http://bdc.casnw.net)的“中國雪深長時間序列數(shù)據(jù)集(1979—2018年)”,空間分辨率為0.25°。2000年之前的歸一化植被指數(shù)(Normalized Difference Vegetation Index, NDVI)數(shù)據(jù)來源于GIMMS NDVIg3數(shù)據(jù)集,在Google Earth Engine上進(jìn)行圖像計(jì)算和提取等處理；2000年之后的NDVI數(shù)據(jù)(分辨率為1 km)、全國土地利用/覆被數(shù)據(jù)(空間分辨率為30 m)、GDP空間分布數(shù)據(jù)(分辨率為1 km)、1:100萬地貌類型、1:100萬土壤類型、1:100萬植被類型數(shù)據(jù),以及30 m高程數(shù)據(jù)均在中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(https://www.resdc.cn)下載。路網(wǎng)數(shù)據(jù)來源于OpenStreetMap(https://www.openstreetmap.org)。本研究所使用的年糧食播種面積、年糧食產(chǎn)量、年肉類產(chǎn)量、年末牲畜數(shù)量均來自各旗縣及其所屬盟、市的相關(guān)年份的統(tǒng)計(jì)年鑒并求得多年平均值,人工造林面積數(shù)據(jù)從《中國林業(yè)統(tǒng)計(jì)年鑒》中獲取。在模型運(yùn)算過程中,本研究所用柵格數(shù)據(jù)均采用Krasovsky_1940_Albers投影,并統(tǒng)一重采樣為30 m分辨率。

2 研究方法

2.1 防風(fēng)固沙服務(wù)物質(zhì)量模擬

在充分考慮氣候、地形、土壤理化性質(zhì)以及植被覆蓋等因素后,本研究采用修正風(fēng)蝕方程RWEQ(Revised Wind Equation)模型,定量模擬渾善達(dá)克地區(qū)風(fēng)蝕量的時空動態(tài)變化。根據(jù)RWEQ模型以及國內(nèi)實(shí)際情況進(jìn)行參數(shù)的調(diào)整[13, 36—38],計(jì)算公式如下：

(1)

Qmax=109.8(WF×EF×SCF×K′×C)

(2)

s=150.71(WF×EF×SCF×K′×C)-0.3711

(3)

(4)

Qrmax=109.8(WF×EF×SCF×K′)

(5)

sr=150.71(WF×EF×SCF×K′)-0.3711

(6)

G=SR-SL

(7)

式中,SL為實(shí)際風(fēng)蝕量(kg m-2a-1)；Qmax為實(shí)際風(fēng)力的最大輸沙能力(kg/m)；s表示實(shí)際關(guān)鍵地塊長度(m)；SR表示潛在風(fēng)蝕量(kg m-2a-1)；Qrmax表示潛在風(fēng)力的最大輸沙能力(kg/m)；sr表示潛在關(guān)鍵地塊長度(m)；G為防風(fēng)固沙的物質(zhì)量(kg m-2a-1)；z為所計(jì)算的下風(fēng)向距離(m),本次計(jì)算取50 m；WF表示氣候因子(kg/m)；EF為土壤可蝕性成分(無量綱)；SCF表示土壤結(jié)皮因子(無量綱)；K′表示地表糙度因子(無量綱)；C表示植被因子(無量綱)。

(1)氣候因子WF

氣候因子主要考慮了降水、氣溫、風(fēng)速、雪蓋等自然因素影響下風(fēng)力對土壤的搬運(yùn)能力,其中,風(fēng)力因子為氣候因子中風(fēng)速數(shù)據(jù)的相關(guān)計(jì)算,主要計(jì)算公式如下：

WF=Wf×(ρ/g)×SW×SD

(8)

Wf=u2(u2-u1)2×Nd

(9)

式中,WF為氣候因子(kg/m)；Wf為風(fēng)力因子(m/s3)；g為重力加速度(m/s2)；ρ為空氣密度(kg/m3)；SW為土壤濕度因子；SD為雪蓋因子(無積雪覆蓋天數(shù)/研究總天數(shù)),定義雪蓋深度<25.4 mm時為無積雪覆蓋；u1為2 m處起沙風(fēng)速,取5 m/s；u2為由10 m轉(zhuǎn)換為2 m處的月均風(fēng)速值(m/s)；Nd為各月風(fēng)速大于5 m/s的天數(shù)。

(2)土壤可蝕性因子EF、土壤結(jié)皮因子SCF

土壤可蝕性因子與土壤結(jié)皮因子分別表示在一定理化條件下土壤受風(fēng)蝕影響的程度與土壤結(jié)皮抵抗風(fēng)蝕能力的大小,二者的表達(dá)式分別為：

(10)

(11)

土壤質(zhì)地?cái)?shù)據(jù)為國際標(biāo)準(zhǔn),采用對數(shù)正態(tài)分布模型將土壤粒徑轉(zhuǎn)換為RWEQ模型需要的美國標(biāo)準(zhǔn)。其中,sa為土壤粗砂含量(%)；si為土壤粉砂含量(%)；cl為土壤粘粒含量(%)；OM為土壤有機(jī)質(zhì)含量(%)；CaCO3為碳酸鈣含量(%),本研究中取0。

(3)植被因子C

植被因子表示一定的植被條件對風(fēng)蝕的抑制程度。

C=e-0.0483(SC×100)

(12)

SC=(NDVI-NDVImin)/(NDVImax-NDVImin)

(13)

式中,SC為植被覆蓋度(%),NDVI、NDVImax、NDVImin分別為NDVI的實(shí)際值、最大值和最小值。

(4)地表粗糙度因子K′

地表糙度因子用于模擬地表粗糙程度對風(fēng)蝕的影響因子。

K′=cosα

(14)

式中,α為坡度,由30 m的DEM數(shù)據(jù)經(jīng)過ArcGIS的slope模塊計(jì)算得到。

2.2 防風(fēng)固沙服務(wù)保有率模擬

RWEQ模型計(jì)算的防風(fēng)固沙量受降水、風(fēng)場強(qiáng)度等氣候變化的影響較大,不能準(zhǔn)確地評估生態(tài)系統(tǒng)自身防風(fēng)固沙的效果。通過計(jì)算防風(fēng)固沙服務(wù)保有率能夠避免氣候因素對服務(wù)功能的影響,進(jìn)一步分析防風(fēng)固沙服務(wù)能力。防風(fēng)固沙服務(wù)功能保有率計(jì)算公式為：

(15)

式中,F表示防風(fēng)固沙服務(wù)功能保有率；G為防風(fēng)固沙量(kg m-2a-1)；SR為潛在風(fēng)蝕量(kg m-2a-1)。

2.3 防風(fēng)固沙服務(wù)變化趨勢分析

本研究采用最小二乘法線性擬合函數(shù)分析渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)物質(zhì)量1980—2018年的變化趨勢及其空間分布變化。據(jù)1980—2018年重大生態(tài)工程與生態(tài)建設(shè)實(shí)施的年份為關(guān)鍵節(jié)點(diǎn),在柵格尺度上分析1980—2018年、1980—2000年、2000—2010年以及2010—2018年四個時間段渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)物質(zhì)量變化趨勢及空間分布特征。

(16)

式中, Trend為防風(fēng)固沙服務(wù)物質(zhì)量或保有率的變化趨勢；n為研究總年數(shù)；vi為第i年防風(fēng)固沙服務(wù)物質(zhì)量年均值。 Trend>0,表明防風(fēng)固沙能力呈增加趨勢； Trend<0,表明防風(fēng)固沙能力呈下降趨勢。

2.4 防風(fēng)固沙服務(wù)驅(qū)動因素分析

2.4.1地理探測器

地理探測器主要包括4個探測工具：因子探測、風(fēng)險(xiǎn)區(qū)探測、生態(tài)探測、交互作用探測[30]。

(1)因子探測器通過層內(nèi)方差之和(Within Sum of Squares, SSW)與全區(qū)總方差(Total Sum of Squares, SSW)計(jì)算各探測因子對防風(fēng)固沙服務(wù)物質(zhì)量空間分異的解釋力,用q值度量,計(jì)算公式如下：

(17)

(2)風(fēng)險(xiǎn)區(qū)探測用于判斷兩個影響因子的子區(qū)域間的防風(fēng)固沙物質(zhì)量均值是否有顯著差別,并通過t統(tǒng)計(jì)量來檢驗(yàn)：

(18)

(3)生態(tài)探測通過統(tǒng)計(jì)量F評估兩個影響因子X1與X2對防風(fēng)固沙物質(zhì)量空間分異的影響是否顯著：

(19)

式中,NX1與NX2分別為影響因子X1與X2的樣本數(shù)目,SSWX1和SSWX2分別表示對應(yīng)影響因子分層的層內(nèi)方差之和；L1與L2表示影響因子X1與X2的分層數(shù)目。

(4)交互探測能夠判斷不同影響因子的共同作用對防風(fēng)固沙服務(wù)物質(zhì)量的解釋力的整體影響(增加、減弱或相互獨(dú)立)。通過比較影響因子X1與X2對防風(fēng)固沙物質(zhì)量各自作用與交互作用的q值(q(X1)、q(X2)、q(X1∩X2))大小,將交互作用分為5類(表 1)。

2.4.2基于地理探測器的驅(qū)動因素分析

防風(fēng)固沙服務(wù)物質(zhì)量的變化同時受到氣候、土壤、植被類型等自然因素和退牧、禁牧等人類活動的影響。基于RWEQ模型以及防風(fēng)固沙服務(wù)或其他服務(wù)驅(qū)動因素研究中所涉及的驅(qū)動因素的出現(xiàn)頻率和研究結(jié)果[23, 26, 29, 39],本研究選取地貌類型、年降水量、土壤類型等11類自然因素和單位面積國內(nèi)生產(chǎn)總值(Gross Domestic Product, GDP)、人口密度等9類社會經(jīng)濟(jì)因素作為探測因子,分析渾善達(dá)克地區(qū)2000、2005、2010、2015、2018年和1980—2018年年均防風(fēng)固沙服務(wù)物質(zhì)量變化的驅(qū)動因素(表 2)?；跍喩七_(dá)克當(dāng)?shù)氐膶?shí)際情況以及數(shù)據(jù)的不同類型,本研究將渾善達(dá)克的地貌類型、土壤類型、植被類型分別分為4類(平原、臺地、丘陵、山地)、7類(針葉林、闊葉林、灌叢、荒漠、草原、草叢和草甸、栽培植被)、7類(淋溶土、半淋溶土、鈣層土、干旱土、初育土、半水成土、水成土、鹽堿土),數(shù)值類變量利用ArcGIS中的自然間斷法,將人口密度、單位面積GDP分為4類,其余變量分為10類,分析2000—2018年渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量時空變化的驅(qū)動因素。

表1 地理探測器交互作用類型

表2 探測因子指標(biāo)

3 結(jié)果分析

3.1 防風(fēng)固沙服務(wù)時空動態(tài)變化

3.1.1防風(fēng)固沙量時空變化

根據(jù)渾善達(dá)克地區(qū)各年單位面積防風(fēng)固沙量平均值(表 3),可以發(fā)現(xiàn),1980—2018年,渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量的變化范圍為13.01—22.90 kg m-2a-1,防風(fēng)固沙總量為2.18×1012—3.81×1012kg/a。1980年,渾善達(dá)克全區(qū)單位面積防風(fēng)固沙量和總量均為最高,隨后二者在一定范圍內(nèi)上下波動,其中,單位面積防風(fēng)固沙量在13.01—19.47 kg m-2a-1之間變化,在2018年小幅度回升。相比于防風(fēng)固沙總量,渾善達(dá)克地區(qū)潛在風(fēng)蝕量和實(shí)際風(fēng)蝕量有較大幅度的變化。1980年潛在風(fēng)蝕量最高,1990年有較大幅度的下降,在1995年保持穩(wěn)定,而后在2000年又有較大幅度下降,之后一直保持穩(wěn)定；1980—1995年實(shí)際風(fēng)蝕量基本保持穩(wěn)定,在2000年有顯著下降,之后保持穩(wěn)定。

表3 1980—2018年渾善達(dá)克防風(fēng)固沙服務(wù)及保有率年均值

空間上,渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)總體表現(xiàn)為西北部荒漠、中部沙地以及南部林地普遍高于其他地區(qū)的空間分布格局(圖2)。其中,正藍(lán)旗、正鑲白旗、蘇尼特左旗的單位面積防風(fēng)固沙量較高,數(shù)值均在20 kg m-2a-1以上；張北縣、沽源縣、康?？h、尚義縣單位防風(fēng)固沙量較低,主要在4—10 kg m-2a-1。此外,單位面積防風(fēng)固沙量的高值區(qū)逐漸從中部偏東地區(qū)向西北地區(qū)移動。1980—2018年,渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)高值分布面積整體呈現(xiàn)波動減小的趨勢。1980年,高值區(qū)域(≥30 kg/m2)面積最大,為4.63×104km2；2005年,面積縮減到最小,為2.14×104km2。此外,2015年單位面積防風(fēng)固沙量的極高值區(qū)域(≥40 kg/m2)面積縮減到最小,僅為0.63×104km2；而2018年區(qū)域面積又有所回彈,高值區(qū)域面積增長至4.41 km2,其中,極高值區(qū)域面積為1.85×104km2。

圖2 1980—2018年渾善達(dá)克防風(fēng)固沙量空間分布格局Fig.2 The spatial pattern of wind erosion prevented amount of the Otindag during 1980—2018

綜合1980—2018年渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量的結(jié)果,可以發(fā)現(xiàn),38年來渾善達(dá)克地區(qū)的風(fēng)蝕程度有所好轉(zhuǎn),潛在風(fēng)蝕與實(shí)際風(fēng)蝕情況均有不同程度的下降,防風(fēng)固沙總量也有小幅下降(表 3)。這可能與內(nèi)蒙古地區(qū)風(fēng)蝕發(fā)生的原動力(即風(fēng)速)的變化,以及影響植被固沙作用的氣候條件的不穩(wěn)定有關(guān),全球氣候變暖和內(nèi)蒙古地區(qū)近30年來風(fēng)速的減小緩解了渾善達(dá)克地區(qū)風(fēng)蝕情況[6, 19,26]。本研究結(jié)果還顯示,2000年以來潛在風(fēng)蝕量、實(shí)際風(fēng)蝕量相比以前有顯著下降并保持穩(wěn)定,這與2000年以來退耕還林還草和京津風(fēng)沙源治理工程的實(shí)施有較大關(guān)系[35]。

3.1.2防風(fēng)固沙保有率時空變化

1980—2018年,渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)保有率自1980年的最低值70.79%波動上升,至2018年達(dá)到最大值94.28%(表 3)。其中,1980—1995年,保有率變化幅度較小,保持在70%左右,2000年達(dá)到89.11%后,出現(xiàn)一定程度的下降,2005—2018年,防風(fēng)固沙保有率以83.31%為起點(diǎn),大致呈現(xiàn)較為穩(wěn)定的上升態(tài)勢。

空間上,渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)保有率整體呈現(xiàn)出自東南向西北遞減的空間分布特征。東南部的豐寧滿族自治縣、圍場滿族蒙古族自治縣和克什克騰旗保有率始終高于其他地區(qū),西北部蘇尼特左旗、蘇尼特右旗的荒漠及中覆蓋度草地地區(qū)的保有率則相對較低(圖3)。1980—2018年,防風(fēng)固沙服務(wù)保有率低值區(qū)(保有率<60%)逐漸向蘇尼特右旗北部與中西部沙地邊緣地區(qū)集中,且面積顯著減小,其中,極低值區(qū)域(保有率<40%)面積自1995年開始顯著減少。1980年,低值區(qū)面積最大,以極低值區(qū)域?yàn)橹?低值區(qū)面積為4.83×104km2,包括極低值面積3.07×104km2；2000年,防風(fēng)固沙服務(wù)保有率低值面積顯著減小,而2005年,面積又?jǐn)U大至3.06×104km2。值得注意的是,極低值的面積沒有出現(xiàn)大幅度回升,且2005年之后的極低值面積變化都不大。2018年,全區(qū)保有率低值區(qū)域面積最小,為0.19×104km2,其中,極低值區(qū)域面積僅為0.02×104km2。

圖3 1980—2018年渾善達(dá)克防風(fēng)固沙保有率空間分布格局Fig.3 The spatial pattern of soil retention rate of the Otindag during 1980—2018

整體來看,渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)保有率整體顯著提高,年平均保有率為82.42%。尤其是2000年之后極低值區(qū)域面積的顯著減少可能與2000年之后實(shí)行的京津風(fēng)沙源治理工程、退耕還林還草等生態(tài)恢復(fù)工程相關(guān),說明植被對風(fēng)沙的抑制作用逐漸顯著，生態(tài)恢復(fù)取得了初步的效果；而2005年,保有率的短暫性降低,可能與當(dāng)年氣候相對干旱,阻礙了植被的生長,削弱了防風(fēng)固沙服務(wù)功能有關(guān)[24, 40]。

3.2 防風(fēng)固沙服務(wù)趨勢變化

3.2.1單位面積防風(fēng)固沙服務(wù)變化趨勢

總體而言,1980—2018年渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量以0.43 kg m-2a-1的速度減少。其中,1980—2000年、2000—2010年,單位面積防風(fēng)固沙量分別以1.38 kg m-2a-1和1.04 kg m-2a-1的速度減少,2010—2018年,單位面積防風(fēng)固沙量則以1.77 kg m-2a-1的速度增加(表4)。

表4 不同時段渾善達(dá)克防風(fēng)固沙量以及保有率年均變化趨勢

空間格局的變化趨勢上來看,1980—2018年,渾善達(dá)克大部分地區(qū)單位面積防風(fēng)固沙量表現(xiàn)為下降趨勢,增加與顯著增加區(qū)域集中在蘇尼特左旗、蘇尼特右旗,以及正鑲白旗、正藍(lán)旗和阿巴嘎旗的部分區(qū)域(圖4)。

圖4 1980—2018年各時段渾善達(dá)克單位面積防風(fēng)固沙量年均空間變化趨勢Fig.4 Spatial change trend of average annual wind erosion prevented amount of the Otindag in different period between 1980 and 2018

1980—2000年單位面積防風(fēng)固沙量變化的整體空間格局與1980—2018年的整體格局大致相同(圖4),說明2000年以來渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量相對穩(wěn)定。2000—2010年除蘇尼特左旗、渾善達(dá)克沙地及周邊地區(qū)為降低或顯著降低趨勢外,大部分地區(qū)單位面積防風(fēng)固沙量為增加或顯著增加趨勢(圖4)。2010—2018年單位面積防風(fēng)固沙量整體恢復(fù)下降的趨勢,增加趨勢集中區(qū)轉(zhuǎn)移至蘇尼特左旗北部、阿巴嘎旗以及渾善達(dá)克沙地地區(qū)(圖4)。在所有時段內(nèi),蘇尼特右旗、蘇尼特左旗以及渾善達(dá)克沙地地區(qū)單位面積防風(fēng)固沙量均出現(xiàn)明顯的增加趨勢,防風(fēng)固沙能力得到了穩(wěn)定且連續(xù)的提高。

圖5 1980—2018年年均降水量與年均風(fēng)速變化Fig.5 Changes in annual average precipitation and wind speed from 1980 to 2018

結(jié)合表 3可以看出,1980—2000年潛在風(fēng)蝕與實(shí)際風(fēng)蝕均呈波動減小趨勢,2000—2010年二者則基本保持穩(wěn)定。兩個時間段單位面積防風(fēng)固沙量變化趨勢呈現(xiàn)明顯差異(圖4),可能由于1980—2000年風(fēng)速顯著下降,降水年際變化幅度大,風(fēng)蝕動力減弱使得大部分地區(qū)單位面積防風(fēng)固沙量為下降趨勢；2000年以后,風(fēng)速基本保持穩(wěn)定,降水增多促進(jìn)了植被的自然恢復(fù),使得大部分地區(qū)在2000—2010年單位面積防風(fēng)固沙量為增加趨勢(圖5)。

3.2.2防風(fēng)固沙保有率變化趨勢

1980—2018年渾善達(dá)克地區(qū)防風(fēng)固沙保有率以3.46%/a的速度增加(表 4)。其中,1980—2000年增速較快,為6.29%/a；2000—2010年轉(zhuǎn)變?yōu)闇p小趨勢且速度相對較慢,為-1.34%/a；2010—2018年保有率恢復(fù)增加趨勢,為2.39%/a。

1980—2018年渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)保有率總體呈現(xiàn)增加趨勢,其中,西北荒漠地區(qū)及其與中部沙地過渡地區(qū)增加趨勢顯著(圖6)。以2000年、2010年為時間節(jié)點(diǎn),發(fā)現(xiàn)不同時段內(nèi),防風(fēng)固沙保有率變化趨勢明顯不同：1980—2000年渾善達(dá)克絕大部分地區(qū)保有率為增加趨勢,其中,蘇尼特右旗、蘇尼特左旗以及正鑲白旗為顯著增加趨勢；2000—2010年以蘇尼特右旗、蘇尼特左旗、鑲黃旗、正鑲白旗為主的北部、中部地區(qū)呈現(xiàn)降低趨勢,說明該時段內(nèi)植被對防風(fēng)固沙服務(wù)的影響減弱；2010—2018年渾善達(dá)克防風(fēng)固沙服務(wù)保有率又整體恢復(fù)增加趨勢,但顯著增加區(qū)域僅出現(xiàn)在蘇尼特右旗中北部的小部分地區(qū)(圖6)。

圖6 1980—2018年各時段渾善達(dá)克防風(fēng)固沙服務(wù)保有率年均空間變化趨勢Fig.6 Spatial change trend of average annual soil retention rate of the Otindag in different period between 1980 and 2018

由此可見,1980—2018年渾善達(dá)克地區(qū)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力有了較為顯著的提高,在2010年之后,顯著增加區(qū)域的大幅減少說明渾善達(dá)克地區(qū)生態(tài)系統(tǒng)防風(fēng)固沙能力已經(jīng)維持在一個相對穩(wěn)定的狀態(tài)。而2000—2010年,由于風(fēng)速持續(xù)增加導(dǎo)致潛在風(fēng)蝕量增加,在單位面積防風(fēng)固沙量相對穩(wěn)定的情況下,該時段內(nèi)防風(fēng)固沙服務(wù)保有率出現(xiàn)短暫下降。

3.3 防風(fēng)固沙服務(wù)驅(qū)動因素分析

3.3.1單因素作用分析

2000—2018年多年平均的單因素作用結(jié)果顯示,對防風(fēng)固沙量空間變化的影響程度(即解釋力q值)較大的探測因子，其影響程度由大到小依次為：土壤類型(75.15%)、年末牲畜數(shù)量(25.20%)、年降水量(24.16%)、人工造林面積(21.11%)、年糧食產(chǎn)量(18.96%)、年糧食播種面積(18.04%)、NDVI(16.80%)、植被覆蓋度(16.74%)、年肉類產(chǎn)量(16.12%)、年均溫(13.98%)、平均風(fēng)速(11.10%)、高程(10.41%)(圖7)。綜合2000、2005、2010、2015、2018年的探測結(jié)果,可以發(fā)現(xiàn),土壤類型始終是防風(fēng)固沙服務(wù)空間變化的主要驅(qū)動因素,且影響程度穩(wěn)定在70%左右。同樣影響較強(qiáng)且持續(xù)時間較久的因子還包括年降水量和平均風(fēng)速,解釋力分別在24.16%和11.10%左右。

由于自然環(huán)境與社會經(jīng)濟(jì)條件不同,各年份防風(fēng)固沙量的主要影響因子也存在差異：2000年退耕還林還草、京津風(fēng)沙源治理等生態(tài)工程的實(shí)施,年糧食播種面積、年末牲畜數(shù)量以及人工造林面積對防風(fēng)固沙服務(wù)空間變化的影響較強(qiáng)。相比于2000年,2005—2015年年降水量與平均風(fēng)速對防風(fēng)固沙服務(wù)空間格局變化的影響程度有較為明顯的增加,年降水量的影響相對穩(wěn)定,解釋力保持在25%左右；平均風(fēng)速的影響程度則逐年增強(qiáng)。這可能與該時間段內(nèi)年降水量相對穩(wěn)定,風(fēng)速的持續(xù)增強(qiáng)有關(guān)。2000—2018年,由于退牧還草、劃區(qū)輪牧、草地補(bǔ)播等生態(tài)工程的實(shí)施,渾善達(dá)克地區(qū)年末牲畜數(shù)量明顯減少并在2010年后保持穩(wěn)定,生態(tài)系統(tǒng)植被也有所恢復(fù)。因此,年末牲畜數(shù)量對防風(fēng)固沙服務(wù)空間變化的影響程度逐年減弱,而NDVI和植被覆蓋度的影響程度則有所增加。

如圖8所示,單位面積防風(fēng)固沙量均值隨不同探測因子的分類數(shù)值的增加而發(fā)生不同的變化(探測因子分類間斷值見表 5)：單位面積防風(fēng)固沙量均值隨NDVI、植被覆蓋度、坡度、起伏度、人口密度、單位面積GDP的增加而波動下降。其中,植被覆蓋度、NDVI在第3分類處對應(yīng)的單位面積防風(fēng)固沙量均值最大,隨后為明顯的波動下降；而人口密度與單位面積GDP的持續(xù)增加,使單位面積防風(fēng)固沙量均值呈現(xiàn)不同幅度的持續(xù)下降。

圖7 各年探測因子q值變化Fig.7 q value of detecting factors in different yearX1：高程 Digital Elevation Model (DEM)；X2：地貌類型 Geomorphic type；X3：坡度 Slope；X4：起伏度 Prominence；X5：年降水量 Annual precipitation；X6：年均溫 Average annual temperature；X7：平均風(fēng)速 Average annual wind speed；X8：植被類型 Vegetation type；X9：土壤類型 Soil type；X10：歸一化植被指數(shù) Normalized Difference Vegetation Index (NDVI)；X11：植被覆蓋度 Fractional Vegetation Cover (FVC)；X12：單位面積國內(nèi)生產(chǎn)總值 Gross Domestic Product (GDP) per unit area；X13：人口密度 Population density,X14：人工造林面積 Artificial afforestation area；X15：與道路的距離 Distance from the road；X16：與城鎮(zhèn)的距離 Distance from the town；X17：年糧食播種面積 Annual total grain sown area；X18：年糧食產(chǎn)量 Annual total grain production；X19：年肉類產(chǎn)量 Annual total meat production；X20：年末牲畜數(shù)量 Number of livestock at the end of the year

圖8 各探測因子單位面積防風(fēng)固沙量均值變化Fig.8 Change of average wind erosion prevented of different detecting factors各探測因子分類具體間斷值見表5

表5 各探測因子的分類間斷值

單位面積防風(fēng)固沙量均值隨平均風(fēng)速、與城鎮(zhèn)的距離、與道路的距離的增加表現(xiàn)為波動上升的趨勢。隨著人工造林面積、年糧食產(chǎn)量、年糧食播種面積以及年肉類產(chǎn)量增加,單位面積防風(fēng)固沙量均值變化較為一致,均為先波動增加,在某一分類處達(dá)到最大之后,出現(xiàn)顯著下降,隨后保持相對穩(wěn)定的上下波動。

因此,在17類數(shù)值型驅(qū)動因素中,單位面積防風(fēng)固沙量均值與平均風(fēng)速、與城鎮(zhèn)的距離、與道路的距離整體表現(xiàn)為正相關(guān)關(guān)系；與坡度、起伏度、人口密度、單位面積GDP大致表現(xiàn)為負(fù)相關(guān)關(guān)系；與植被覆蓋度和NDVI表面為先緩慢增加,隨后顯著降低的關(guān)系。

在2000—2018年多年平均生態(tài)探測的結(jié)果中,年降水量、土壤類型、人口密度與其他所有影響因子之間均存在顯著差異,年均溫、NDVI、植被覆蓋度、人工造林面積和肉類產(chǎn)量也與其他大部分影響因子之間存在顯著性差異,(表 6)。結(jié)合因子探測的結(jié)果,可以判斷土壤類型、年降水量、人工造林面積是影響渾善達(dá)克地區(qū)單位面積防風(fēng)固沙量空間變化的兩個最主要因素。

表6 探測因子之間統(tǒng)計(jì)顯著性

3.3.2多因素交互作用分析

根據(jù)多年平均交互探測器結(jié)果(表 7),發(fā)現(xiàn)各因子之間對渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)的影響均存在增強(qiáng)型交互作用。除地貌類型、年均溫、植被類型、與城鎮(zhèn)的距離與其他大部分因子之間為非線性增強(qiáng)關(guān)系外,大部分因子與其他因子之間普遍存在雙因子增強(qiáng)關(guān)系。土壤類型、年降水量、人工造林面積、牲畜數(shù)量與其他因子之間的雙因子增強(qiáng)作用對防風(fēng)固沙服務(wù)空間變化的影響程度較大。尤其是土壤類型與人工造林面積、糧食產(chǎn)量、糧食播種面積之間的雙因子增強(qiáng)作用最強(qiáng),對防風(fēng)固沙量變化的影響程度(q值)均達(dá)到最大值82%。這說明土壤類型也加強(qiáng)了其他因素對防風(fēng)固沙服務(wù)空間分布的影響程度。

而年均溫與除地貌與土壤類型外的其他因子之間的非線性增強(qiáng)關(guān)系,在所有因子之間的非線性增強(qiáng)作用中對防風(fēng)固沙服務(wù)空間變化的影響程度最高,在30%左右。可以看出,盡管年均溫對防風(fēng)固沙服務(wù)空間變化的單因素作用的影響較弱,但通過改變區(qū)域自然環(huán)境狀態(tài)或人類活動,使得兩因子間的交互作用對防風(fēng)固沙的影響水平高于二者單獨(dú)作用之和,從而間接影響了防風(fēng)固沙服務(wù)的空間變化。因此,年均溫對防風(fēng)固沙量的變化也存在較為重要的間接影響。

表7 各因子之間交互作用類型與q值

4 結(jié)論

本研究采用修正風(fēng)蝕模型模擬了渾善達(dá)克地區(qū)1980—2018年防風(fēng)固沙服務(wù)及保有率的時空變化,以重大生態(tài)建設(shè)實(shí)施等年份為時間節(jié)點(diǎn),分析不同時段內(nèi)二者的變化趨勢,并通過地理探測器探究年降水量、人口密度等20項(xiàng)自然與人為因子及因子之間的交互作用對該地區(qū)防風(fēng)固沙服務(wù)空間格局變化的影響。結(jié)果表明：

(1)1980—2018年,渾善達(dá)克地區(qū)風(fēng)速的下降,劃區(qū)輪牧、草地補(bǔ)播等生態(tài)恢復(fù)政策的實(shí)施,加之適宜植被生長的氣候變化,使得防風(fēng)固沙服務(wù)能力整體得到一定的提高,單位面積防風(fēng)固沙量波動減少,防風(fēng)固沙服務(wù)空間分布差距逐漸縮小；

(2)土壤類型、年末牲畜數(shù)量、年降水量、人工造林面積是造成渾善達(dá)克生態(tài)功能區(qū)防風(fēng)固沙服務(wù)空間差異及分布變化的主要影響因素。分別在初育土、52.84—86.86萬頭、404.21 mm、5167 hm2時,各自對應(yīng)的單位面積防風(fēng)固沙量均值最大。

(3)各驅(qū)動因素間的交互作用都會放大單因子對渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)空間分布的影響。其中,土壤類型與其他因子的交互作用對防風(fēng)固沙服務(wù)空間變化的影響,均大于單因子分別對防風(fēng)固沙服務(wù)空間變化的影響,整體對防風(fēng)固沙服務(wù)空間變化影響程度最大。

(4)年均溫雖然對單位面積防風(fēng)固沙量空間分布的影響較小,但當(dāng)它與除地貌類型與土壤類型的其他因子共同作用時,對防風(fēng)固沙服務(wù)空間變化影響均高于兩個驅(qū)動因素各自對防風(fēng)固沙服務(wù)的影響之和,對單位面積防風(fēng)固沙量空間分布具有較強(qiáng)的間接影響作用。

5 討論

在長時間序列下,渾善達(dá)克防風(fēng)固沙重點(diǎn)生態(tài)功能區(qū)的防風(fēng)固沙能力得到了一定程度的提高,風(fēng)蝕程度有所緩解,單位面積防風(fēng)固沙量高值區(qū)從中部、東部沙地向西北荒漠地區(qū)轉(zhuǎn)移。特別是在2000年之后,防風(fēng)固沙能力與單位面積防風(fēng)固沙量均穩(wěn)定在相對較低的區(qū)間范圍內(nèi)。結(jié)合地理探測器的分析結(jié)果,可以得出,1980—2018年渾善達(dá)克生態(tài)功能區(qū)的上述變化,主要與該段時間內(nèi)土壤類型,年均風(fēng)速、年降水量、人工造林面積以及牲畜數(shù)量的時空格局變化有關(guān)。1980—2000年,年均風(fēng)速、年降水量等因素年際變化不穩(wěn)定,京津風(fēng)沙源治理工程尚未開始,是導(dǎo)致渾善達(dá)克地區(qū)防風(fēng)固沙服務(wù)能力整體較弱,風(fēng)蝕程度較高的主要原因。2000年之后,圍欄豐育、劃區(qū)輪牧、草地補(bǔ)播等生態(tài)保護(hù)政策的實(shí)施,使得人工造林面積、牲畜數(shù)量在時空分布上較為穩(wěn)定,加之年均風(fēng)速、年降水量等氣候條件促進(jìn)植被恢復(fù),使得渾善達(dá)克生態(tài)功能區(qū)防風(fēng)固沙服務(wù)能力整體提高,區(qū)域間防風(fēng)固沙服務(wù)能力的差異也逐漸縮小。此外,在各驅(qū)動因素的交互作用中,土壤類型的時空分布決定了區(qū)域內(nèi)地理?xiàng)l件的差異,從而放大其他因子對防風(fēng)固沙服務(wù)空間分布的影響,使二者對防風(fēng)固沙服務(wù)分布的影響大于單因子影響但小于二者各自的影響之和。年均溫能夠影響除地貌類型和土壤類型外其他因子的時空分布,當(dāng)年均溫與其他因子共同作用于防風(fēng)固沙服務(wù)空間分布時,其影響程度高于兩個單因子獨(dú)立的影響程度之和,具有較強(qiáng)的間接作用。

但是,兩種以上因子之間的交互作用對防風(fēng)固沙服務(wù)能力及時空變化的影響如何,周邊地區(qū)對渾善達(dá)克生態(tài)功能區(qū)防風(fēng)固沙服務(wù)是否產(chǎn)生影響,還需要進(jìn)一步地深入研究,為渾善達(dá)克地區(qū)生態(tài)建設(shè)與生態(tài)恢復(fù)提供科學(xué)依據(jù)。