米朝娟,周自翔,劉 婷,武 佳,孫彥旭

西安科技大學(xué)測(cè)繪科學(xué)與技術(shù)學(xué)院,西安 710054

生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)是干旱和半干旱地區(qū)生態(tài)系統(tǒng)最重要的生態(tài)服務(wù)功能之一,實(shí)施防風(fēng)固沙工程是抑制和固定風(fēng)沙區(qū)沙塵不被再次移動(dòng)以實(shí)現(xiàn)荒漠化逆轉(zhuǎn),有利于區(qū)域經(jīng)濟(jì)和人類(lèi)福祉的可持續(xù)發(fā)展[1]。生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)通常從風(fēng)蝕角度為防風(fēng)固沙服務(wù)提供研究。學(xué)術(shù)界對(duì)風(fēng)蝕的評(píng)估始于1940年代,風(fēng)蝕模型是評(píng)估土壤風(fēng)蝕狀況的主要技術(shù)手段,包括德克薩斯模型(Texas erosion analysis model,TEAM)[2—3]、風(fēng)蝕預(yù)報(bào)系統(tǒng)(WEPS)[4—5]、風(fēng)蝕方程(WEQ)[6]、修正風(fēng)蝕方程模型(Revised Wind Erosion Equation,RWEQ)[7—9]等。由于一些模型主要集中在風(fēng)蝕發(fā)生的可能性,無(wú)法評(píng)估潛在風(fēng)蝕量。而大多數(shù)學(xué)者考慮到REWQ模型具有良好的適用性,涉及氣候條件、植被覆蓋、土壤易蝕性、土壤結(jié)皮、地表粗糙度等綜合因素來(lái)計(jì)算植被造成的沙粒滯留量所發(fā)生的風(fēng)力侵蝕[10—12]。部分學(xué)者不斷探索和驗(yàn)證該模型并通過(guò)調(diào)整參數(shù),利用RWEQ模型對(duì)寧夏[7—12]、內(nèi)蒙古[13—14]、京津冀[15—17]等中國(guó)干旱半干旱風(fēng)沙區(qū)的風(fēng)蝕狀況進(jìn)行了評(píng)估,取得了良好的效果。

植被是影響防風(fēng)固沙生態(tài)功能的關(guān)鍵指標(biāo),也是檢驗(yàn)防風(fēng)固沙區(qū)生態(tài)保護(hù)成效的重要依據(jù)[11]。恢復(fù)植被有利于保持土壤水分、改變地表狀態(tài),從而降低風(fēng)蝕發(fā)生的潛力[12]。風(fēng)蝕的動(dòng)態(tài)變化是自然氣候變化和人類(lèi)活動(dòng)共同作用的結(jié)果,退耕還林使地表覆蓋發(fā)生顯著變化,進(jìn)而與生態(tài)系統(tǒng)服務(wù)能力的變化相關(guān)度最高[18]。應(yīng)用地理探測(cè)器模型分析自然和人為因子對(duì)寧夏地區(qū)防風(fēng)固沙服務(wù)功能空間格局形成及演變過(guò)程中的貢獻(xiàn)與交互作用[1]發(fā)現(xiàn)土壤類(lèi)型和植被類(lèi)型在防風(fēng)固沙服務(wù)功能空間分布格局中發(fā)揮最重要的作用。而以歸一化植被指數(shù)(NDVI)計(jì)算植被覆蓋度,結(jié)合氣象要素從不同尺度來(lái)評(píng)估寧夏靈武白芨灘自然保護(hù)區(qū)防風(fēng)固沙功能時(shí)[12],植被覆蓋度的增加與研究區(qū)的防風(fēng)固沙服務(wù)成正比。朱趁趁等[14]利用RWEQ模型實(shí)現(xiàn)土地利用變化對(duì)內(nèi)蒙古荒漠草原的防風(fēng)固沙能力進(jìn)行評(píng)估。結(jié)果表明：土地利用方式以林地恢復(fù)、建設(shè)用地?cái)U(kuò)張、不同覆蓋度草地轉(zhuǎn)化變化為主對(duì)防風(fēng)固沙服務(wù)有一定的增強(qiáng)作用。徐潔等[19]基于RWEQ和HYSPLIT模型模擬了生態(tài)功能區(qū)防風(fēng)固沙服務(wù)和防風(fēng)固沙服務(wù)空間流動(dòng)路徑,從生態(tài)系統(tǒng)服務(wù)流動(dòng)的角度建立了生態(tài)功能區(qū)及其防風(fēng)固沙服務(wù)受益區(qū)之間的時(shí)空聯(lián)系。而以往的研究大多是大區(qū)域尺度,很少評(píng)估土壤風(fēng)蝕對(duì)風(fēng)沙區(qū)邊緣地帶的影響,而且采用的氣象數(shù)據(jù)精度不夠。由于目前風(fēng)沙邊緣地區(qū)的氣象環(huán)境變化、植被恢復(fù)與防風(fēng)固沙能力之間的相互作用機(jī)制與如何保證防風(fēng)固沙服務(wù)功能的可持續(xù)性還需進(jìn)一步研究。

