锡林郭勒盟风沙源治理区防风固沙功能变化评估

张 彪,王 爽

(1.中国科学院地理科学与资源研究所,北京 100101;2.中国科学院大学资源与环境学院,北京 100049)

土地沙化是全球面临的重大生态环境问题,其与地形、土壤、植被、气候和人类活动等因素有关[1-2]。近年来我国土地退化态势较严重。根据2015年发布的《中国荒漠化和沙化状况公报》,截至2014年,我国仍有荒漠化土地261.16万km2和沙化土地172.12万km2。为遏制土地沙化趋势,我国部署实施了退耕还林还草、三北防护林体系、天然林保护和京津风沙源治理等重大生态建设工程。由于生态建设工程工期长，地域广,加上生态系统恢复漫长复杂,及时监测、评估工程区生态状况及功能变化尤为重要[3]。京津风沙源治理工程是为遏制北京及周边地区土地沙化趋势、改善京津地区大气环境质量而实施的一项重大生态工程,建设范围西起内蒙古达尔罕茂明安联合旗,东至河北省平泉县,南至山西省代县,北至内蒙古东乌珠穆沁旗,涉及国土面积45.8万km2。自2000年京津风沙源治理工程实施以来,工程区植被恢复与土壤风蚀趋势受到重点关注。李庆旭等[3]基于MOD13Q1数据测算2000—2015年京津风沙源区植被覆盖度变化范围在35%～45%之间,且以年均0.4%速率波动增加。石莎等[4]利用传统生态学方法对京津风沙源治理工程区2001—2005年植被恢复效果进行野外调查,发现植被覆盖度整体上升。严恩萍等[5]采用MODIS02B数据源分析了重点治理区(涉及山西和内蒙古共48个县域)植被覆盖变化特征,表明2000—2012年治理区植被覆盖总体呈上升趋势。卫洁等[6]利用MODIS数据监测发现,2000—2014年京津风沙源区84.51%的地区植被覆盖呈增加趋势。生态工程区植被状况变化必然引起区域防风固沙功能的相应变化。高尚玉等[7]调查评估发现,2005年京津风沙源区土壤风蚀量比2001年减少1.96亿t。吴丹等[8]综合评估表明,2000—2010年京津风沙源治理工程区土壤风蚀以微度和轻度为主,且风蚀量总体呈下降态势。迟文峰等[9]以遥感手段监测发现,1990—2015年间工程区土壤风蚀模数整体下降,且在生态工程实施后改善趋势明显。1990—2015年内蒙古49.23%的区域土壤风蚀强度下降,但仍有7.11%的区域风蚀强度增加[10]。然而,相比于工程区植被覆盖状况的监测,区域防风固沙功能的动态变化研究仍相对薄弱。

锡林郭勒盟位于京津风沙源治理工程区西北部,是工程区内面积最大的地级行政区。锡林郭勒盟中部是东西长约400 km的浑善达克沙地,其西缘为固定、半固定沙丘广布的嘎亥额勒苏沙地,其生态治理效果成为评估京津风沙源治理工程成效的关键。巩国丽等[11]运用修正风蚀方程(RWEQ)评估发现,20世纪90年代以来锡林郭勒盟土壤风蚀强度总体减弱。丁文广等[12]利用Landsat卫星影像分析发现,1975—2014年锡林郭勒盟沙漠化面积先增后减,沙漠化重心先东移后西返。王艳琦等[13]采用Landsat-8 OLI遥感影像分析认为,2000—2015年锡林郭勒盟沙漠化程度总体呈现逆转趋势。此外,浑善达克沙漠化防治生态功能区的防风固沙功能受到关注[14]。但是,以往研究侧重锡林郭勒盟土壤风蚀状况的变化趋势及驱动力分析,而对区域防风固沙功能的动态变化及其空间差异揭示不足。

随着京津风沙源治理二期工程(2013—2022年)的实施,分区施策与精准修复已成为生态治理工程的重要内容[15]。为此,以锡林郭勒盟为例,采用REWQ模型与GIS空间分析技术,评估分析2000—2015年锡林郭勒盟防风固沙功能的时空变化及其区域差异,并探讨植被覆盖度和降水量对区域防风固沙功能的影响,可为风沙源治理工程的分区施策治理与成效监测评估提供参考依据。

