呼伦贝尔森林-草原生态交错带土壤风蚀量时空变化及驱动力分析

胡梦甜,张 慧,2①,乔亚军,刘 坤,王 智②,徐网谷

(1.生态环境部南京环境科学研究所,江苏 南京 210042;2.南京信息工程大学江苏省大气环境与装备技术协同创新中心/ 地理科学学院,江苏 南京 210044)

风蚀对地表土壤的大量搬运和堆积,是导致干旱半干旱区土地沙化和沙地化进程最重要、最直接的作用过程之一[1-2]。土壤风蚀研究方法主要有野外调查观测、风洞模拟、元素示踪法、遥感和GIS等[3],但这些方法无法满足大范围区域风蚀量快速动态估算的需求[4]。为了定量模拟大范围区域风蚀量,自20世纪60年代起,国外先后开发了大量的风蚀模型,主要包括风蚀方程(wind erosion equation,WEQ)[5]、修正风蚀方程(revised wind erosion equation,RWEQ)[6]、德克萨斯模型(taxas erosion analysis model,TEAM)、风蚀评价模型(wind erosion assessment model,WEAM)[7]、风蚀预报系统(wind erosion prediction system,WEPS)[8]等。其中,WEPS系统较为先进,但该系统建模过程复杂,数据要求繁杂,且缺乏相关参数的本地化工作,故该模型在国内应用较少[9]。相比先进的WEPS系统,RWEQ模型因其参数较易获取、操作更为简单的特点,在我国不少地区的土壤风蚀和防风固沙功能计算与评估中得到了广泛应用,如在内蒙古自治区锡林郭勒盟、青海省浑善达克沙地、黑河流域[10-11]等地都开展了相关工作。该模型也被作为《生态保护红线划定指南》和《全国生态状况评估技术规范》中防风固沙功能重要性的评估方法。

对于风蚀量变化的驱动因素,学者大多采用相关性分析、主成分或聚类分析等传统方法进行研究[12],缺少从地理分异的角度对风蚀量变化的定量归因和驱动因素的空间差异性研究。地理探测器是一种强大的、能直接量化驱动因素及其交互作用影响的方法,它不必遵循传统统计方法的假设,且不涉及复杂的参数设置过程[13]。近年来,地理探测器在土地利用[14]、生态服务功能[15]等地理现象的空间驱动力分析领域应用广泛。

呼伦贝尔森林-草原生态交错带是我国生态系统结构保存完整、健康状况良好的林草交错带,是我国北方地区重要的生态屏障[16]。研究区冬春季寒冷多风,该区的草原是以栗钙土、风沙土为主的干草原,易形成风蚀[17],因此探究该地区土壤风蚀量的变化及其驱动因素对维护区域生态安全具有重要意义。近年来,虽然有学者研究了呼伦贝尔森林-草原生态交错带土壤风蚀,但主要集中于区域的生态效益评估[18-19],鲜有对区域风蚀量变化的驱动力开展相关研究。笔者利用RWEQ风蚀修正模型研究呼伦贝尔森林-草原生态交错带2000、2010和2018年的土壤风蚀量时空分布变化,并识别其驱动因素,以期为该区域科学防治土壤风蚀、遏制区域沙化趋势提供科学支撑。

1 研究区概况

呼伦贝尔森林-草原生态交错带处于大兴安岭西麓山地向亚洲中部蒙古高原东北部过渡的区域,东北区域为林区,海拔700～1 700 m,由东北向西南依次为农田、森林草原、草甸草原和干旱草原,海拔550～1 000 m。研究区位于内蒙古自治区呼伦贝尔市中部,地理位置处于北纬46°10′～53°26′,东经117°33′～122°55′ 之间,行政区域涉及呼伦贝尔市的额尔古纳市、根河市、牙克石市、陈巴尔虎旗、海拉尔区、鄂温克自治旗和新巴尔虎左旗共7个旗市(图1)。该地区的经济活动主要包括种植业和畜牧业,其中种植业以冬小麦、油菜和苜蓿等种植为主;畜牧业以养羊、牛为主。研究区处在温带-寒温带气候区,气候较干燥,多大风,研究区年平均气温在 -2.2～2.4 ℃之间,年降水量为290～450 mm[16],年均蒸发量590～960 mm,年日照时数为2 600～2 800 h,年均风速1.8～2.45 m·s-1。研究区土地利用类型以森林、草地和湿地为主,总占比超过90%。其中,森林集中分布在研究区的东北部山区;草地主要分布在西南部;农田集中分布在林草交错带,主要位于额尔古纳市、呼伦贝尔市、牙克石市区域内;湿地主要分布在额尔古纳河、海拉尔河、根河、辉河区域。