基于生態(tài)工程的防風(fēng)固沙、植被恢復(fù)、土壤保育等生態(tài)功能監(jiān)測(cè)與評(píng)價(jià),可為促進(jìn)地區(qū)經(jīng)濟(jì)和社會(huì)的可持續(xù)發(fā)展提供實(shí)踐指導(dǎo)[20]。在我國(guó)土地沙化依然在擴(kuò)展的趨勢(shì)下,陜北長(zhǎng)城沿線(xiàn)處于毛烏素沙漠的邊緣地區(qū),長(zhǎng)期受到自然因素和人類(lèi)活動(dòng)擾動(dòng)影響,是最脆弱的生態(tài)系統(tǒng)之一[21]。目前,榆林大地基本上消滅了沙漠化土地現(xiàn)象,在全球沙漠治理的意義上率先實(shí)現(xiàn)了荒漠化的逆轉(zhuǎn)。但該地在干旱少雨、生態(tài)脆弱的自然條件下,導(dǎo)致頻發(fā)干旱和風(fēng)蝕現(xiàn)象。因此,從生態(tài)系統(tǒng)服務(wù)的角度探討植被恢復(fù)、氣象環(huán)境變化與防風(fēng)固沙服務(wù)能力之間相互作用機(jī)質(zhì)以及如何保證防風(fēng)固沙服務(wù)的可持續(xù)性。本文研究?jī)?nèi)容包括：(1)分析陜北長(zhǎng)城沿線(xiàn)2000—2018年植被覆蓋時(shí)空變化。(2)利用RWEQ模型模擬長(zhǎng)城沿線(xiàn)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)物理量的時(shí)空模式,實(shí)現(xiàn)防風(fēng)固沙服務(wù)的量化。(3)利用情景分析評(píng)估風(fēng)蝕因子改變的情況下,植被恢復(fù)對(duì)防風(fēng)固沙服務(wù)物理量和服務(wù)能力的影響。擬解決的關(guān)鍵問(wèn)題：(1)以定量化研究的結(jié)果為長(zhǎng)城沿線(xiàn)2000—2018年的風(fēng)蝕現(xiàn)象的變化情況提供科學(xué)依據(jù)。(2)探討風(fēng)蝕因子改變,植被恢復(fù)對(duì)風(fēng)沙邊緣地區(qū)風(fēng)蝕量的影響及生態(tài)系統(tǒng)防風(fēng)固沙物質(zhì)量和能力的影響。明確氣象環(huán)境改變時(shí),植被恢復(fù)與防風(fēng)固沙服務(wù)的作用機(jī)制,并推進(jìn)土地資源的保護(hù)和建設(shè),實(shí)現(xiàn)可持續(xù)高質(zhì)量發(fā)展，也為我國(guó)陜北地區(qū)防風(fēng)固沙治理工作進(jìn)程提供一定的科學(xué)參考。

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

1.1 研究區(qū)概況

陜北長(zhǎng)城沿線(xiàn)(圖1)位于陜西省最北端,包括定邊、靖邊、橫山、榆陽(yáng)、神木、府谷六個(gè)市縣(區(qū))[21—22],北鄰鄂爾多斯草原,西鄰寧夏自治區(qū)。位于36°45′—39°58′N(xiāo),107°35′—111°29′E之間。作為陜北黃土高原丘陵溝壑區(qū),地處干旱與半干旱地帶,地勢(shì)東低西高、年降水量偏少、土壤類(lèi)型主要以沙土為主,為風(fēng)力侵蝕提供了“有利”的發(fā)生條件。作為干旱頻發(fā)區(qū),面臨著植被退化、水土流失等生態(tài)問(wèn)題,生態(tài)系統(tǒng)脆弱易受人類(lèi)活動(dòng)干擾和氣候變化的影響。

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

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

