張 彪,王 爽

(1.中國科學院地理科學與資源研究所,北京 100101;2.中國科學院大學資源與環(huán)境學院,北京 100049)

土地沙化是全球面臨的重大生態(tài)環(huán)境問題,其與地形、土壤、植被、氣候和人類活動等因素有關[1-2]。近年來我國土地退化態(tài)勢較嚴重。根據2015年發(fā)布的《中國荒漠化和沙化狀況公報》,截至2014年,我國仍有荒漠化土地261.16萬km2和沙化土地172.12萬km2。為遏制土地沙化趨勢,我國部署實施了退耕還林還草、三北防護林體系、天然林保護和京津風沙源治理等重大生態(tài)建設工程。由于生態(tài)建設工程工期長，地域廣,加上生態(tài)系統(tǒng)恢復漫長復雜,及時監(jiān)測、評估工程區(qū)生態(tài)狀況及功能變化尤為重要[3]。京津風沙源治理工程是為遏制北京及周邊地區(qū)土地沙化趨勢、改善京津地區(qū)大氣環(huán)境質量而實施的一項重大生態(tài)工程,建設范圍西起內蒙古達爾罕茂明安聯(lián)合旗,東至河北省平泉縣,南至山西省代縣,北至內蒙古東烏珠穆沁旗,涉及國土面積45.8萬km2。自2000年京津風沙源治理工程實施以來,工程區(qū)植被恢復與土壤風蝕趨勢受到重點關注。李慶旭等[3]基于MOD13Q1數據測算2000—2015年京津風沙源區(qū)植被覆蓋度變化范圍在35%～45%之間,且以年均0.4%速率波動增加。石莎等[4]利用傳統(tǒng)生態(tài)學方法對京津風沙源治理工程區(qū)2001—2005年植被恢復效果進行野外調查,發(fā)現植被覆蓋度整體上升。嚴恩萍等[5]采用MODIS02B數據源分析了重點治理區(qū)(涉及山西和內蒙古共48個縣域)植被覆蓋變化特征,表明2000—2012年治理區(qū)植被覆蓋總體呈上升趨勢。衛(wèi)潔等[6]利用MODIS數據監(jiān)測發(fā)現,2000—2014年京津風沙源區(qū)84.51%的地區(qū)植被覆蓋呈增加趨勢。生態(tài)工程區(qū)植被狀況變化必然引起區(qū)域防風固沙功能的相應變化。高尚玉等[7]調查評估發(fā)現,2005年京津風沙源區(qū)土壤風蝕量比2001年減少1.96億t。吳丹等[8]綜合評估表明,2000—2010年京津風沙源治理工程區(qū)土壤風蝕以微度和輕度為主,且風蝕量總體呈下降態(tài)勢。遲文峰等[9]以遙感手段監(jiān)測發(fā)現,1990—2015年間工程區(qū)土壤風蝕模數整體下降,且在生態(tài)工程實施后改善趨勢明顯。1990—2015年內蒙古49.23%的區(qū)域土壤風蝕強度下降,但仍有7.11%的區(qū)域風蝕強度增加[10]。然而,相比于工程區(qū)植被覆蓋狀況的監(jiān)測,區(qū)域防風固沙功能的動態(tài)變化研究仍相對薄弱。

錫林郭勒盟位于京津風沙源治理工程區(qū)西北部,是工程區(qū)內面積最大的地級行政區(qū)。錫林郭勒盟中部是東西長約400 km的渾善達克沙地,其西緣為固定、半固定沙丘廣布的嘎亥額勒蘇沙地,其生態(tài)治理效果成為評估京津風沙源治理工程成效的關鍵。鞏國麗等[11]運用修正風蝕方程(RWEQ)評估發(fā)現,20世紀90年代以來錫林郭勒盟土壤風蝕強度總體減弱。丁文廣等[12]利用Landsat衛(wèi)星影像分析發(fā)現,1975—2014年錫林郭勒盟沙漠化面積先增后減,沙漠化重心先東移后西返。王艷琦等[13]采用Landsat-8 OLI遙感影像分析認為,2000—2015年錫林郭勒盟沙漠化程度總體呈現逆轉趨勢。此外,渾善達克沙漠化防治生態(tài)功能區(qū)的防風固沙功能受到關注[14]。但是,以往研究側重錫林郭勒盟土壤風蝕狀況的變化趨勢及驅動力分析,而對區(qū)域防風固沙功能的動態(tài)變化及其空間差異揭示不足。