该图基于审图号为蒙S(2020)027号的标准地图制作。

2 数据来源与研究方法

2.1 数据来源

研究区采用的遥感数据源为2000和2010年的Landsat 5 TM卫星影像(空间分辨率30 m)、2018年的Landsat 8 TM 卫星影像(空间分辨率30 m)(https:∥glovis.usgs.gov/)。依据《全国30米分辨率土地利用分类系统》并结合研究区实际情况,考虑到大兴安岭地区在1987年发生过特大火灾,2000年过火林区的植被还未恢复,因此在土地利用分类中增加火烧迹地类型,将土地利用类型分为森林、草地、农田、湿地、城镇用地、沙地和火烧迹地7类。采用人工目视解译方法分类,精度达95%以上,土地利用变化检测总体精度达85%以上。气象数据包括月平均风速,数据来源于国家青藏高原科学数据中心(http:∥data.tpdc.ac.cn/zh-hans/data),空间分辨率为1 km。雪盖因子采用中国雪深长时间序列数据集(http:∥data.tpdc.ac.cn/zh-hans/data/df40346a-0202-4ed2-bb07-b65dfcda9368/);土壤特性因子数据来自世界土壤数据库(HWSD,http:∥webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database),包括土壤粗砂、细砂、黏粒及有机质含量等理化性质;NDVI数据采用资源环境科学与数据中心发布的数据产品(http:∥www.resdc.cn/),空间分辨率为1 km;土壤湿度数据采用干旱植被指数法计算得到,其中地表温度采用NASA发布的MOD11A2(https:∥ladsweb.modaps.eosdis.nasa.gov)数据;地表糙度因子采用GDEM DEM 30 m数字高程数据(http:∥www.gscloud.cn/)计算得到。将所有因子统一为250 m分辨率的栅格数据,保证栅格计算过程准确。

2.2 风蚀量计算模型

风蚀量计算采用RWEQ模型, RWEQ作为美国农田风蚀模型,其计算参数均来源于美国本土[20]。迟文峰等[21]在内蒙古高原利用137Cs示踪技术检验RWEQ模型的模拟效果,发现该模型的风蚀量计算结果拟合性较好(R2=0.83,P<0.01)。《生态保护红线划定指南》[22]将该模型作为风蚀量和防风固沙功能的推荐评估方法,RWEQ方法计算公式为

(1)

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

(2)

Qmax=109.8×FW×FE×FSC×K′×C。

(3)

式(1)～(3)中,SL为实际土壤侵蚀量,t·km-2·a-1;S为区域防风固沙系数;Qmax为风沙滞留量,kg·m-1;Z为最大风蚀出现距离,m;FW为气候因子,kg·m-1;FE为土壤可蚀性因子;FSC为土壤结皮因子;K′为地表粗糙度因子;C为植被因子。

气候因子(FW)计算公式为

(4)

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

(5)

(6)

ρ=348×(1.013-0.118 3×LE+0.004 8×LE2)/T。

(7)

式(4)～(7)中,fW为风力因子,m3·s-3;ρ为空气密度,kg·m-3;g为重力加速度,取9.8 m·s-2;WS为各月土壤湿度因子,表征土壤湿度抑制土壤风蚀效果的大小,土壤湿度越大,土壤湿度因子越小,越不易起沙,该因子以往大多采用气象站点数据插值计算得到,由于气象站点分布不均,插值结果往往有“牛眼”存在,因此采用干旱植被指数法表征土壤湿度因子[23-24];DS为雪盖因子;u1为起沙风速,参照江凌[25]的计算方法,取5 m·s-1;j为1个月内日平均风速≥u1的天数,j=1,2,…,m;uj为第j天的日平均风速,m·s-1;TLS,i为i评估区域的地表温度,℃;TLS,dry为评估区域NDVI对应的最高地表温度,℃,即干边;TLS,wet为评估区域NDVI对应的最低地表温度,℃,即湿边;LE为海拔高度,km;T为绝对温度,K,即在各月平均气温数据t的基础上加常数273.15。