本文使用的數(shù)據(jù)包括NDVI數(shù)據(jù)、數(shù)字高程模型(Digital Elevation Model,DEM)、土壤數(shù)據(jù)、雪蓋數(shù)據(jù)和氣象數(shù)據(jù)。(1)NDVI數(shù)據(jù)來(lái)源于美國(guó)地質(zhì)勘探局(http://glovis.usgs.gov/)數(shù)據(jù)產(chǎn)品,選用2000—2018年時(shí)間段,空間分辨率為1km,利用MRT(MODIS Reprojection Tools)最大合成年數(shù)據(jù)。(2)DEM數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院資源與環(huán)境數(shù)據(jù)云平臺(tái)(http://www.resdc.cn)空間分辨率為1km。(3)土壤數(shù)據(jù)和積雪覆蓋數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院旱區(qū)寒區(qū)科學(xué)數(shù)據(jù)中心(http://westdc.westgis.ac.cn),土壤數(shù)據(jù)為世界土壤數(shù)據(jù)庫(kù)(Harmonized World Soil Database,HWSD)所提供的中國(guó)土壤數(shù)據(jù)為1∶100萬(wàn)土壤數(shù)據(jù)庫(kù)。中國(guó)雪深長(zhǎng)時(shí)間序列數(shù)據(jù)集提供1979年1月1日到2020年12月31日逐日的中國(guó)范圍的積雪厚度分布數(shù)據(jù),空間分辨率為0.25°。(4)氣象數(shù)據(jù)(風(fēng)速、氣溫、降水)來(lái)源于時(shí)空三級(jí)環(huán)境大數(shù)據(jù)平臺(tái)(https://data.tpdc.ac.cn/zh-hans/)提供的中國(guó)區(qū)域地面要素?cái)?shù)據(jù)集(China meteorological forcing dataset,CMFD)[23—25],該數(shù)據(jù)時(shí)間范圍為1979—2018年,時(shí)間分辨率為3小時(shí),空間分辨率為0.1°,數(shù)據(jù)格式為NETCDF格式,運(yùn)用Python語(yǔ)言提取出研究區(qū)內(nèi)349個(gè)規(guī)則點(diǎn)氣象要素,大于實(shí)際研究區(qū)4個(gè)國(guó)家氣象站點(diǎn),足以從柵格尺度上真實(shí)的反映出研究區(qū)內(nèi)風(fēng)速大于5m/s天數(shù)；(5)土壤濕度數(shù)據(jù)來(lái)源于世界氣候中心(https://www.worldweather.cn/),空間分辨率為0.1°。最終將所有數(shù)據(jù)統(tǒng)一為CGCS_2000_Albers坐標(biāo)系,利用ArcGIS平臺(tái)重采樣為1km分辨率。

2 研究方法

2.1 Theil-Sen斜率估算與M-K顯著性檢驗(yàn)結(jié)合

為了進(jìn)一步體現(xiàn)出研究區(qū)植被變化的顯著性,采用斜率估算和顯著性檢驗(yàn)相結(jié)合的方法。而Theil-Sen斜率估算和Mann-Kendall顯著性檢驗(yàn)是兩種非參數(shù)檢驗(yàn)方法,目前已經(jīng)得到廣泛的應(yīng)用[26]。因其計(jì)算效率高,能夠有效地測(cè)定時(shí)間變化趨勢(shì)的起始位置,并具有檢測(cè)范圍廣、定量化程度高等優(yōu)點(diǎn),比一元線(xiàn)性回歸趨勢(shì)檢驗(yàn)空間變化更為敏感。Theil-Sen斜率估算通過(guò)計(jì)算序列中連續(xù)兩期時(shí)間序列數(shù)據(jù)之間的斜率,將所有像元的數(shù)值對(duì)斜率的中值做時(shí)間序列的總變化趨勢(shì)。公式如下：

(1)

式中,β是所有數(shù)據(jù)對(duì)斜率的中值。其大小表示平均變化率,正負(fù)表示時(shí)間序列的變化趨勢(shì)。Kendall()為取中值函數(shù)。xj和xi分別表示時(shí)間序列中的第j項(xiàng)和第i項(xiàng)的值。但由于該斜率估算方法不能實(shí)現(xiàn)時(shí)間序列的趨勢(shì)顯著性判斷。通常使用Mann-Kendall顯著性檢驗(yàn)來(lái)判斷時(shí)間序列的趨勢(shì)顯著性。

Mann-Kendall方法是由彭曼首次提出并使用,并由Kendall進(jìn)行改進(jìn),對(duì)于序列X=(X1,X2,X3...Xn)先確定所有對(duì)偶值,之后比較與檢驗(yàn)統(tǒng)計(jì)量S的大小關(guān)系。該方法的原始假設(shè)是：H0：時(shí)間序列中的數(shù)據(jù)隨機(jī)排列,即不存在顯著趨勢(shì)；H1：時(shí)間序列數(shù)據(jù)存在上升或下降的單調(diào)趨勢(shì)。公式如下：

(2)

(3)

由于時(shí)間序列維度大小不同,顯著性檢驗(yàn)統(tǒng)計(jì)量的選取有所不同。當(dāng)時(shí)間序列小于10時(shí),直接使用統(tǒng)計(jì)量S進(jìn)行雙邊趨勢(shì)檢驗(yàn)。在給定顯著性水平α下,如果|S|≥Sα/2則拒絕原假設(shè),則存在顯著趨勢(shì)。反之則不存在顯著性趨勢(shì)。當(dāng)時(shí)間序列≥10時(shí),統(tǒng)計(jì)量S近似服從標(biāo)準(zhǔn)正態(tài)分布,可將S標(biāo)準(zhǔn)化后得到檢驗(yàn)統(tǒng)計(jì)量Z,Z值計(jì)算公式如下：

(4)

(5)

式中,n是時(shí)間序列中數(shù)據(jù)個(gè)數(shù)；m是序列中重復(fù)出現(xiàn)的數(shù)據(jù)組個(gè)數(shù)；ti是第i組重復(fù)數(shù)據(jù)組中的重復(fù)數(shù)據(jù)個(gè)數(shù)。同樣,采用雙邊趨勢(shì)檢驗(yàn),在給定顯著性水平α下,在正態(tài)分布表中查得臨界值Z1-α/2。如果|Z|≤Z1-α/2時(shí),則接受原假設(shè)即趨勢(shì)不顯著；反之,則相反。

2.2 生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)物質(zhì)量評(píng)價(jià)

植被所引起的風(fēng)蝕減少量可視為防風(fēng)固沙服務(wù)的物理量,其表示在裸土條件下潛在風(fēng)蝕量與植被覆蓋條件下實(shí)際風(fēng)蝕量之間的差[27]。在充分考慮氣候條件、植被覆蓋、土壤等要素情況下,可根據(jù)RWEQ模型定量評(píng)估潛在風(fēng)蝕量和實(shí)際風(fēng)蝕量。利用兩者的差值進(jìn)行生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)物理量的評(píng)價(jià)[28]。其計(jì)算公式如下：

(6)

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

(7)

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

(8)

(9)

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

(10)

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

(11)

G=SLr-SL

(12)

式中,SLr表示單位面積年潛在風(fēng)蝕量(kg/m2)；Qrmax為潛在風(fēng)力最大輸沙能力(kg/m)；Sr為潛在關(guān)鍵地塊長(zhǎng)度(m)；SL表示單位面積年實(shí)際風(fēng)蝕量(kg/m2)；Qmax為風(fēng)力最大輸沙能力(kg/m);S為關(guān)鍵地塊長(zhǎng)度(m)；G表示單位面積防風(fēng)固沙量(kg/m2);z表示下風(fēng)向距離(m),本文計(jì)算取50m；WF表示風(fēng)蝕因子(kg/m)；EF表示土壤可蝕性因子(無(wú)量綱)；SCF表示為土壤結(jié)皮因子(無(wú)量綱)；K′表示為土壤粗糙度因子(無(wú)量綱)；C表示植被因子(無(wú)量綱)。具體表達(dá)式如下：

風(fēng)蝕因子WF在考慮風(fēng)速作用的前提下,還將積雪覆蓋對(duì)地表的保護(hù)和土壤濕度納入考量,計(jì)算公式為：

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