隨著京津風沙源治理二期工程(2013—2022年)的實施,分區(qū)施策與精準修復已成為生態(tài)治理工程的重要內容[15]。為此,以錫林郭勒盟為例,采用REWQ模型與GIS空間分析技術,評估分析2000—2015年錫林郭勒盟防風固沙功能的時空變化及其區(qū)域差異,并探討植被覆蓋度和降水量對區(qū)域防風固沙功能的影響,可為風沙源治理工程的分區(qū)施策治理與成效監(jiān)測評估提供參考依據。

1 研究區(qū)概況

錫林郭勒盟位于京津風沙源治理工程區(qū)西北部、內蒙古中部(圖1),是典型農牧交錯地區(qū),也是北京和華北地區(qū)的重要生態(tài)屏障[16]。錫林郭勒盟土地總面積為202 580 km2,現轄9旗(東烏珠穆沁旗、西烏珠穆沁旗、阿巴嘎旗、蘇尼特左旗、蘇尼特右旗、正藍旗、正鑲白旗、鑲黃旗、太仆寺旗)、2市(錫林浩特市、二連浩特市)和1縣(多倫縣)。

圖1 錫林郭勒盟位置與組成

錫林郭勒盟地形平坦開闊,地勢自西南向東北方向傾斜。海拔高度處于800～1 800 m之間[16],海拔最高處位于太仆寺旗東部,最低處位于東烏珠穆沁旗和西烏珠穆沁旗。錫林郭勒盟氣候屬于中溫帶干旱、半干早大陸性季風氣候區(qū),年均氣溫為1～2 ℃,年平均降水量為200～350 mm,且集中在6—8月[17]。錫林郭勒盟主導風向為西北風,全年大風日數達50～80 d,其中3—5月大風日數占全年的40%～50%[18]。錫林郭勒盟地帶性植被為草原,其面積約占草地總面積的71.88 %,類型包括典型草原、草甸草原和荒漠草原3大亞型[19]。據錫林郭勒盟2018年統(tǒng)計年鑒,2017年末全盟常住人口為105.16萬人,自然增長率為4.1%。其中,城鎮(zhèn)常住人口為68.49萬人,占總人口比例為65.13%;農村牧區(qū)常住人口為36.67萬人,占總人口比例為34.87%[20]。

2 材料與方法

2.1 RWEQ模型的構建

土壤風蝕是風力作用導致表土物質脫離原空間位置的過程[21]。風蝕模型是評估土壤風蝕狀況的主要技術手段,并已發(fā)展演化出通用風蝕方程(WEQ)、德克薩斯侵蝕分析模型(TEAM)和風蝕預報系統(tǒng)(WEPS)等多種模型[22]。其中,修正風蝕方程(RWEQ)綜合反映了氣候、植被、土壤性質和地表粗糙度等自然因素的影響,已被廣泛應用于我國土壤風蝕狀況評估[9-11,14]。防風固沙功能是在地表無植被狀況下土壤風蝕量與植被覆蓋條件下土壤風蝕量的差值[10,23],采用RWEQ模型[22]定量評估錫林郭勒盟潛在風蝕量(LS,p)與實際風蝕量(LS,r),以兩者之差來表示區(qū)域防風固沙量(SRQ,QSR)與平均固沙能力(SRA,ASR)。具體計算公式為

Qp，max=109.8×(FW×FE×FSC×K′)，

(1)

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

(2)

(3)

Qr，max=109.8×(FW×FE×FSC×K′×C)，

(4)

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

(5)

(6)

QSR=10×(LS，p-LS，r)×A，

(7)

ASR=QSR/A。

(8)

式(1)～(8)中,LS,p為潛在風蝕量,kg·m-2;Qp,max為潛在風力最大輸沙能力,kg·m-1;rs為潛在關鍵地塊長度,m;LS,r為實際風蝕量,kg·m-2;Qr,max為實際風力最大輸沙能力,kg·m-1;ps為實際關鍵地塊長度,m;z為下風向距離,取50 m;FW為氣象因子,kg·m-1;FE和FSC分別為土壤可蝕性因子和土壤結皮因子;K′和C分別為土壤糙度因子與植被因子;QSR為研究區(qū)防風固沙量,t·a-1;A為研究區(qū)面積,hm2;ASR為防風固沙能力,t·hm-2·a-1。