土壤可蚀性因子(FE)和土壤结皮因子(FSC)计算公式为

(8)

(9)

式(8)～(9)中,wsa为土壤粗砂含量,%;wsi为土壤粉砂含量,%;wcl为土壤黏粒含量,%;wOM为土壤有机质含量,%;wCa为碳酸钙含量,%,此次计算未予考虑,其值取0。

植被覆盖因子(C)计算公式为

C=e-CS×ai。

(10)

式(10)中,CS为植被覆盖度,%,由12个月的NDVI指数计算得到年均植被覆盖度;ai为不同植被类型的系数,因不同植被类型的土壤风蚀效果不同,森林、草地、农田、沙地(包括裸地和沙地)分别取值0.153 5、0.115 1、0.043 8、0.071 3[25]。

地表糙度因子(K′)计算公式为

K′=e1.86Kr-0.127Crr-2.41Kr0.934,

(11)

(12)

式(11)～(12)中,Kr为地形粗糙度长度,m;Crr为随机糙度,在区域尺度的计算中可以忽略不计;L为地势起伏参数,m;△H为距离L范围内的海拔高程差,m。

2.3 驱动因素选择

土壤风蚀的产生受到土壤、地形等自然本底因素的制约,同时气温、降水、风速等气候因素以及人类活动也会对其产生影响[10]。根据研究区的实际情况,从气候因素和人类活动方面选取5个影响因子,包括降水变化(X1)、气温变化(X2)、植被覆盖度变化(X3)、风速变化(X4)和土地利用类型(X5),作为探测研究区土壤风蚀变化量的驱动因素。以研究区范围为基础创建渔网,参照各影响因子空间分辨率将渔网大小设置为1 000 m×1 000 m。采用ArcGIS软件Spatial Analysis Tool工具包中的Zonal Statistic工具,将各影响因子按照平均值统计到各个渔网。因为地理探测器中自变量为类型量,对于顺序量需要进行离散化[13]。在R语言环境下,采用等间隔法、分位数法、自然断点法和标准差法比较降水变化、气温变化、植被覆盖度变化、风速变化4个因子的离散化效果,采用自然断点法将植被覆盖度变化、降水变化、气温变化、风速变化分别分为8、6、9、9类。土地利用之间转化类型达38类,因此将土地利用变化因子分为38类。为深入探究呼伦贝尔森林-草原生态交错带土壤风蚀量时空变化驱动机制,针对2000—2018年土壤风蚀量变化的主要影响因素开展研究。

2.4 地理探测器模型

采用地理探测器模型进行土壤风蚀量侵蚀变化的驱动力分析。地理探测器模型主要包括因子探测器和交互探测器[26]。因子探测器用于分析自变量对因变量的解释程度,解释力的强弱通过比较q值的大小反映[13],q取值在0～1之间,q值越大,表示该影响因子对土壤风蚀变化量的影响越大。交互作用探测器可探索2个自变量的联合效应是否会增加、减少对因变量的解释力[27]。通过比较因子单独作用时的q值〔q(x1)和q(x2)〕与交互作用时的q值〔q(x1∩x2)〕,对2个因子之间的关系进行界定。q(x1∩x2)>q(x1)+q(x2)、q(x1∩x2)同时大于q(x1)和q(x2)、q(x1∩x2)处于q(x1)和q(x2)之间、q(x1∩x2)同时小于q(x1)和q(x2)、q(x1∩x2)=q(x1)+q(x2)分别表征非线性增强、双因子增强、单因子非线性减弱、非线性减弱、独立。

3 结果与分析

3.1 呼伦贝尔森林-草原生态交错带土壤风蚀量及强度的时空分布变化

从风蚀量的年际变化来看,2000、2010和2018年研究区土壤风蚀总量分别为9.74×107、1.33×108、8.51×107t,总体呈现先上升后下降的趋势,2000—2018年土壤风蚀总量减少1.23×107t,平均年降速为68.33万t,年均下降率为0.70%(图1)。

