2000—2015年锡林郭勒盟防风固沙服务功能变化驱动因素分析

张 玥,许端阳,王子玉,张晓宇

1 中国科学院地理科学与资源研究所,北京 100101 2 中国科学院大学,北京 100049 3 北京林业大学,北京 100083

中国是世界上受风蚀危害最严重的国家之一,国家林业局第五次沙化土地监测结果表明我国沙化土地面积达172.12×104km2,占国土面积的17.93%[1]。在干旱和半干旱沙区,防风固沙服务功能是沙区植被生态系统的存在对风蚀作用的削弱功能,保障沙区人民可持续地进行生产生活[2- 5]。近年来我国相继实施了“三北”防护林、退耕还林还草、天然林保护、京津风沙源治理、草原禁牧、退牧还草等生态治理工程,使得我国北方沙区防风固沙服务功能得到有效提升。与此同时,我国北方地区在全球变暖背景下出现的气候暖干化趋势以及城镇化发展中人为干扰强度的增大导致的局部地区土壤风蚀危害不容忽视。因此,开展区域尺度的防风固沙服务功能演变以及驱动机制研究对于区域生态科学治理与逐渐改善具有重要意义。

近年来,许多学者应用多源数据、采用不同方法对防风固沙服务功能演变及驱动因素进行了研究[6- 11]。在影响防风固沙服务功能的气候因素方面,一些学者利用相关系数、趋势分析方法探讨风速、温度、降雨量、植被盖度等因子与土壤风蚀的关系,如Garbrecht等在美国大平原南部农田区域的研究表明温度升高和降雨减少会导致土壤干燥、降低农作物产量,从而加剧土壤风蚀,也有一些研究分析未来气候情境下防风固沙服务功能演变过程[12- 15]；在人为因素方面,多元回归分析、随机过程模型、叠加分析等方法被用来评估放牧、开垦土地、人工造林、土地利用变化等因子对防风固沙服务功能的影响,如Li等基于多元回归分析认为造林面积增加是内蒙古地区防风固沙服务功能提升的主要因素,Sharratt等发现哥伦比亚高原农作物产量的提高能够减轻土壤风蚀危害[16- 19]。上述已有研究对认识和理解防风固沙服务功能空间分布格局的驱动机制具有重要意义,但这些方法的不足是无法处理类别变量,受到较多假设条件制约,在充分挖掘防风固沙服务功能与其驱动因素之间的关联信息,特别是不同驱动因素交互作用方面存在局限性[20- 21]。

地理探测器是基于统计学方法开发的揭示地理现象空间分布格局与其影响因素之间关系的模型,数值型数据与定性数据均可带入模型运算以及探测因子交互作用的特点使得地理探测器在气象、环境污染、生态、人类健康、区域规划等领域被广泛运用[22- 28]。防风固沙服务功能是沙区多种自然和人为因素相互作用的体现,与多元回归分析、主成分分析等传统方法相比,应用地理探测器在分析其驱动因素、气候和人为因子耦合研究方面具有独特优势。因此,本研究选择内蒙古锡林郭勒盟为研究区,在评估2000—2015年该区域防风固沙服务功能时空变化基础上,应用地理探测器模型分析自然和人为因子对锡林郭勒盟防风固沙服务功能空间分布格局的影响,并结合趋势分析法检测结果探讨不同影响因子在防风固沙服务功能演变过程中的交互作用,旨在为科学认识区域防风固沙生态系统服务功能形成机制、实现对区域生态治理的精准管控以及沙区生态补偿机制建立提供有力支撑。

1 研究区与数据

1.1 研究区概况