(13)

wf=u2(u2-u1)2×Nd

(14)

式中,wf表示風(fēng)力系數(shù)(m/s3)；g表示重力加速度(m/s2),本文取9.8 m/s2；ρ表示空氣密度(kg/m3),本文中采用20℃時(shí)的空氣密度值,即1.205kg/m3；SW表示土壤濕度因子(無(wú)量綱)；SD表示積雪覆蓋因子(無(wú)量綱),即無(wú)積雪覆蓋的天數(shù)與研究總天數(shù)之比,積雪深度小于25.4mm時(shí),表示無(wú)積雪；u1為最小起沙風(fēng)速,模型中默認(rèn)為5m/s;u2表示2m高度處的風(fēng)速(m/s)；Nd為每月風(fēng)速超過(guò)閾值風(fēng)速的天數(shù)。

土壤可蝕性因子EF指在一定土壤理化條件下土壤受風(fēng)蝕影響大小。土壤結(jié)皮因子SCF指在一定理化條件下土壤結(jié)皮抵抗風(fēng)蝕能力的大小。其表達(dá)式為：

(15)

(16)

式中,sa為土壤砂粒含量(%);si為土壤粉粒含量(%)；cl為土壤粘粒含量(%)；OM為土壤有機(jī)質(zhì)含量(%)；CaCO3為碳酸鈣含量(%)。假定土壤可蝕性因子和土壤結(jié)皮因子隨時(shí)間保持不變,在應(yīng)用世界土壤數(shù)據(jù)庫(kù)的基礎(chǔ)上,結(jié)合陜北長(zhǎng)城沿線(xiàn)風(fēng)沙邊緣區(qū)土壤質(zhì)地特征,參考RWEQ模型內(nèi)嵌土壤資料表,計(jì)算出土壤可蝕性因子和土壤結(jié)皮因子。

地表粗糙度因子K′是由地形引起的土地表面粗糙程度對(duì)土壤風(fēng)蝕的影響,公式如下：

K′=cosα

(17)

式中,α表示為坡度。

植被因子C表示一定的植被覆蓋對(duì)土壤風(fēng)蝕的抑制程度,公式如下：

C=e-0.0483(FVC)

(18)

(19)

式中,FVC代表植被覆蓋度(%),NDVIveg和NDVIsoil分別表示植被覆蓋的最大和最小值。本文利用植被指數(shù)法,基于像元二分模型計(jì)算得出95%和5%的累積頻率對(duì)應(yīng)下的植被覆蓋最大和最小值[29—30]。

生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)功能可以表明植被的實(shí)際固沙的能力[31]。利用研究區(qū)面積,定量化出防風(fēng)固沙物理量,以此來(lái)評(píng)估區(qū)域生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力的大小。計(jì)算公式為：

QSR=10×G×A

(20)

ASR=QSR/A

(21)

式中,QSR為研究區(qū)防風(fēng)固沙量(t/a)；A為研究區(qū)面積(hm2)；ASR為生態(tài)系統(tǒng)防風(fēng)固沙能力(t/hm2)。

3 結(jié)果與分析

3.1 Theil-Sen斜率估算與M-K檢驗(yàn)

采用Theil-Sen斜率估算和M-K顯著性檢驗(yàn)相結(jié)合的方法,利用NDVI數(shù)據(jù)得到的陜北長(zhǎng)城沿線(xiàn)風(fēng)沙邊緣區(qū)植被變化情況(圖2)。結(jié)果表明：2000—2018年研究區(qū)內(nèi)98.74%的區(qū)域植被變化表現(xiàn)為改善。其中,約有90.34%植被覆蓋區(qū)域的植被狀況呈現(xiàn)出顯著改善。結(jié)合哨兵二號(hào)影像數(shù)據(jù)對(duì)比來(lái)看,榆陽(yáng)區(qū)、神木市等因城市擴(kuò)展,導(dǎo)致植被出現(xiàn)退化。不顯著改善集中表現(xiàn)在定邊縣占比較大,說(shuō)明風(fēng)沙邊緣地區(qū)的植被還需進(jìn)一步保護(hù)。

雖然Theil-Sen斜率估算和M-K檢驗(yàn)方法有計(jì)算效率高、檢測(cè)范圍廣、定量化程度高等優(yōu)點(diǎn),比一元線(xiàn)性回歸趨勢(shì)檢驗(yàn)空間變化更為敏感。但本文檢測(cè)到的陜北長(zhǎng)城沿線(xiàn)風(fēng)沙邊緣地區(qū)植被變化空間分布模式與基于一元線(xiàn)性回歸分析方法檢測(cè)到的植被空間變化模式大體上一致。

3.2 風(fēng)蝕模型因子時(shí)空格局

3.2.1風(fēng)蝕因子結(jié)果

從表1可以看出,2000—2018年風(fēng)蝕因子呈現(xiàn)出減少-增多-減少的趨勢(shì),且風(fēng)蝕因子均值從24.02m/s波動(dòng)式降至12.26m/s,整體上相對(duì)減小。在2010年風(fēng)蝕因子均值達(dá)到最大,為74.13m/s,2005年最低,為6.32m/s。從各點(diǎn)大于5m/s的數(shù)值和天數(shù)來(lái)看,2010年相對(duì)多于其他年份,致使2010年的風(fēng)蝕因子出現(xiàn)最大值。

表1 2000—2018年風(fēng)蝕因子變化

在空間分布來(lái)看(圖3),風(fēng)蝕因子分布差異明顯,且風(fēng)蝕因子高值分布在榆陽(yáng)區(qū)和神木縣周邊。對(duì)比分析,2010年較其他年份差異性更為顯著。在2005年定邊縣風(fēng)蝕因子高于其他五個(gè)縣區(qū),統(tǒng)計(jì)分析這是由于當(dāng)年定邊縣周邊氣象要素點(diǎn)大于起沙風(fēng)速閾值的天數(shù)和個(gè)數(shù)偏多。