评估期内,研究区单位面积风蚀量在7.03～10.88 t·hm-2之间波动,单位面积风蚀量的年均下降率为0.54%,研究区单位面积风蚀量的下降幅度低于土壤风蚀总量的下降幅度。

从研究区土壤风蚀量及强度的空间分布来看,土壤风蚀量整体呈现由东北部林区向西南部草原区逐渐增加的特征(图1)。依据《土壤侵蚀分类分级标准》将研究区土壤风蚀强度分为5级(剧烈、极强烈、强烈、中度、轻度和微度)(表1)。结果表明,2018年剧烈风蚀区域面积占研究区总面积的7.61%,该区域年均单位面积风蚀量达378.91 t·hm-2,主要包括新巴尔虎左旗西部、陈巴尔虎旗南部和鄂温克族自治旗西部区域;土壤极强烈风蚀区域面积占研究区总面积的5.86%,该区域年均单位面积风蚀量达105.91 t·hm-2,集中分布在新巴尔虎左旗中部,零散分布在鄂温克族自治旗;土壤强烈和中度风蚀区域面积占比相近,分别为4.13%和3.97%,分布于新巴尔虎左旗、陈巴尔虎旗和鄂温克族自治旗交界处;研究区有13.86%的土壤属轻度风蚀区,主要位于林草过渡带;土壤微度风蚀区面积占研究区总面积的64.57%,主要分布在林区。

表1 2000—2018年呼伦贝尔森林-草原生态交错带土壤风蚀变化

从研究区风蚀量及强度空间变化来看,2000—2018年,研究区有90.86%的区域土壤风蚀强度未发生变化;有9.04%的区域土壤风蚀强度减轻(表1),主要分布在西南部的新巴尔虎左旗、陈巴尔虎旗和鄂温克族自治旗,单位面积风蚀量减少了27 t·hm-2;有0.10%的区域土壤风蚀强度加重,主要分布在新巴尔虎器西北角和鄂温克族自治旗正北角,单位面积风蚀量增加了20 t·hm-2。可见,2000—2018年研究区土壤风蚀程度整体较为稳定,部分地区加重(图2),其中西南部地区土壤风蚀程度减轻,而新巴尔虎器西北角和鄂温克族自治旗正北角部分地区土壤风蚀强度恶化。

该图基于审图号为蒙S(2020)027号的标准地图制作。

3.2 土壤风蚀量变化的驱动力分析

依据地理探测器分析得出,单个影响因子变化量对土壤风蚀量变化量的解释力排序为土地利用变化>降水变化>风速变化>植被覆盖度变化>气温变化。整体来看,各因子变化量对土壤侵蚀量的解释力q值均很小,最大值亦不超过0.1,表明单一因子变化量对土壤风蚀量变化驱动作用有限。与风速变化和气温变化因子相比,降水变化对土壤风蚀量变化的影响力较高,说明降水量增加对区域风蚀量减少发挥着重要作用[11]。进一步探究各因子变化量的交互作用,发现任意2种因子交互作用的解释力高于单个因子(表2),且降水变化与其他因子变化均呈双因子非线性增强作用。

表2 各驱动因子对土壤风蚀量变化量的解释力(q值)

其中,降水变化和土地利用变化的协同作用对土壤风蚀量变化的解释力最大,解释力q值达0.22。同时,植被覆盖度变化与其他因子变化也均呈双因子非线性增强作用,且植被覆盖度变化协同降水变化解释力最强,q值为0.14;植被覆盖度变化协同土地利用变化次之,q值为0.13。

3.3 土壤风蚀量变化的人为影响因素分析

利用研究区2000—2018年土地利用现状图和土地利用转换图(图3～4),定量分析土地利用变化对风蚀量的影响。2000—2018年大兴安岭林草交错带各土地利用类型发生了复杂的相互转换,农田转出面积变大,草地面积增加,湿地面积萎缩,退耕还草和城市化是农田的主要转出方向。2000—2010年,退耕还草是农田面积减少、草地面积增加的主要转换方式,退耕还草面积达到201.71 km2,以鄂温克族自治旗林草交错带内的农田转换为主。2010年以后,城镇用地增加成为农田主要的转出方向,期间城镇用地增加占用农田40.35 km2,城镇用地占用农田扩张趋势明显。2000—2018年间,湿地面积净减少117.66 km2,尤其是2000—2010年,湿地大面积退化为草地,多发生在新巴尔虎旗境内的呼伦湖、海拉尔流域(表3)。2000—2018年,土地利用类型转换分别造成了24.32 t风蚀量的增加和12.31×106t风蚀量的减少。