锡林郭勒盟位于内蒙古自治区中东部北侧,地理位置处于北纬42°02′—47°77′和东经111°21′—120°12′之间,面积为20.3万km2,总人口105.16万,是隔断蒙古沙源与京津冀区域的重要植被沙障(图1)。海拔在800—1800 m之间,地势自西南向东北倾斜,其中浑善达克沙地位于锡林郭勒盟中部,由西北向东南呈条带状分布,属于半固定沙地。冬寒夏燥,春秋多风,属于中温带干旱、半干旱大陆性气候,年均温为0—4 ℃；雨季短促,年降水量250—350 mm,由东南向西北递减,且降雨多集中在夏秋季节；太阳辐射强,蒸发量大(1500—2700 mm)且由东向西递增。草地是锡林郭勒盟的主要土地覆被类型,该区域草原类型空间分布与水热条件密切相关,由西南到东北依次为荒漠草原、典型草原和草甸草原,土壤类型以黑钙土、栗钙土、风化土和棕钙土为主。锡林郭勒盟是京津风沙源治理工程的主体区域,在京津风沙源治理一期工程结束时该区域完成沙源治理任务21500 km2,其中造林面积6200 km2,流动半流动沙地面积减少4500 km2。

图1 研究区地理位置Fig.1 Location of study area

1.2 数据来源与处理

本研究使用的温度、降雨、风速、太阳辐射等气象数据源自中国气象数据网；雪盖数据、土壤质地数据来源于寒区旱区科学数据中心,其中雪盖数据从中国雪深长时间序列数据集(1979—2016年)中提取,空间分辨率为0.25°,土壤质地数据来源于中国土壤特征数据集,空间分辨率为1 km；归一化植被指数(Normalized Vegetation Index,NDVI)数据来源于美国国家航空航天局(NASA)戈达德航天中心的GIMMS NDVI3g数据集,时间分辨率为半月,并进行几何校正、图形增强等处理；1∶100万土壤类型、1∶100万植被类型、1∶100万地貌类型数据、土地利用数据以及GDP和人口密度数据在中国科学院资源环境科学数据中心网站下载,空间分辨率均为1 km；DEM以及坡度数据来源于地理空间数据云。本研究所使用的牲畜数量数据来源于相关年份锡林郭勒盟统计年鉴,人工造林面积数据查阅相关年份中国林业统计年鉴。为体现人为因子作用空间异质性,本研究综合考虑数据可获取性、锡林郭勒盟实际情况并参考相关文献资料,将旗县尺度的牲畜数量和人工造林面积结合土地利用数据赋予不同权重离散到栅格尺度上,其中牲畜数量数据按照高覆盖度草地、中覆盖度草地、低覆盖度草地权重分别为0.5、0.3、0.2进行离散；人工造林面积数据按照裸土地、坡度大于6°耕地、沼泽地、盐碱地、沙地权重分别为0.25、0.25、0.2、0.15、0.15进行离散。为便于模型运算,本研究所用栅格数据统一重采样为1 km×1 km,采用Krasovsky_1940_Albers投影。

2 研究方法

2.1 防风固沙服务功能评估

(1)

PD=109.8×(CF×SCF×SEF×SRF)

(2)

PL=150.71×(WF×SCF×SEF×SRF)-0.3711

(3)

(4)

RD=109.8×(CF×SCF×SEF×SRF×VF)

(5)

RL=150.71×(CF×SCF×SEF×SRF×VF)-0.3711

(6)

SF=PWE-RWE

(7)

式中,PWE代表潜在土壤风蚀量(kg/m2)；PD指潜在风力的最大输沙能力(kg/m)；PL为潜在关键地块长度(m)；RWE代表实际土壤风蚀量(kg/m2)；RD指实际风力的最大输沙能力(kg/m)；RL为实际关键地块长度(m)；SF表示单位面积防风固沙量(kg/m2)；l代表下风向距离(m)；CF表示气候因子(kg/m)；SCF表示土壤结皮因子；SEF表示土壤可蚀性因子；SRF为土壤粗糙度因子；VF为植被因子。

2.2 防风固沙服务功能变化趋势分析

本研究采用线性趋势法分析锡林郭勒盟防风固沙量2000—2015年变化趋势。线性趋势法在栅格尺度上对防风固沙量变化趋势进行模拟,能够反映研究时段内防风固沙量变化幅度的空间分布[32]。此外,为识别防风固沙服务功能发生显著变化的区域,利用t检验法(P<0.05)进行显著性检验。趋势分析法的计算公式如下：

(8)

式中,m代表拟计算时间长度；SFj为第j年的年均防风固沙量。当Trend为正时,表明该区域防风固沙服务功能具有增强趋势；当Trend为负时表明具有减弱趋势。