1 研究区概况

锡林郭勒盟位于京津风沙源治理工程区西北部、内蒙古中部(图1),是典型农牧交错地区,也是北京和华北地区的重要生态屏障[16]。锡林郭勒盟土地总面积为202 580 km2,现辖9旗(东乌珠穆沁旗、西乌珠穆沁旗、阿巴嘎旗、苏尼特左旗、苏尼特右旗、正蓝旗、正镶白旗、镶黄旗、太仆寺旗)、2市(锡林浩特市、二连浩特市)和1县(多伦县)。

图1 锡林郭勒盟位置与组成

锡林郭勒盟地形平坦开阔,地势自西南向东北方向倾斜。海拔高度处于800～1 800 m之间[16],海拔最高处位于太仆寺旗东部,最低处位于东乌珠穆沁旗和西乌珠穆沁旗。锡林郭勒盟气候属于中温带干旱、半干早大陆性季风气候区,年均气温为1～2 ℃,年平均降水量为200～350 mm,且集中在6—8月[17]。锡林郭勒盟主导风向为西北风,全年大风日数达50～80 d,其中3—5月大风日数占全年的40%～50%[18]。锡林郭勒盟地带性植被为草原,其面积约占草地总面积的71.88 %,类型包括典型草原、草甸草原和荒漠草原3大亚型[19]。据锡林郭勒盟2018年统计年鉴,2017年末全盟常住人口为105.16万人,自然增长率为4.1%。其中,城镇常住人口为68.49万人,占总人口比例为65.13%;农村牧区常住人口为36.67万人,占总人口比例为34.87%[20]。

2 材料与方法

2.1 RWEQ模型的构建

土壤风蚀是风力作用导致表土物质脱离原空间位置的过程[21]。风蚀模型是评估土壤风蚀状况的主要技术手段,并已发展演化出通用风蚀方程(WEQ)、德克萨斯侵蚀分析模型(TEAM)和风蚀预报系统(WEPS)等多种模型[22]。其中,修正风蚀方程(RWEQ)综合反映了气候、植被、土壤性质和地表粗糙度等自然因素的影响,已被广泛应用于我国土壤风蚀状况评估[9-11,14]。防风固沙功能是在地表无植被状况下土壤风蚀量与植被覆盖条件下土壤风蚀量的差值[10,23],采用RWEQ模型[22]定量评估锡林郭勒盟潜在风蚀量(LS,p)与实际风蚀量(LS,r),以两者之差来表示区域防风固沙量(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为研究区防风固沙量,t·a-1;A为研究区面积,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)

土壤可蚀性受土壤颗粒粒径以及有机质、黏土和碳酸钙等物质含量的影响,表层坚硬结皮也能有效防止风蚀的发生[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)。由于区域尺度评估中耕作产生的随机糙度难以获取,采用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/),在锡林郭勒盟辖区分布有10个气象站,选取2000—2015年共16 a气象观测数据,采用月均风速、降水、气温和日照时数等数据插值得到风力因子和土壤湿度因子。雪盖因子利用中国西部环境与生态科学数据中心(http:∥westdc.westgis.ac.cn)的中国雪深长时间序列数据集计算。土壤数据来源于中国西部环境与生态科学数据中心提供的1∶100万土壤图及所附的土壤属性表和空间数据,并分别采用锡林郭勒盟不同土壤类型的相应物质含量估算可蚀性因子与结皮因子。NDVI数据来自美国地球资源观测系统数据中心的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%的变化率波动增加,直观表现为区域土壤风蚀强度的整体下降(图3)。但区域防风固沙能力增加幅度稍低于防风固沙量,主要原因是防风固沙量受到防风固沙能力与区域面积的双重影响。

图2 2000—2015年锡林郭勒盟防风固沙功能变化

3.2 防风固沙功能空间差异特征

为揭示锡林郭勒盟防风固沙功能空间差异,依据全盟固沙能力值域及风力侵蚀强度分级[25]设置相应防风固沙功能分级标准(表1)。