2.2 RWEQ模型參數處理

(1)氣象因子(WF)

自然條件下土壤風蝕受風速、溫度、降雨、太陽輻射和降雪等氣象因素影響,氣象因子為各類氣象因素對風蝕的綜合影響,計算公式為

FW=fW×(ρ/g)×SW×SD,

(9)

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

(10)

式(9)～(10)中,fW為風力因子,m3·s-3;ρ為空氣密度,當氣溫為20 ℃時為1.205 kg·m-3;g為重力加速度,取9.8 m·s-2;SW和SD分別為土壤濕度因子和雪蓋因子,計算方法參照文獻[9-11];u1為起沙風速,參照以往研究[11,16,22]取5 m·s-1;u2為氣象站月均風速,m·s-1;Nd為各月風速大于5 m·s-1的時間,d。

(2)土壤可蝕性因子(EF)和結皮因子(SCF)

土壤可蝕性受土壤顆粒粒徑以及有機質、黏土和碳酸鈣等物質含量的影響,表層堅硬結皮也能有效防止風蝕的發(fā)生[23]。因此,可從土壤理化條件判別土壤可蝕性因子,而土壤結皮因子可反映一定土壤理化條件下土壤結皮抵抗風蝕的能力[22],其計算公式分別為

FE=[29.09+0.31wsa+0.17wsi+0.33(wsa/

wcl)-2.59wOM-0.95wCaCO3]/100,

(11)

FSC=1/(1+0.006 6wcl2+0.021wOM2)。

(12)

式(11)～(12)中,wsa為土壤粗砂質量分數,%;wsi為土壤粉砂質量分數,%;wcl為土壤黏粒質量分數,%;wOM為土壤有機質質量分數,%;wCaCO3為碳酸鈣質量分數,%。

(3)植被覆蓋因子(C)

地表植被不僅增加地表糙度而增大起沙風速,且對土壤顆粒移動有一定阻礙作用,因此對土壤風蝕具有重要抑制作用[23]。植被覆蓋因子代表植被條件對風蝕的抑制程度,計算公式為

C=e-0.048 3CS,

(13)

CS=(INDV-INDV，min)/(INDV，max-INDV，min)。

(14)

式(13)～(14)中,CS為植被覆蓋度,%;INDV,max與INDV,min分別為歸一化植被指數(NDVI，INDV)最大值和最小值。

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

坡度、坡向等地形因子可影響植被生長狀況,對風蝕過程也存在顯著影響。地表糙度指農田因耕作產生塊狀土以及土壟存在而對土壤風蝕產生的影響[23],包括隨機糙度(Crr)和土壟糙度(Kr)。由于區(qū)域尺度評估中耕作產生的隨機糙度難以獲取,采用smith-carson方程計算土壟造成的地形起伏度來替代[24],計算公式為

Kr=0.2(ΔH)2/L，

(15)

K′=e(1.86Kr-2.41Kr0.934-0.127Crr)。

(16)

式(15)～(16)中,Kr為土壟糙度,cm;Crr為隨機糙度,取0;K′為地表糙度因子,cm;L為地勢起伏參數;ΔH為距離L范圍內海拔高度差,cm。

2.3 數據來源

氣象數據來源于中國氣象科學數據共享服務網(http:∥cdc.cma.gov.cn/),在錫林郭勒盟轄區(qū)分布有10個氣象站,選取2000—2015年共16 a氣象觀測數據,采用月均風速、降水、氣溫和日照時數等數據插值得到風力因子和土壤濕度因子。雪蓋因子利用中國西部環(huán)境與生態(tài)科學數據中心(http:∥westdc.westgis.ac.cn)的中國雪深長時間序列數據集計算。土壤數據來源于中國西部環(huán)境與生態(tài)科學數據中心提供的1∶100萬土壤圖及所附的土壤屬性表和空間數據,并分別采用錫林郭勒盟不同土壤類型的相應物質含量估算可蝕性因子與結皮因子。NDVI數據來自美國地球資源觀測系統(tǒng)數據中心的MOD13Q1產品,該數據已經過幾何精糾正、輻射校正和大氣校正等預處理,時空分辨率分別為16 d和250 m。對該數據集去除噪聲干擾后,利用MRT投影轉換工具進行投影和格式轉換批處理,并采用最大值合成法獲得2000—2015年NDVI數據。為保證不同精度數據的一致性,以上參數因子均重采樣為100 m×100 m的柵格單元參與模型計算。