圖2 2000—2018年植被變化趨勢(shì)及顯著性檢驗(yàn)Fig.2 Vegetation change trend and significance test from 2000 to 2018

圖3 2000—2018年風(fēng)蝕因子空間分布Fig.3 Spatial distribution of wind erosion factors from 2000 to 2018

3.2.2土壤及地形因子結(jié)果

土壤可蝕性因子表征土壤抵抗侵蝕能力的大小,土壤對(duì)風(fēng)蝕作用敏感程度越大,可蝕性越高。陜北長(zhǎng)城沿線(xiàn)風(fēng)沙邊緣區(qū)EF取值(圖4)在0.26—0.61,平均值為0.48。在空間分布上,榆陽(yáng)地區(qū)土壤對(duì)風(fēng)蝕作用較敏感,EF取值較高,而定邊縣和府谷縣EF取值較低。從土壤理化性質(zhì)出發(fā),紅砂土的EF值最高,這是因?yàn)槠渖沉：窟_(dá)88%以上,容易發(fā)生風(fēng)蝕現(xiàn)象。潛育土EF值最低且為0.266,因其粘粒和有機(jī)質(zhì)含量達(dá)到43.7%相比于其他土壤類(lèi)型不易起沙。

土壤結(jié)皮因子與土壤易受風(fēng)力侵蝕威脅成反比,即土壤表層結(jié)皮越堅(jiān)固,土壤結(jié)皮因子SCF值越大,土壤越不易受到風(fēng)蝕威脅。整個(gè)研究區(qū)SCF的取值(圖4)在0.06—1.00,平均值為0.56。SCF越高抵抗風(fēng)蝕能力越差。從空間分布來(lái)看,SCF高值分布在長(zhǎng)城沿線(xiàn)分割線(xiàn)以北,與內(nèi)蒙古毛烏素沙漠相鄰且抵抗風(fēng)蝕能力較差。從土壤理化性質(zhì)來(lái)看,紅沙土粘粒含量和有機(jī)質(zhì)含量最低,導(dǎo)致SCF最高,為0.837,抗風(fēng)蝕能力最差。始成土的粘粒和有機(jī)質(zhì)含量相對(duì)較高,SCF最低,說(shuō)明其抗風(fēng)蝕能力最強(qiáng),不易發(fā)生風(fēng)蝕。

地表粗糙度因子K′ 反映了地形對(duì)風(fēng)蝕強(qiáng)度大小的影響。長(zhǎng)城沿線(xiàn)K′(圖4)取值范圍為0.004—1.00,均值為0.65?？臻g分布上,長(zhǎng)城沿線(xiàn)分割線(xiàn)以北地表粗糙度占大部分,主要分布在定邊縣北部、靖邊縣、榆陽(yáng)區(qū)、神木縣。根據(jù)海拔高度和地形來(lái)看,定邊縣和靖邊縣海拔北低南高,榆陽(yáng)區(qū)、神木縣西北高東南低,且兩兩縣區(qū)高值區(qū)和低值區(qū)相近,造成地表粗糙度高,更易于侵蝕的發(fā)生。且定邊縣和靖邊縣南部山地海拔較高和起伏較大,不利于侵蝕的發(fā)生。

圖4 土壤因子分布Fig.4 Distribution of soil factors

3.2.3植被因子結(jié)果

植被因子表征植被對(duì)風(fēng)蝕抑制程度的大小,植被生長(zhǎng)越好對(duì)風(fēng)蝕的抑制程度越大進(jìn)而影響生態(tài)系統(tǒng)防風(fēng)固沙能力。研究區(qū)C因子取值范圍為0.007—1.00,植被整體呈現(xiàn)好轉(zhuǎn)趨勢(shì)但均值呈輕微下降趨勢(shì),局部存在波動(dòng)?？臻g分布上(圖5),C因子較高值主要分布在神木市和榆陽(yáng)區(qū)等地,表現(xiàn)出易被風(fēng)力侵蝕,C因子較低值則相反。

圖5 2000—2018植被因子Fig.5 The vegetation factor from 2000 to 2018

3.3 防風(fēng)固沙物質(zhì)量及時(shí)空格局

3.3.1潛在與實(shí)際風(fēng)蝕時(shí)空變化

潛在風(fēng)蝕能夠反映出土壤可能的最大風(fēng)蝕量,即土壤風(fēng)蝕危害的“危險(xiǎn)程度”。從表2可看出,2000—2018年潛在風(fēng)蝕量在4569.18—64164.44萬(wàn)t,2010年達(dá)到潛在風(fēng)蝕量峰值,2005年最少。與2000年相比,2018年減少了13160.93萬(wàn)t?？梢?jiàn)陜北長(zhǎng)城沿線(xiàn)土壤風(fēng)蝕危害的“危險(xiǎn)程度”在呈波動(dòng)式降低。與2000年相比,2018年平均單位面積潛在風(fēng)蝕量下降了3.92kg/m2,總體呈“減少-增加-減少”的波動(dòng)式變化。2000—2018年實(shí)際風(fēng)蝕量在1056.31—24211.89萬(wàn)t,在2010年達(dá)到實(shí)際風(fēng)蝕量峰值(24211.89萬(wàn)t),2005年最少。與2000年相比,2018年增加了894.94萬(wàn)t。與2010年相比,2018年減少了20161.58萬(wàn)t,平均單位面積風(fēng)蝕量減少了6.0kg/m2,減少幅度較大,整體呈“減少-增加-減少”的變化特點(diǎn)。