结果表明,锡林郭勒盟东南部防风固沙能力最强,年均防风固沙超过120 t·hm-2,主要包括正蓝旗、正镶白旗、苏尼特左旗东南部等地区,其面积占锡林郭勒盟总面积的9.91%。防风固沙能力较高区面积最大,占全盟面积的53.81%,集中分布在苏尼特右旗、阿巴嘎旗以及苏尼特左旗西北部等区域,多年平均防风固沙能力介于60～120 t·hm-2之间。防风固沙能力一般区的空间分布零散,其面积约占锡林郭勒盟总面积的22.45%,年均防风固沙能力为20～60 t·hm-2。另外,锡林郭勒盟分布有8.53%的防风固沙能力较低区,主要位于东乌珠穆沁旗、西乌珠穆沁旗、太仆寺旗、多伦县和镶黄旗以及锡林浩特市等,年均防风固沙能力处于5～20 t·hm-2之间。而防风固沙能力低值区的年均防风固沙能力不足5 t·hm-2,其面积约占锡林郭勒盟总面积的5.30%,零散分布于阿巴嘎旗和苏尼特左旗西北部、正镶白旗东南部、太仆寺旗中部以及苏尼特右旗西部等地区。因此,锡林郭勒盟防风固沙功能整体呈现由东南向西北与东北方向递减的空间分布特征(图4),主要原因是东北部地区风力侵蚀潜在风险较小,实际存在的植被防风固沙能力较低。

图3 2000—2015年锡林郭勒盟土壤风蚀强度空间分布

表1 2000—2015年锡林郭勒盟防风固沙能力变化

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

防风固沙能力分区固沙能力(ASR)分级依据/(t·hm-2)分区面积占比/%防风固沙能力变化分区变化值(Δ)分级依据/(t·hm-2)变化区面积占比/% 高值区ASR≥1209.91明显增加区Δ≥203.72 较高区60≤ASR<12053.81一般增加区5≤Δ<2028.91 一般区20≤ASR<6022.45无变化区-5≤Δ<551.82 较低区5≤ASR<208.53一般降低区-20≤Δ<-59.16 低值区ASR<55.30明显降低区Δ<-206.39

3.3 防风固沙功能区域变化特征