3 結果與分析

3.1 防風固沙功能年際變化特征

評估結果表明,2000—2015年錫林郭勒盟年均防風固沙量為12.39億～16.3億t,并分別于2001和2006年達到最低值和最高值,多年防風固沙量平均為14.56億t。評估期內,錫林郭勒盟年均防風固沙量波動增加趨勢顯著(圖2),平均年增速為0.071億t,年均變化率為7.1%。

錫林郭勒盟防風固沙能力2000—2015年在67.13～81.86 t·hm-2之間變動,多年平均值為74.78 t·hm-2。且隨著評估年份的變化,錫林郭勒盟防風固沙能力以平均每年6.3%的變化率波動增加,直觀表現為區(qū)域土壤風蝕強度的整體下降(圖3)。但區(qū)域防風固沙能力增加幅度稍低于防風固沙量,主要原因是防風固沙量受到防風固沙能力與區(qū)域面積的雙重影響。

圖2 2000—2015年錫林郭勒盟防風固沙功能變化

3.2 防風固沙功能空間差異特征

為揭示錫林郭勒盟防風固沙功能空間差異,依據全盟固沙能力值域及風力侵蝕強度分級[25]設置相應防風固沙功能分級標準(表1)。

結果表明,錫林郭勒盟東南部防風固沙能力最強,年均防風固沙超過120 t·hm-2,主要包括正藍旗、正鑲白旗、蘇尼特左旗東南部等地區(qū),其面積占錫林郭勒盟總面積的9.91%。防風固沙能力較高區(qū)面積最大,占全盟面積的53.81%,集中分布在蘇尼特右旗、阿巴嘎旗以及蘇尼特左旗西北部等區(qū)域,多年平均防風固沙能力介于60～120 t·hm-2之間。防風固沙能力一般區(qū)的空間分布零散,其面積約占錫林郭勒盟總面積的22.45%,年均防風固沙能力為20～60 t·hm-2。另外,錫林郭勒盟分布有8.53%的防風固沙能力較低區(qū),主要位于東烏珠穆沁旗、西烏珠穆沁旗、太仆寺旗、多倫縣和鑲黃旗以及錫林浩特市等,年均防風固沙能力處于5～20 t·hm-2之間。而防風固沙能力低值區(qū)的年均防風固沙能力不足5 t·hm-2,其面積約占錫林郭勒盟總面積的5.30%,零散分布于阿巴嘎旗和蘇尼特左旗西北部、正鑲白旗東南部、太仆寺旗中部以及蘇尼特右旗西部等地區(qū)。因此,錫林郭勒盟防風固沙功能整體呈現由東南向西北與東北方向遞減的空間分布特征(圖4),主要原因是東北部地區(qū)風力侵蝕潛在風險較小,實際存在的植被防風固沙能力較低。

圖3 2000—2015年錫林郭勒盟土壤風蝕強度空間分布

表1 2000—2015年錫林郭勒盟防風固沙能力變化

Table 1 Changes of sand-fixing capacity in the Xilin Gol League from 2000 to 2015

防風固沙能力分區(qū)固沙能力(ASR)分級依據/(t·hm-2)分區(qū)面積占比/%防風固沙能力變化分區(qū)變化值(Δ)分級依據/(t·hm-2)變化區(qū)面積占比/% 高值區(qū)ASR≥1209.91明顯增加區(qū)Δ≥203.72 較高區(qū)60≤ASR<12053.81一般增加區(qū)5≤Δ<2028.91 一般區(qū)20≤ASR<6022.45無變化區(qū)-5≤Δ<551.82 較低區(qū)5≤ASR<208.53一般降低區(qū)-20≤Δ<-59.16 低值區(qū)ASR<55.30明顯降低區(qū)Δ<-206.39