2000—2018年潛在風(fēng)蝕出現(xiàn)顯著的空間差異性。根據(jù)各個(gè)因子空間分布,在風(fēng)蝕因子分布較大、地表較粗糙、土壤結(jié)皮較大的地區(qū)潛在風(fēng)蝕強(qiáng)度最為劇烈且變化幅度最為明顯,如榆陽(yáng)區(qū)和神木市。而定邊縣、靖邊縣和府谷縣潛在風(fēng)蝕強(qiáng)度較輕。實(shí)際風(fēng)蝕年際空間變化上,也存在明顯的物質(zhì)量變化。依據(jù)水利部《土壤風(fēng)蝕分類(lèi)標(biāo)準(zhǔn)(SL190—2007)》將土壤風(fēng)蝕分為微度(<200t/a)、輕度(200—2500t/a)、中度(2500—5000t/a)、強(qiáng)烈(5000—8000t/a)、極強(qiáng)烈(8000—15000t/a)、劇烈(>15000t/a)[32—33]。根據(jù)分類(lèi)統(tǒng)計(jì),微度土壤風(fēng)蝕面積增加了16.16%,輕度風(fēng)蝕面積減少了15.94%,2010年強(qiáng)烈、極強(qiáng)烈、劇烈風(fēng)蝕面積占比最大,與此對(duì)應(yīng)該年的風(fēng)蝕因子最大,導(dǎo)致風(fēng)蝕現(xiàn)象嚴(yán)重。2010到2018年,風(fēng)蝕因子逐漸減小,土壤風(fēng)蝕強(qiáng)度也逐漸有了明顯的改善。定邊縣、靖邊縣和府谷縣微度風(fēng)蝕強(qiáng)度占比較多；榆陽(yáng)區(qū)和神木市潛在風(fēng)蝕最為劇烈且變化幅度最為明顯,也是單位面積實(shí)際風(fēng)蝕量較嚴(yán)重的區(qū)域(圖6)。從時(shí)間統(tǒng)計(jì)和空間分布比較來(lái)看,每年的潛在風(fēng)蝕量>防風(fēng)固沙量>實(shí)際風(fēng)蝕量,說(shuō)明在最大危害程度內(nèi),防風(fēng)固沙服務(wù)發(fā)揮明顯作用。而從整體情況來(lái)看,在2000—2018年陜北長(zhǎng)城沿線(xiàn)防風(fēng)固沙治理過(guò)程中,土壤風(fēng)蝕現(xiàn)象已得到顯著遏制。相比于潛在風(fēng)蝕,實(shí)際風(fēng)蝕多考慮植被因子,由于其他因定不變,主要是風(fēng)蝕和植被因子對(duì)風(fēng)蝕的影響最密切。

表2 2000—2018年防風(fēng)固沙物質(zhì)量統(tǒng)計(jì)/(kg/m2)

圖6 2000—2018年風(fēng)蝕強(qiáng)度空間分布Fig.6 Spatial distribution of wind erosion intensity from 2000 to 2018

3.3.2防風(fēng)固沙服務(wù)物質(zhì)量時(shí)空變化

潛在與實(shí)際風(fēng)蝕差值表征土壤風(fēng)蝕模數(shù)(單位面積防風(fēng)固沙物理量)。研究表明：防風(fēng)固沙物質(zhì)量從2000年的21944.56萬(wàn)t減少到2018年的7987.44萬(wàn)t。2010年變化波動(dòng)增長(zhǎng)較大且達(dá)到最大固沙量。依據(jù)張彪等[31]將防風(fēng)固沙服務(wù)能力ASR分為：低值區(qū)(ASR≤5t/hm2)、較低區(qū)(5≤ASR≤20t/hm2)、一般區(qū)(20≤ASR≤60t/hm2)、較高區(qū)(60≤ASR≤120t/hm2)和高值區(qū)(ASR≥120t/hm2)。從防風(fēng)固沙服務(wù)能力分區(qū)統(tǒng)計(jì)來(lái)看,如表3所示,2000—2018年生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力高值區(qū)和較高值區(qū)出現(xiàn)減小-增多-減小的變化特點(diǎn),且整體上是減小的趨勢(shì),且對(duì)比觀(guān)察不同年份的防風(fēng)固沙強(qiáng)度的分布規(guī)律即存在相似又存在差異。長(zhǎng)城沿線(xiàn)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力表現(xiàn)出明顯的空間異質(zhì)性,防風(fēng)固沙服務(wù)能力高值區(qū)主要集中在神木市和榆陽(yáng)區(qū),如圖7所示。根據(jù)模型分析,如果風(fēng)蝕因子偏大且植被因子降低,即風(fēng)場(chǎng)強(qiáng)度增加導(dǎo)致防風(fēng)固沙量顯著減少。風(fēng)速的不穩(wěn)定性會(huì)導(dǎo)致風(fēng)蝕因子結(jié)果出現(xiàn)較大偏差,進(jìn)而對(duì)防風(fēng)固沙服務(wù)的空間分布變化產(chǎn)生較大的影響。在抑制土壤風(fēng)蝕現(xiàn)象發(fā)生的基礎(chǔ)上,固沙能力也在減小,后期需要加以重點(diǎn)保護(hù)植被。

表3 2000—2018年防風(fēng)固沙能力空間分布/km2

圖7 2000—2018年防風(fēng)固沙能力空間分布Fig.7 Spatial distribution of wind prevention and sand fixation capacity from 2000 to 2018

3.4 植被恢復(fù)對(duì)土壤風(fēng)蝕及防風(fēng)固沙服務(wù)影響