从防风固沙功能变化区域来看,2000—2015年锡林郭勒盟有32.63%的区域防风固沙能力增高(表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%的地区防风固沙功能无明显变化,主要分布在东乌珠穆沁旗、阿巴嘎旗和多伦县东南部、锡林浩特市和镶黄旗等(图5)。可见,2000—2015年锡林郭勒盟防风固沙功能较稳定,而东南部和中北部地区防风固沙能力明显提升,西部地区防风固沙能力有明显下降现象,需要重点加以关注。

3.4 防风固沙功能影响因素分析

区域防风固沙功能不仅受到植被覆盖状况影响,同时受降水量、土壤、地形和风力等其他因素影响。统计结果表明,2000—2015年锡林郭勒盟植被覆盖度变动在40.12%～59.05%之间,且以每年0.32%的平均增速波动增加(图6)。利用SPSS软件进行Kendall非参数相关性检验,结果发现植被覆盖度对区域防风固沙能力的提升具有显著正面影响(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),说明该区域降水量对防风固沙能力的变化产生了积极影响。从防风固沙能力变化较大年份来看,2005年降水量与植被覆盖度比2004年均明显降低,造成防风固沙功能显著下降;2012年降水量与植被覆盖度大幅增加,明显提高了防风固沙能力。因此,近年来降水量变化和生态工程建设均对区域防风固沙功能具有明显影响。

图4 锡林郭勒盟多年平均防风固沙能力空间分布

图5 2000—2015年锡林郭勒盟防风固沙能力变化区

图6 2000—2015年锡林郭勒盟防风固沙能力与植被覆盖变化

图7 2000—2015年锡林郭勒盟防风固沙能力与降水量变化

4 讨论

基于锡林郭勒盟的植被、土壤及气象数据,采用修正风蚀方程(RWEQ)和GIS空间分析技术,测算评估了2000—2015年锡林郭勒盟防风固沙功能变化及其区域差异。笔者研究发现,锡林郭勒盟多年平均防风固沙能力为74.48 t·hm-2,稍高于巩国丽等[16]测算1990—2010年内蒙古典型草原区防风固沙能力50～70 t·hm-2的结果,原因是两项研究评估阶段不同,且近年来该区域防风固沙功能明显增加。笔者研究中防风固沙能力值高于江凌等[23]估算的2000—2010年内蒙古生态系统年均防风固沙能力为48.80 t·hm-2的结论,这主要是因为锡林郭勒盟位于内蒙古风沙治理关键地带,是发挥防风固沙功能的重要区域。此外,笔者研究发现锡林郭勒盟防风固沙能力空间变化特征与巩国丽等[16]、王艳琦等[13]研究结果一致,植被覆盖度变化与李庆旭等[3]、邵艳莹等[26]研究结果一致。不过,佟斯琴等[27]、孙斌等[28]认为,1980—2010年间内蒙古降水量呈现波动下降趋势,而笔者研究发现锡林郭勒盟降水量在2000—2015年呈现波动增加趋势,原因在于评估时间段与研究区域的不同。

土壤风蚀强度易受地形、土壤、植被、气候和人类活动的综合影响[1-2]。其中,气候因素主要表现为风速、降水和温度等方面的影响[10,29],降水在直接影响土壤湿度和粘附能力的同时,间接调节植被生长状况进而影响风蚀[30]。植被覆盖度通过减缓风速和减少沉积物可有效保护表层土壤[31],因此笔者研究选取降水量和植被覆盖度对研究区防风功能的影响进行分析。结果表明,锡林郭勒盟防风固沙功能的变化,不仅受到风沙源治理工程措施的影响,而且与风力、气温等气象条件有关,这与ZHANG等[31]对内蒙古土壤风蚀动态的影响因素研究结果相一致。因此未来研究应重点关注防风固沙功能与工程措施以及气象要素波动的关系解析。

PI等[32]应用REWQ经验模型对比分析了中国和美国干旱半干旱地区的防风固沙能力,结果表明作为农田风蚀模型的RWEQ具有一定误差,笔者将其应用于以草原为主的研究区也存在一定局限性。因此,虽然RWEQ模型方便定量评估土壤侵蚀模数,但其计算参数均来源于美国大平原的统计值,属于缺乏一定理论与物理过程基础的经验模型[33],且具有很强地域性,具体应用该模型时要对计算参数进行修正。巩国丽等[34]对RWEQ中土壤结皮和可蚀性因子进行改进,并将其应用到我国北方风沙土区。不过,迟文峰等[35]利用内蒙古高原137Cs示踪技术监测成果检验RWEQ模型反演结果发现具有较好拟合性(R2=0.83,P<0.01)。笔者研究在综合考虑与借鉴前人对公式与参数修正结果的基础上,为满足模型对数据空间和时间分辨率的要求,对部分输入参数进行了插值处理,可能导致评估结果存在一定程度的误差。此外,RWEQ模型在风蚀因子分类与相互影响方面存在不足,以至于风蚀影响因子仅是特定区域的经验表达,不具有普适性的风蚀动力学理论基础[21]。因此,构建并应用具有理论基础与广泛适用性的土壤风蚀模型更能准确反映区域防风固沙功能变化状况。

5 结论

研究表明,2000—2015年锡林郭勒盟防风固沙功能明显提升,年均防风固沙量为14.56亿t,防风固沙能力达74.78 t·hm-2,且整体呈现由东南向西北与东北方向递减趋势。相比2000年,2015年锡林郭勒盟东南部防风固沙能力增加明显,全盟有32.63%的区域防风固沙能力增高,51.82%的区域防风固沙能力稳定,另有西部地区15.55%的区域防风固沙能力降低,是区域防风固沙功能需要重点关注区域。此外,锡林郭勒盟防风固沙功能变化与植被覆盖度和降水量变化呈显著相关,说明近年来降水变化和生态工程建设均对区域生态功能的提升有积极作用。风沙源治理工程的推进实施应综合考虑气候变化、生态工程和人类活动的复合影响,利用发挥气候变化带来的正面效应,注重生态工程布局建设的区域适宜性,并对生态系统进行科学的保护与修复。