2.3 防风固沙服务功能驱动因素交互作用分析

本研究应用地理探测器模型分析锡林郭勒盟防风固沙服务功能变化过程中不同驱动因素的交互作用,其基于空间方差分析,兼容类型数据(如土壤类型)和连续数据(如降雨量),包括分异及因子探测、交互作用探测、风险区探测和生态探测四部分(http://www.sssampling.org/GeogDetector/)[20,28]。

(1)因子探测器中潜在影响因子对防风固沙服务功能的解释力计算公式如下：

(9)

(2)交互作用探测用于判断不同影响因子共同作用于防风固沙服务功能时的交互效果,识别影响因子共同作用对防风固沙服务功能的解释力是否有增强或减弱作用,或者是彼此独立的。交互作用类型如图2所示。

(3)风险区探测用于判断两个影响因子分区间的属性均值是否有显著的差别,可用于揭示防风固沙服务功能较强区域的因子分区特征,用t统计量来检验：

(10)

(4)生态探测用于评估两影响因子Xa比Xb对防风固沙量空间分布的影响是否具有显著性差异,以γ统计量来判断：

(11)

(12)

式中,SXa与SXb分别为影响因子Xa与Xb的样本数目；VXa和VXb分别表示因子Xa与Xb分层之间的层内方差和；Na和Nb代表因子Xa与Xb的分层数目。

图2 地理探测器原理示意与交互作用类型[28]Fig.2 Geographical detector principle diagram and interaction typesM：某个研究区域；SF1、SF2、SF3：用于探测的因变量,例如在本文中为防风沙量；C1、C2、C3、C4：影响因子C的子分区；D1、D2、D3、D4：影响因子D的子分区；PD,SF：影响因子对防风固沙服务功能的解释力

防风固沙服务功能的产生受到土壤、地形等自然本底因素的制约,同时气温、风速等气候因素以及禁牧造林等人类活动也会对其产生影响,因此本研究选取包括自然和人为两大类的12个影响因子分析锡林郭勒盟防风固沙服务功能变化驱动机制(表1)。在ArcGIS 10.2中利用自然断点法将年均降水量、年均温度、年均风速、高程、坡度、牲畜数量、GDP、人口密度、人工造林面积划分6类；同时将土壤类型、植被类型、地貌类型分别划分为7类、7类和5类。为深入探究锡林郭勒盟防风固沙服务功能时空变化驱动机制及不同影响因素作用差异,本研究从防风固沙服务功能空间格局形成及时间变化两个维度进行地理探测；其中,变化区域探测针对防风固沙服务功能增加区域和降低区域分别开展。

表1 探测因子指标

3 结果分析

3.1 防风固沙服务功能时空变化分析

锡林郭勒盟2000—2015年单位面积防风固沙量在空间分布上差异明显,总体上表现为东南部温带落叶阔叶林及灌丛区域防风固沙服务功能高、东北部草甸草原地区防风固沙服务功能较高、西部荒漠区域防风固沙服务功能低的分布格局(图3)。从行政区划来看,正蓝旗平均防风固沙量最高,达14.68 kg/m2,正镶白旗、多伦县和东乌珠穆沁旗平均防风固沙量较高,均在6 kg/m2以上,镶黄旗、二连浩特市和苏尼特右旗最低,不足4 kg/m2。

通过对锡林郭勒盟2000—2015年单位面积防风固沙量变化趋势及显著性分析,可以发现2000—2015年防风固沙量总体上呈增加趋势,增加区域面积为115104 km2,其中显著增加的面积为27839 km2,占总增加面积的24.19%,主要分布在苏尼特左旗中东部、多伦县和东乌珠穆沁旗部分区域。与此同时,防风固沙量呈减少趋势的面积为84756 km2,其中显著减少的面积为24985 km2,占总减少面积的29.48%,集中分布在太仆寺旗、正镶白旗、西乌珠穆沁旗和镶黄旗南部区域(图4)。

图3 2000和2015年锡林郭勒盟防风固沙服务功能空间格局Fig.3 Spatial pattern of windbreak and sand-fixing function in Xilingol League in 2000 and 2015

图4 2000—2015年锡林郭勒盟防风固沙服务功能变化趋势Fig.4 Change trends of windbreak and sand-fixing function in study area between 2000 and 2015

3.2 防风固沙服务功能空间格局形成的驱动因素交互作用分析

本研究运用因子探测器探索影响因子在锡林郭勒盟防风固沙服务功能空间分布格局中的重要程度,因子探测器能够量化潜在影响因子对防风固沙量的解释力。各因子对防风固沙量的解释力(PD,SF)排序为：土壤类型(64.07%)>植被类型(19.97%)>人工造林面积(18.49%)>牲畜数量(14.44%)>年均风速(12.65%)>年均降水量(10.5%)>地貌类型(10.45%)>人口密度(9.45%)>GDP(8.96%)>年均温度(7.73%)>高程(4.93%)>坡度(0.22%)。

根据各探测因子的PD,SF值,可以看出土壤类型是影响防风固沙服务功能最重要的因子,植被类型、人工造林面积、牲畜数量、年均风速对防风固沙量空间分布也有较强的影响,解释力均在12%以上,年均温度、GDP等因子对防风固沙量空间分布影响较小,高程与坡度两个因子几乎没有影响。在生态探测结果中,大多数影响因子之间均具有统计显著性差异(表2)。因子探测与生态探测结果相结合表明土壤类型在锡林郭勒盟防风固沙服务功能空间分布格局中发挥最重要的作用,其PD,SF值最高且与其它因子具有显著性差异。

表2 探测因子之间统计显著性

风险探测器用于探测各因子分区中防风固沙服务功能最高的区域,从而揭示防风固沙功能空间格局形成的内在机制。风险探测结果如表3所示,从中可以看出随年均降雨量的增加,单位面积防风固沙量逐渐增大,在锡林郭勒盟年均降雨量为357—490 mm时,防风固沙量均值达到最大值,为8.823 kg/m2,表明年均降雨量的增加通过促进植被生长、增大土壤湿度等方式使得防风固沙服务功能得到有效提升；防风固沙量均值随着年均温度和年均风速分别呈现先逐渐增高、后快速减少的趋势和波动变化,在2.45—3.22 ℃和2.78—3.08 m/s范围时达到最大值,分别为8.03 kg/m2和8.397 kg/m2。随土壤、植被、地貌等类型的不同,单位面积防风固沙量均值波动变化,荒漠风沙土、草原风沙土、石质土、酸性粗骨土、钙质粗骨土等土壤类型的防风固沙量最高，均值为18.08 kg/m2，温带落叶阔叶林、温带落叶小叶疏林等植被类型以及低海拔丘陵、中海拔丘陵等地貌类型的防风固沙量最高,均值分别为17.289 kg/m2和8.585 kg/m2。随高程、坡度的增加,防风固沙量呈现先增加后减少的趋势,分别在1249—1371 m和1.22—2.75 °范围时达到最大值9.328 kg/m2和6.642 kg/m2。在影响锡林郭勒盟防风固沙服务功能的人为因子中,防风固沙量基本随着各人为因子增加呈现出先波动增大后减少的趋势,在牲畜数量为122—270头、人工造林面积为8.51—11.35 hm2时达到最大值,分别为12.035 kg/m2和10.548 kg/m2；防风固沙量随人口密度和GDP的增加表现为先增加后减少再增加到最大值的波动变化,防风固沙量分别在人口密度和GDP为4—16人/km2和42—101万元/km2范围时达到最大值8.972 kg/m2和8.807 kg/m2。

表3 防风固沙量最大均值因子分区特征(置信水平95%)

3.3 防风固沙服务功能变化的驱动因素交互作用分析

3.3.1防风固沙服务功能增强区域

在锡林郭勒盟2000—2015年间防风固沙服务功能增强区域,各因子对防风固沙量的解释力(PD,SF)排序为：土壤类型(56.2%)>年均风速(15.55%)>人工造林面积(12.03%)>植被类型(11.63%)>牲畜数量(10.72%)>地貌类型(10.63%)>年均降水量(9.42%)>年均温度(7.16%)>人口密度(5.34%)>GDP(5.06%)>高程(1.24%)>坡度(0.43%)。根据PD,SF值可以看出土壤类型是影响该区域防风固沙服务功能空间布局的主要因子,年均风速、人工造林面积、植被类型、牲畜数量、地貌类型也具有较强的影响,解释力均在10%以上,且自然因子的影响总体高于人为因子。

本研究对该区域2000年、2005年、2010年、2015年的防风固沙量与各影响因子图层相应属性数据分别进行地理探测(图5),可以发现土壤类型、植被类型的解释力相对稳定,呈现微弱增加趋势；年均降水量、年均温度、年均风速等其它自然因子的解释力基本表现为减少趋势；而人工造林面积、城镇化率、农作物耕种面积与牲畜数量的解释力呈现出先增加后减少的趋势,除牲畜数量外,其他人为因子的解释力均在2010年达到最大值,在2015年则有较大幅度减少。

图5 2000—2015年防风固沙服务功能变化区域PD,SF值动态变化Fig.5 Dynamic changes ofPD,SF value in the area of windbreak and sand-fixing function changes between 2000 and 2015X1、X2、X3、X4、X5、X6、X7、X8、X9、X10、X11、X12表示探测因子,具体见表1

防风固沙服务功能增强区域交互探测器结果显示各因子对锡林郭勒盟防风固沙功能影响存在交互作用,且任意两种探测因子交互作用的解释力高于单个因子(表4)。从表4可以看出除地貌类型外,土壤类型与其它因子的交互作用关系均为非线性增强,且其PD,SF值最高,表明土壤类型在该区域防风固沙服务功能空间格局中发挥着最重要的作用；除高程和坡度因子外,地貌类型与其它因子的交互作用均为双因子增强,主要是由于该区域地貌类型大多是平原和低海拔台地,空间异质性不太明显；值得注意的是虽然高程的PD,SF较小,但与其它因子的交互作用类型均表现非线性增强关系,这表明高程通过影响植被生长和放牧活动等因素在该区域防风固沙服务功能分布格局形成中也具有重要作用；人为因子间的交互作用主要为双因子增强关系,结合人为因子的PD,SF值可以发现人为因子在该区域对防风固沙服务功能的影响略低于自然因子。

表4 防风固沙服务功能增强区域因子交互作用类型

3.3.2防风固沙服务功能减弱区域

在锡林郭勒盟2000—2015年间防风固沙服务功能减弱区域,各因子对防风固沙量的解释力(PD,SF)排序为：土壤类型(76.1%)>植被类型(24.55%)>人工造林面积(22.51%)>牲畜数量(16.35%)>高程(12.72%)>年均温度(12.15%)>人口密度(11.98%)>GDP(11.29%)>年均风速(10.98%)>地貌类型(10.82%)>年均降水量(10.15%)>坡度(0.96%)。根据PD,SF值可以看出各影响因子的解释力总体较高,这表明地理探测器在该区域具有更强的适用性,土壤类型的解释力高达76.1%,在影响该区域防风固沙服务功能空间分布方面具有绝对优势,植被类型和人工造林面积的PD,SF值均高于22%,表明这两个因子也具有相当重要的影响,且人为因子的影响总体高于自然因子。

从各影响因子的解释力时间变化图来看(图5),土壤类型、植被类型、地貌类型等自然本底背景因子的解释力相对稳定；年均降水量、风速、温度等气候因子的解释力表现为增加趋势；在人为因子中,除GDP的解释力呈现先减少后增加外,其它3个因子的解释力均为大幅度增加趋势,尤其是牲畜数量和人工造林面积,2015年PD,SF值分别为38.51%和23.73%。

该区域交互探测器结果亦显示出各因子之间交互作用显著(表5),任意两因子对防风固沙服务功能的影响均高于单个因子；牲畜数量与其它因子间的交互作用均为非线性增强关系,表明放牧活动对该区域防风固沙量具有重要影响；人为因子间的交互作用主要表现为非线性增强关系,结合人为因子的PD,SF值进一步证明人为因子对该区域防风固沙服务功能的影响高于自然因子。

表5 防风固沙服务功能减弱区域因子交互作用类型

4 讨论与结论

在影响锡林郭勒盟防风固沙服务功能演变的探测因子中,与年均温度、风速等气候因子相比,年均降水量对防风固沙量的影响相对较高,证实了降水在较干旱地区促进植被生长从而减少风蚀发生的重要作用,这与前人的研究成果一致[33-36]；温度升高有利于植物发育降低土壤侵蚀,然而随温度升高带来的水分蒸发、降低土壤湿度的影响有可能会抵消对防风固沙功能提升的积极作用；锡林郭勒盟风速最低的多伦县由于耕地面积较多,翻种耕地造成土壤疏松带来的风蚀作用较大。因此,适宜的气候因子范围对于防风固沙服务功能提升具有重要影响,同时不同影响因子之间的精准契合有助于单位面积防风固沙量的增加。

在不同的探测尺度上,土壤类型、植被类型对防风固沙服务功能空间分布的解释力均较高,但在防风固沙服务功能增强区域,自然因子的解释力总体相对较高,在防风固沙服务功能减弱区域,人为因子的解释力相对较高,说明人为因子的空间分布与防风固沙量减少区域分布更加吻合。在功能增强区域,人为因子的PD,SF值表现为先增加后减少的趋势,可能的解释是近年来禁牧造林等政策,比如该区域2002年以来实施的京津源风沙治理工程有利于防风固沙服务功能得到较大幅度提升,但人类活动的作用发挥到一定程度后,仍然是自然因子发挥着主导作用。在功能减弱区域,人为因子的解释力基本表现为大幅度增加趋势,对比该区域2000年和2015年的土地覆盖数据,发现变化较大的是低覆盖度草原面积增加了761 km2,沙地和建设用地面积增加360 km2,同时高覆盖度和中覆盖度草原面积减少了1189 km2,这表明禁牧造林等保护政策的实施对防风固沙服务功能的提升作用可能比较微弱,但若不加以保护,沙区土地退化速度将大幅度增加。这从防风固沙服务功能增强和减弱区域人为因子之间的交互作用关系也可以看出,在增强区域人为因子主要表现为双因子增强,而在减弱区域主要表现为非线性增强,交互作用更加显著。

本研究以锡林郭勒盟为研究区,采用修正风蚀模型评估其2000—2015年防风固沙服务功能时空变化,应用地理探测器分析自然和人为因子对该区域防风固沙服务功能空间格局形成及变化过程中的贡献及交互作用。结果表明防风固沙服务功能空间分布差异显著,防风固沙量总体上呈增加趋势；土壤类型在该区域防风固沙服务功能空间分布格局中发挥最重要的作用,植被类型、人工造林面积、牲畜数量、年均风速等因子对防风固沙量空间分布影响较强；防风固沙服务功能变化区域的2000—2015年动态地理探测结果显示土壤类型、植被类型均具有较高的解释力且相对稳定,功能增强区域气候因子的PD,SF值呈现出减少趋势,人为因子表现为先增加后减少趋势,功能减弱区域变化趋势大体与之相反,探测因子之间交互作用显著。本研究认为土壤类型、植被类型、牲畜数量等是影响研究区防风固沙服务功能空间分布格局形成的主要因素,这与前人研究成果不太一致,例如大多研究认为风速、植被盖度、温度等是影响防风固沙服务功能的主要因素,主要是由于本研究中应有的地理探测器模型兼容类型数据和连续数据,而相关系数、多元回归、主成分分析等方法不适用于分析类型量,如土壤类型、植被类型、地貌类型等,这些自然本底背景在影响防风固沙服务功能区域分异中具有明显作用。需要指出的是,本文考虑到数据尺度匹配问题将旗县尺度的牲畜数量和人工造林面积结合土地利用数据赋予不同权重离散到栅格尺度上,这可能会对结果产生一定影响,但该处理方式能表征因子的相对大小,因此应用地理探测器分析仍具有合理性和科学性。本文对防风固沙服务功能变化的驱动因素应用地理探测器模型进行探索性研究,土壤类型、地貌类型等是相对稳定的自然因素,而温度、降雨、风速等变率较大的自然因子与人为因子的耦合研究将是下一步有待解决的关键科学问题,例如在未来气候变化情境下牲畜数量、人工造林面积的最适宜范围确定将对区域生态环境科学治理恢复与经济高质发展之间的协调提供有效参考。