為了進(jìn)一步探究風(fēng)蝕因子變化,植被恢復(fù)對(duì)風(fēng)蝕現(xiàn)象及生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)的影響,利用情景模擬法,以最小(2005年)、最大風(fēng)蝕(2010年)因子為基準(zhǔn),分設(shè)兩種情景(即情景1和2)。由于計(jì)算潛在風(fēng)蝕時(shí)沒(méi)有考慮到植被覆蓋,其他因子固定不變,即：2005年的潛在風(fēng)蝕為情景1條件下的潛在風(fēng)蝕,2010年的潛在風(fēng)蝕為情景2條件下的潛在風(fēng)蝕。從物質(zhì)定量化分析角度出發(fā)(表4),橫向?qū)Ρ惹榫?得出：最小風(fēng)蝕因子情況下,植被長(zhǎng)時(shí)間序列恢復(fù),植被因子降低,抵抗風(fēng)蝕能力增強(qiáng),實(shí)際風(fēng)蝕量從1066.56萬(wàn)t減少到825.38萬(wàn)t,相比減少了241.18萬(wàn)t,防風(fēng)固沙量在逐年提升的同時(shí)防風(fēng)固沙服務(wù)能力也顯著提升。與情景1相比,情景2中風(fēng)蝕因子的改變使得土壤風(fēng)蝕量從2000年的11764.29萬(wàn)t波浪式增加到2018年的22649.44萬(wàn)t,說(shuō)明植被不再是影響風(fēng)蝕的主要因素,計(jì)算風(fēng)蝕因子的土壤濕度、風(fēng)力系數(shù)等因子的改變使得地表更不易抵抗風(fēng)蝕。但縱向?qū)Ρ惹榫?和2：隨著風(fēng)蝕因子的顯著變化,土壤風(fēng)蝕危害的“危險(xiǎn)程度”達(dá)到最大(潛在風(fēng)蝕),以2000年為例,實(shí)際風(fēng)蝕量從1066.56萬(wàn)t顯著增長(zhǎng)到11764.29萬(wàn)t,防風(fēng)固沙量從3449.6萬(wàn)t變到52175.12}。由此可見(jiàn)：其他因素不變,氣象條件惡劣的情況下,極易發(fā)生土壤風(fēng)蝕。但兩種情境下的每年防風(fēng)固沙量>實(shí)際風(fēng)蝕量,利用植被固沙的效果也越來(lái)越明顯,生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)功能也在增加。

表4 2000—2018年分情景下防風(fēng)固沙物質(zhì)量統(tǒng)計(jì)/(kg/m2)

空間分布上(圖8),在風(fēng)蝕因子條件下(情景1),植被恢復(fù),2000—2018年風(fēng)蝕現(xiàn)象逐漸退化。在最大風(fēng)蝕因子下(情景2),造成風(fēng)蝕現(xiàn)象嚴(yán)重,且主要集中在榆陽(yáng)區(qū)和神木縣,導(dǎo)致定邊縣和靖邊縣受到殃及。這四個(gè)縣區(qū)的防風(fēng)固沙量也出現(xiàn)明顯差異?？梢?jiàn)風(fēng)蝕因子是主導(dǎo)土壤風(fēng)蝕的主導(dǎo)因素。整體來(lái)看：采取“鎖綠”來(lái)防沙治沙是最為有效的,且隨著風(fēng)蝕因子的改變水土流失治理成效越為明顯,生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力也越強(qiáng)。植被恢復(fù)對(duì)長(zhǎng)城沿線(xiàn)的土壤風(fēng)蝕起到有效的抑制,對(duì)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)起到了一定的增強(qiáng)作用。但風(fēng)蝕因子中風(fēng)速屬于不可防控的氣象因素,只能依靠恢復(fù)植被來(lái)增強(qiáng)生態(tài)系統(tǒng)防風(fēng)固沙能力。

圖8 分情境下風(fēng)蝕空間分布Fig.8 The amount of wind erosion under different scenarios

4 結(jié)論

本文基于RWEQ模型分析了陜北長(zhǎng)城沿線(xiàn)防風(fēng)固沙服務(wù)物理量的時(shí)空模式,實(shí)現(xiàn)風(fēng)沙邊緣地區(qū)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)的量化,并評(píng)估風(fēng)蝕和植被因子變化下對(duì)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)物理量和服務(wù)能力的影響。得到以下結(jié)論：(1)通過(guò)分析2000—2018年陜北長(zhǎng)城沿線(xiàn)六個(gè)縣區(qū)年際NDVI的時(shí)空變化,風(fēng)沙區(qū)內(nèi)植被呈明顯增加趨勢(shì)；(2)2000—2018年整個(gè)風(fēng)沙區(qū)實(shí)現(xiàn)了從“沙進(jìn)人退”到“綠進(jìn)沙退”的轉(zhuǎn)變,土壤風(fēng)蝕危害的“危險(xiǎn)程度”也在逐步降低,實(shí)際風(fēng)蝕和生態(tài)系統(tǒng)防風(fēng)固沙總量2005年和2010年分別為4569.18萬(wàn)t和64164.44萬(wàn)t,且潛在風(fēng)蝕量>防風(fēng)固沙量>實(shí)際風(fēng)蝕量,生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)發(fā)揮顯著作用。(3)利用情景變化定量分析,探索植被恢復(fù)可以一定程度上抑制風(fēng)蝕現(xiàn)象。在最小風(fēng)蝕因子情況下,不受外力條件變化,植被防風(fēng)固沙服務(wù)的作用發(fā)揮不明顯,且一定程度抑制風(fēng)蝕現(xiàn)象的發(fā)生。風(fēng)蝕因子改變加重研究區(qū)的土壤風(fēng)蝕,植被所發(fā)揮的生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)越為明顯。而兩因子相結(jié)合(即實(shí)際情況下的風(fēng)蝕量),風(fēng)蝕因子占主導(dǎo)地位,對(duì)土壤風(fēng)蝕和防風(fēng)固沙量影響作用較大。不論是真實(shí)情況還是分情景下的風(fēng)蝕,生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)空間差異顯著,榆陽(yáng)區(qū)和神木市的土壤風(fēng)蝕現(xiàn)象最為嚴(yán)重,也是生態(tài)系統(tǒng)防風(fēng)固沙服物質(zhì)量和能力高值集中區(qū)。總之,恢復(fù)植被是實(shí)現(xiàn)區(qū)域性荒漠化逆轉(zhuǎn)和改善生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)與能力的必要條件,恢復(fù)植被不僅可以減少土壤風(fēng)蝕,也可以增強(qiáng)生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)能力。

