陜北長城沿線植被恢復與生態系統防風固沙服務模擬分析

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

西安科技大學測繪科學與技術學院,西安 710054

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

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

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

1 研究區及數據源

1.1 研究區概況

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

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

1.2 數據來源

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

2 研究方法

2.1 Theil-Sen斜率估算與M-K顯著性檢驗結合

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

(1)

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

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

(2)

(3)

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

(4)

(5)

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

2.2 生態系統防風固沙服務物質量評價

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

(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表示單位面積年潛在風蝕量(kg/m2)；Qrmax為潛在風力最大輸沙能力(kg/m)；Sr為潛在關鍵地塊長度(m)；SL表示單位面積年實際風蝕量(kg/m2)；Qmax為風力最大輸沙能力(kg/m);S為關鍵地塊長度(m)；G表示單位面積防風固沙量(kg/m2);z表示下風向距離(m),本文計算取50m；WF表示風蝕因子(kg/m)；EF表示土壤可蝕性因子(無量綱)；SCF表示為土壤結皮因子(無量綱)；K′表示為土壤粗糙度因子(無量綱)；C表示植被因子(無量綱)。具體表達式如下：

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

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

(13)

wf=u2(u2-u1)2×Nd

(14)

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

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

(15)

(16)

式中,sa為土壤砂粒含量(%);si為土壤粉粒含量(%)；cl為土壤粘粒含量(%)；OM為土壤有機質含量(%)；CaCO3為碳酸鈣含量(%)。假定土壤可蝕性因子和土壤結皮因子隨時間保持不變,在應用世界土壤數據庫的基礎上,結合陜北長城沿線風沙邊緣區土壤質地特征,參考RWEQ模型內嵌土壤資料表,計算出土壤可蝕性因子和土壤結皮因子。

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

K′=cosα

(17)

式中,α表示為坡度。

植被因子C表示一定的植被覆蓋對土壤風蝕的抑制程度,公式如下：

C=e-0.0483(FVC)

(18)

(19)

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

生態系統防風固沙服務功能可以表明植被的實際固沙的能力[31]。利用研究區面積,定量化出防風固沙物理量,以此來評估區域生態系統防風固沙服務能力的大小。計算公式為：

QSR=10×G×A

(20)

ASR=QSR/A

(21)

式中,QSR為研究區防風固沙量(t/a)；A為研究區面積(hm2)；ASR為生態系統防風固沙能力(t/hm2)。

3 結果與分析

3.1 Theil-Sen斜率估算與M-K檢驗

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

雖然Theil-Sen斜率估算和M-K檢驗方法有計算效率高、檢測范圍廣、定量化程度高等優點,比一元線性回歸趨勢檢驗空間變化更為敏感。但本文檢測到的陜北長城沿線風沙邊緣地區植被變化空間分布模式與基于一元線性回歸分析方法檢測到的植被空間變化模式大體上一致。

3.2 風蝕模型因子時空格局

3.2.1風蝕因子結果

從表1可以看出,2000—2018年風蝕因子呈現出減少-增多-減少的趨勢,且風蝕因子均值從24.02m/s波動式降至12.26m/s,整體上相對減小。在2010年風蝕因子均值達到最大,為74.13m/s,2005年最低,為6.32m/s。從各點大于5m/s的數值和天數來看,2010年相對多于其他年份,致使2010年的風蝕因子出現最大值。

表1 2000—2018年風蝕因子變化

在空間分布來看(圖3),風蝕因子分布差異明顯,且風蝕因子高值分布在榆陽區和神木縣周邊。對比分析,2010年較其他年份差異性更為顯著。在2005年定邊縣風蝕因子高于其他五個縣區,統計分析這是由于當年定邊縣周邊氣象要素點大于起沙風速閾值的天數和個數偏多。

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

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

3.2.2土壤及地形因子結果

土壤可蝕性因子表征土壤抵抗侵蝕能力的大小,土壤對風蝕作用敏感程度越大,可蝕性越高。陜北長城沿線風沙邊緣區EF取值(圖4)在0.26—0.61,平均值為0.48。在空間分布上,榆陽地區土壤對風蝕作用較敏感,EF取值較高,而定邊縣和府谷縣EF取值較低。從土壤理化性質出發,紅砂土的EF值最高,這是因為其沙粒含量達88%以上,容易發生風蝕現象。潛育土EF值最低且為0.266,因其粘粒和有機質含量達到43.7%相比于其他土壤類型不易起沙。

土壤結皮因子與土壤易受風力侵蝕威脅成反比,即土壤表層結皮越堅固,土壤結皮因子SCF值越大,土壤越不易受到風蝕威脅。整個研究區SCF的取值(圖4)在0.06—1.00,平均值為0.56。SCF越高抵抗風蝕能力越差。從空間分布來看,SCF高值分布在長城沿線分割線以北,與內蒙古毛烏素沙漠相鄰且抵抗風蝕能力較差。從土壤理化性質來看,紅沙土粘粒含量和有機質含量最低,導致SCF最高,為0.837,抗風蝕能力最差。始成土的粘粒和有機質含量相對較高,SCF最低,說明其抗風蝕能力最強,不易發生風蝕。

地表粗糙度因子K′ 反映了地形對風蝕強度大小的影響。長城沿線K′(圖4)取值范圍為0.004—1.00,均值為0.65。空間分布上,長城沿線分割線以北地表粗糙度占大部分,主要分布在定邊縣北部、靖邊縣、榆陽區、神木縣。根據海拔高度和地形來看,定邊縣和靖邊縣海拔北低南高,榆陽區、神木縣西北高東南低,且兩兩縣區高值區和低值區相近,造成地表粗糙度高,更易于侵蝕的發生。且定邊縣和靖邊縣南部山地海拔較高和起伏較大,不利于侵蝕的發生。

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

3.2.3植被因子結果

植被因子表征植被對風蝕抑制程度的大小,植被生長越好對風蝕的抑制程度越大進而影響生態系統防風固沙能力。研究區C因子取值范圍為0.007—1.00,植被整體呈現好轉趨勢但均值呈輕微下降趨勢,局部存在波動。空間分布上(圖5),C因子較高值主要分布在神木市和榆陽區等地,表現出易被風力侵蝕,C因子較低值則相反。

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

3.3 防風固沙物質量及時空格局

3.3.1潛在與實際風蝕時空變化

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

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

表2 2000—2018年防風固沙物質量統計/(kg/m2)

圖6 2000—2018年風蝕強度空間分布Fig.6 Spatial distribution of wind erosion intensity from 2000 to 2018

3.3.2防風固沙服務物質量時空變化

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

表3 2000—2018年防風固沙能力空間分布/km2

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

3.4 植被恢復對土壤風蝕及防風固沙服務影響

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

表4 2000—2018年分情景下防風固沙物質量統計/(kg/m2)

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

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

4 結論

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

5 討論

(1)當前土壤風力侵蝕的定量化驗證還存在一定的困難。實測的土壤風蝕量是開展研究結果驗證的基礎,由于沒有實測的土壤風蝕量數據,無法對利用模型模擬的實際風蝕量結果進行驗證,后續可進一步采用沙塵暴發生頻率來驗證本文的可靠性。鑒于此,文中利用情景模擬法可以間接進行驗證,對比分析植被恢復對風蝕現象有一定的遏制作用,與劉碩等[16]的研究結果一致。在改變風力因子時,風蝕現象依舊明顯,植被發揮的生態系統防風固沙服務也顯著提高。

(2)提高模型的優化性。在實現陜北長城沿線2000—2018年生態系統防風固沙量以及服務能力,由于實測氣象站點少,而采用中國區域地面要素數據集,但該數據集時間范圍限制,未能及時做到最新年份；考慮到該數據集時間分辨率較大,且有349個規則分布的點,在選用空間分辨率時采用1km空間分辨率,但對于土壤風蝕研究而言,1km分辨率足以反映其空間差異。在實現雪蓋因子計算的過程中,計算的是年積雪覆蓋天數的概率,得到年尺度風蝕量。在后期研究細化時間尺度,從月尺度考慮生態系統防風固沙物質量和服務能力。

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

本文研究特色在于：采用中國區域地面要素數據集,更好的實現長城沿線風蝕因子的空間分布。但還存在不足之處：在考慮土壤理化性質隨時間不發生變化的情況下,但實際情況中,由于陜北長城沿線六個縣區耕地面積占比較大,長期使用化肥致使土壤理化性質發生輕微改變,導致文章實驗結果的可靠性存在一定的誤差,土壤理化性質的改變對土壤風蝕具有一定的影響,對風沙區土壤理化性質的測定研究也是今后的研究的一個重要科學問題。在可控的人為范圍內,做好治沙工程,堅持不懈地防沙治沙、植樹造林使得生態環境持續改善,為榆林經濟社會的快速發展提供了強有力的基礎保障。