3.3 防風固沙功能區(qū)域變化特征

從防風固沙功能變化區(qū)域來看,2000—2015年錫林郭勒盟有32.63%的區(qū)域防風固沙能力增高(表1)。其中,正藍旗、太仆寺旗和西烏珠穆沁旗等防風固沙功能明顯提升,單位面積防風固沙增加量超過20 t·hm-2,其面積占錫林郭勒盟總面積的3.72%;阿巴嘎旗、正鑲白旗、蘇尼特左旗東南部、東烏珠穆沁旗和多倫縣西南部防風固沙能力有所增加,增高幅度介于5～20 t·hm-2之間,其面積約占錫林郭勒盟總面積的28.91%。相比2000年,2015年錫林郭勒盟有15.55%的地域防風固沙功能下降。其中,蘇尼特右旗東部和蘇尼特左旗西北部防風固沙能力明顯降低,其面積占錫林郭勒盟總面積的6.39%,單位面積防風固沙能力降低幅度超過20 t·hm-2;東烏珠穆沁旗中部防風固沙能力變化為一般降低水平,其面積約占錫林郭勒盟總面積的9.16%。其余51.82%的地區(qū)防風固沙功能無明顯變化,主要分布在東烏珠穆沁旗、阿巴嘎旗和多倫縣東南部、錫林浩特市和鑲黃旗等(圖5)。可見,2000—2015年錫林郭勒盟防風固沙功能較穩(wěn)定,而東南部和中北部地區(qū)防風固沙能力明顯提升,西部地區(qū)防風固沙能力有明顯下降現象,需要重點加以關注。

3.4 防風固沙功能影響因素分析

區(qū)域防風固沙功能不僅受到植被覆蓋狀況影響,同時受降水量、土壤、地形和風力等其他因素影響。統(tǒng)計結果表明,2000—2015年錫林郭勒盟植被覆蓋度變動在40.12%～59.05%之間,且以每年0.32%的平均增速波動增加(圖6)。利用SPSS軟件進行Kendall非參數相關性檢驗,結果發(fā)現植被覆蓋度對區(qū)域防風固沙能力的提升具有顯著正面影響(Sig為0.034,r=0.531)。2000—2015年錫林郭勒盟降水量變動在196.01～389.44 mm之間,多年平均值為264.51 mm。然而,評估期間錫林郭勒盟降水量整體以4.94 mm·(15 a)-1的增速波動增加,且與防風固沙能力變化趨勢呈現較好的一致性(圖7)。經相關性檢驗,降水量與防風固沙能力也存在顯著正相關(Sig為0.037,r=0.551),說明該區(qū)域降水量對防風固沙能力的變化產生了積極影響。從防風固沙能力變化較大年份來看,2005年降水量與植被覆蓋度比2004年均明顯降低,造成防風固沙功能顯著下降;2012年降水量與植被覆蓋度大幅增加,明顯提高了防風固沙能力。因此,近年來降水量變化和生態(tài)工程建設均對區(qū)域防風固沙功能具有明顯影響。

圖4 錫林郭勒盟多年平均防風固沙能力空間分布

圖5 2000—2015年錫林郭勒盟防風固沙能力變化區(qū)

圖6 2000—2015年錫林郭勒盟防風固沙能力與植被覆蓋變化

圖7 2000—2015年錫林郭勒盟防風固沙能力與降水量變化

4 討論

基于錫林郭勒盟的植被、土壤及氣象數據,采用修正風蝕方程(RWEQ)和GIS空間分析技術,測算評估了2000—2015年錫林郭勒盟防風固沙功能變化及其區(qū)域差異。筆者研究發(fā)現,錫林郭勒盟多年平均防風固沙能力為74.48 t·hm-2,稍高于鞏國麗等[16]測算1990—2010年內蒙古典型草原區(qū)防風固沙能力50～70 t·hm-2的結果,原因是兩項研究評估階段不同,且近年來該區(qū)域防風固沙功能明顯增加。筆者研究中防風固沙能力值高于江凌等[23]估算的2000—2010年內蒙古生態(tài)系統(tǒng)年均防風固沙能力為48.80 t·hm-2的結論,這主要是因為錫林郭勒盟位于內蒙古風沙治理關鍵地帶,是發(fā)揮防風固沙功能的重要區(qū)域。此外,筆者研究發(fā)現錫林郭勒盟防風固沙能力空間變化特征與鞏國麗等[16]、王艷琦等[13]研究結果一致,植被覆蓋度變化與李慶旭等[3]、邵艷瑩等[26]研究結果一致。不過,佟斯琴等[27]、孫斌等[28]認為,1980—2010年間內蒙古降水量呈現波動下降趨勢,而筆者研究發(fā)現錫林郭勒盟降水量在2000—2015年呈現波動增加趨勢,原因在于評估時間段與研究區(qū)域的不同。