5 討論

(1)當(dāng)前土壤風(fēng)力侵蝕的定量化驗(yàn)證還存在一定的困難。實(shí)測(cè)的土壤風(fēng)蝕量是開(kāi)展研究結(jié)果驗(yàn)證的基礎(chǔ),由于沒(méi)有實(shí)測(cè)的土壤風(fēng)蝕量數(shù)據(jù),無(wú)法對(duì)利用模型模擬的實(shí)際風(fēng)蝕量結(jié)果進(jìn)行驗(yàn)證,后續(xù)可進(jìn)一步采用沙塵暴發(fā)生頻率來(lái)驗(yàn)證本文的可靠性。鑒于此,文中利用情景模擬法可以間接進(jìn)行驗(yàn)證,對(duì)比分析植被恢復(fù)對(duì)風(fēng)蝕現(xiàn)象有一定的遏制作用,與劉碩等[16]的研究結(jié)果一致。在改變風(fēng)力因子時(shí),風(fēng)蝕現(xiàn)象依舊明顯,植被發(fā)揮的生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)也顯著提高。

(2)提高模型的優(yōu)化性。在實(shí)現(xiàn)陜北長(zhǎng)城沿線(xiàn)2000—2018年生態(tài)系統(tǒng)防風(fēng)固沙量以及服務(wù)能力,由于實(shí)測(cè)氣象站點(diǎn)少,而采用中國(guó)區(qū)域地面要素?cái)?shù)據(jù)集,但該數(shù)據(jù)集時(shí)間范圍限制,未能及時(shí)做到最新年份；考慮到該數(shù)據(jù)集時(shí)間分辨率較大,且有349個(gè)規(guī)則分布的點(diǎn),在選用空間分辨率時(shí)采用1km空間分辨率,但對(duì)于土壤風(fēng)蝕研究而言,1km分辨率足以反映其空間差異。在實(shí)現(xiàn)雪蓋因子計(jì)算的過(guò)程中,計(jì)算的是年積雪覆蓋天數(shù)的概率,得到年尺度風(fēng)蝕量。在后期研究細(xì)化時(shí)間尺度,從月尺度考慮生態(tài)系統(tǒng)防風(fēng)固沙物質(zhì)量和服務(wù)能力。

(3)本文的研究結(jié)論僅從生態(tài)系統(tǒng)防風(fēng)固沙量和服務(wù)能力角度出發(fā),在情景1條件下植被的固沙作用效果不明顯,但在最大風(fēng)蝕因子條件下,植被對(duì)抑制風(fēng)沙吹動(dòng),防風(fēng)固沙效果越來(lái)越明顯。三種情景每年的潛在風(fēng)蝕量>防風(fēng)固沙量>實(shí)際風(fēng)蝕量。風(fēng)蝕因子增大,潛在風(fēng)蝕增加、實(shí)際風(fēng)蝕增加、防風(fēng)固沙量增加,生態(tài)系統(tǒng)防風(fēng)固沙服務(wù)功能也增強(qiáng)。這與治沙研究者提出防沙治沙、植樹(shù)造林工程以來(lái)的預(yù)想一致。從“沙進(jìn)人退”到“綠進(jìn)沙退”的轉(zhuǎn)變是數(shù)十年治沙的功勞,“榆林綠”鎖住了毛烏素沙地。作為自然地貌,毛烏素沙地不是真的消失,只是流動(dòng)沙丘全部得到固定。但這些治理成果還屬人工干預(yù),尚需加快沙區(qū)生態(tài)系統(tǒng)的正向演替,實(shí)現(xiàn)自我循環(huán)發(fā)展,從根本上改變沙區(qū)生態(tài)面貌。

本文研究特色在于：采用中國(guó)區(qū)域地面要素?cái)?shù)據(jù)集,更好的實(shí)現(xiàn)長(zhǎng)城沿線(xiàn)風(fēng)蝕因子的空間分布。但還存在不足之處：在考慮土壤理化性質(zhì)隨時(shí)間不發(fā)生變化的情況下,但實(shí)際情況中,由于陜北長(zhǎng)城沿線(xiàn)六個(gè)縣區(qū)耕地面積占比較大,長(zhǎng)期使用化肥致使土壤理化性質(zhì)發(fā)生輕微改變,導(dǎo)致文章實(shí)驗(yàn)結(jié)果的可靠性存在一定的誤差,土壤理化性質(zhì)的改變對(duì)土壤風(fēng)蝕具有一定的影響,對(duì)風(fēng)沙區(qū)土壤理化性質(zhì)的測(cè)定研究也是今后的研究的一個(gè)重要科學(xué)問(wèn)題。在可控的人為范圍內(nèi),做好治沙工程,堅(jiān)持不懈地防沙治沙、植樹(shù)造林使得生態(tài)環(huán)境持續(xù)改善,為榆林經(jīng)濟(jì)社會(huì)的快速發(fā)展提供了強(qiáng)有力的基礎(chǔ)保障。