表3 研究区2000—2018年土地利用转移面积表

该图基于审图号为蒙S(2020)027号的标准地图制作。

地理探测器研究表明,土地利用变化协同植被覆盖度变化将增大对土壤风蚀量变化的影响,因此将土地利用变化与植被覆盖度变化协同分析土壤风蚀量变化的人为影响因素。将植被覆盖度按照马志勇等[28]提出的植被覆盖度分级标准,按0～0.45、>0.45～0.75和>0.75～1.00分别将研究区植被覆盖度分为低、中、高3类,研究不同土地利用类型、不同植被盖度之间的转化对土壤风蚀量的影响(植被覆盖度分类中湿地为沼泽湿地,不包含河流和湖泊)。

表4列出了居前16位的不同植被覆盖度的土地利用类型下土壤风蚀量变化,研究区因为草地覆盖度增加减少了8.53×106t的土壤风蚀量,主要分布在陈巴尔虎旗、新巴尔虎左旗和鄂温克族自治旗大部分地区;因草地覆盖度降低增加了24.32 t土壤风蚀量,主要分布在鄂温克族自治旗西部巴彦乌拉苏木地区。沙地转化为草地、湿地、森林后土壤风蚀量减少了2.38×106t,主要分布在新巴尔虎旗东部、陈巴尔虎旗中部、鄂温克族自治旗中西部这3条沙带区[29];湿地覆盖度提高后,土壤风蚀量减少了7.34×105t,主要分布在额尔古纳河、海拉尔河等沼泽湿地区域。

表4 研究区不同植被覆盖度的土地利用类型下土壤风蚀量变化

该图基于审图号为蒙S(2020)027号的标准地图制作。

草地覆盖度增加、沙化土地封育(沙地质量改善、沙地转草地、沙地转湿地、沙地转森林)、生态退耕(农田转草地、林地和湿地)、天然林保护(森林覆盖度增加)对风蚀量减少的贡献率分别为69.32%、19.37%、0.06%和1.81%,占土壤风蚀减少总量的90.56%。其中,草地覆盖度增加导致的风蚀量降低最为明显,其次是沙化土地封育,虽然草地覆盖度增加、沙化土地封育与研究区降水增加有一定关系,但是这两者与研究区的围封禁牧、沙化土地封育政策关系更加密切,这些保护政策增加了草地生物量、盖度及高度,减少了草地的风蚀作用,提高了草地的保水能力[30]。因此,今后应继续采取围封禁牧、休牧、轮牧、改良牧草场等措施,增加其覆盖度,可以有效减少该区域的风蚀量。

4 结论

(1)研究区土壤风蚀量整体呈现东北部林区向西南部草原区递增的空间分布特征。2000—2018年土壤风蚀量波动下降,土壤风蚀总量共减少1.23×107t,年均下降率为0.70%;从风蚀强度变化来看,2000—2018年研究区有90.86%的区域土壤风蚀强度未发生变化;有9.04%的区域土壤风蚀强度减轻;有0.10%的区域土壤风蚀强度恶化。

(2)2000—2018年间,研究区土壤风蚀量变化的主要驱动因子解释力表现为土地利用变化>降水变化>风速变化>植被覆盖度变化>气温变化,整体来看,各因子变化量解释力q值均很小,表明单一因子变化量对土壤风蚀量变化的驱动作用有限。对各因子变化量的交互作用分析发现,降水量变化与其他因子变化均呈双因子非线性增强作用,其中,降水变化协同土地利用变化与降水变化协同植被覆盖度变化可显著增强对土壤风蚀量变化的驱动作用。

(3)研究区草地覆盖度增加、沙化土地封育、生态退耕、天然林保护等转换方式对土壤风蚀量降低的贡献率为90.56%,说明该区域实施的沙化土地封育、天然林保护、退耕还草、退耕还林等一系列生态保护措施,对该地区的生态环境改善具有重要的推动作用。今后还需进一步采取禁牧、休牧、轮牧等措施,恢复和提升草地覆盖度,从而更有效地减少区域风蚀量。