土壤風蝕強度易受地形、土壤、植被、氣候和人類活動的綜合影響[1-2]。其中,氣候因素主要表現為風速、降水和溫度等方面的影響[10,29],降水在直接影響土壤濕度和粘附能力的同時,間接調節(jié)植被生長狀況進而影響風蝕[30]。植被覆蓋度通過減緩風速和減少沉積物可有效保護表層土壤[31],因此筆者研究選取降水量和植被覆蓋度對研究區(qū)防風功能的影響進行分析。結果表明,錫林郭勒盟防風固沙功能的變化,不僅受到風沙源治理工程措施的影響,而且與風力、氣溫等氣象條件有關,這與ZHANG等[31]對內蒙古土壤風蝕動態(tài)的影響因素研究結果相一致。因此未來研究應重點關注防風固沙功能與工程措施以及氣象要素波動的關系解析。

PI等[32]應用REWQ經驗模型對比分析了中國和美國干旱半干旱地區(qū)的防風固沙能力,結果表明作為農田風蝕模型的RWEQ具有一定誤差,筆者將其應用于以草原為主的研究區(qū)也存在一定局限性。因此,雖然RWEQ模型方便定量評估土壤侵蝕模數,但其計算參數均來源于美國大平原的統(tǒng)計值,屬于缺乏一定理論與物理過程基礎的經驗模型[33],且具有很強地域性,具體應用該模型時要對計算參數進行修正。鞏國麗等[34]對RWEQ中土壤結皮和可蝕性因子進行改進,并將其應用到我國北方風沙土區(qū)。不過,遲文峰等[35]利用內蒙古高原137Cs示蹤技術監(jiān)測成果檢驗RWEQ模型反演結果發(fā)現具有較好擬合性(R2=0.83,P<0.01)。筆者研究在綜合考慮與借鑒前人對公式與參數修正結果的基礎上,為滿足模型對數據空間和時間分辨率的要求,對部分輸入參數進行了插值處理,可能導致評估結果存在一定程度的誤差。此外,RWEQ模型在風蝕因子分類與相互影響方面存在不足,以至于風蝕影響因子僅是特定區(qū)域的經驗表達,不具有普適性的風蝕動力學理論基礎[21]。因此,構建并應用具有理論基礎與廣泛適用性的土壤風蝕模型更能準確反映區(qū)域防風固沙功能變化狀況。

5 結論

研究表明,2000—2015年錫林郭勒盟防風固沙功能明顯提升,年均防風固沙量為14.56億t,防風固沙能力達74.78 t·hm-2,且整體呈現由東南向西北與東北方向遞減趨勢。相比2000年,2015年錫林郭勒盟東南部防風固沙能力增加明顯,全盟有32.63%的區(qū)域防風固沙能力增高,51.82%的區(qū)域防風固沙能力穩(wěn)定,另有西部地區(qū)15.55%的區(qū)域防風固沙能力降低,是區(qū)域防風固沙功能需要重點關注區(qū)域。此外,錫林郭勒盟防風固沙功能變化與植被覆蓋度和降水量變化呈顯著相關,說明近年來降水變化和生態(tài)工程建設均對區(qū)域生態(tài)功能的提升有積極作用。風沙源治理工程的推進實施應綜合考慮氣候變化、生態(tài)工程和人類活動的復合影響,利用發(fā)揮氣候變化帶來的正面效應,注重生態(tài)工程布局建設的區(qū)域適宜性,并對生態(tài)系統(tǒng)進行科學的保護與修復。