风蚀是土壤颗粒在风力作用下侵蚀、搬运和堆积的过程, 也是土地退化和荒漠化的主要原因^[1]。我国北方干旱半干旱地区是受风蚀现象影响较为严重的地区, 是沙漠化防治的重点区域^[2]。防风固沙服务作为风蚀地区生态系统提供的一项重要防护型服务, 对风蚀发生地与周边地区风沙灾害治理、生态环境恢复具有重要意义。浑善达克地区位于我国内蒙古高原干旱半干旱地区, 是我国北方农牧交错带, 生态环境脆弱, 腹地有浑善达克沙地, 是京津地区风沙主要来源地^[2—3]。2000年以来, 随着京津风沙源治理工程、三北防护林工程、退耕还林还草、退牧禁牧政策的实施, 浑善达克地区生态系统状况呈现改善趋势, 防风固沙服务供给能力也有一定的提升, 风蚀现象有所好转^[4—6]。然而, 气候变化、人类活动与生态系统之间频繁的相互影响, 使得环境敏感区的生态系统结构并不稳定, 其提供的生态系统服务也相对有限^[7—8]。因此, 开展长时间序列的浑善达克地区防风固沙服务时空变化及驱动因素分析对于区域生态系统恢复以及生态建设、土地管理具有重要意义。

防风固沙服务一般通过潜在风蚀量和实际风蚀量之差计算获得。潜在风蚀量和实际风蚀量可通过实地观测和模型模拟确定^[9—12]。Hu等^[9]利用¹³⁷Cs和²¹⁰Pb示踪法模拟浑善达克南部地区的风蚀状况。然而, 这种方法需要花费巨大的人力、物力, 且仅能代表风蚀测定样点小范围的结果, 难以用于区域风蚀量确定。基于模型的风蚀量估算能显著提高计算效率, 且可用于大面积区域的风蚀量计算^[13—14]。修正风蚀方程模型(Revised Wind Erosion Equation, RWEQ)因子参数相对全面、模型构成较为简单、数据更易获取^[13—15], 是应用最为广泛的防风固沙服务计算方法。目前, 已经有较多基于RWEQ模型的防风固沙量的模拟研究。这些区域防风固沙服务模拟以多年研究为主^[16—18], 以避免单个年份的气象数据带来的随机误差。Li等^[19]基于RWEQ模型分析了1990—2015年我国西北干旱地区土地利用变化对防风固沙服务的影响。徐洁等^[20]和张彪等^[4]利用RWEQ分别分析了防风固沙型重点生态功能区的防风固沙服务能力和京津风沙源工程对防风固沙服务产生的效益。王俊枝等^[21]通过RWEQ模型模拟了内蒙古自治区境内的浑善达克生态功能区1990—2015年防风固沙服务的时空分布特征。长时间序列的防风固沙服务模拟有助于全面分析一个区域防风固沙服务的真实状况及其时空格局变化, 但现有的研究主要集中在1990年以后的长时间序列, 更长时间序列的研究还不多见。

由于风蚀现象受到自然和人为等多种因素综合作用影响, 开展防风固沙影响因素分析可以为提高防风固沙服务提供科学依据。目前, 大部分研究采用线性相关性模型^[22—24]、结构方程模型^[25]、约束效应^[26—27]等方法分析植被覆盖度、风速、降水、畜牧业发展、土地利用等因素与防风固沙服务之间的关系^[28—29]。关于浑善达克地区防风固沙服务驱动因素分析的分析相对较少。申陆等^[24]通过主成分分析法探究浑善达克生态功能区单位面积防风固沙量变化的驱动因素。高晓霞等^[23]则通过Logistic回归模型探究浑善达克沙地动态变化的驱动因素。虽然上述研究为探讨防风固沙服务时空格局变化机制提供很好的参考, 但多数研究为时间序列的研究, 而且忽略了多因素之间的共同作用对防风固沙服务的影响。从空间分布格局变化的角度探究防风固沙服务变化驱动机制的研究也相对较少, 类型变量对防风固沙服务的影响也难以衡量。

地理探测器是一种基于空间分异性, 揭示某一变量驱动机制的空间统计学方法^[30], 它能够在探测类型与数值两种不同属性的变量对因变量空间分布影响的同时, 分析自变量两两之间的交互作用及具体作用类型。目前, 在资源环境、社会经济等多领域的空间格局变化过程分析中得到了广泛的应用^[31—34], 但用于分析防风固沙驱动因素的研究还较少。应用地理探测器分析防风固沙服务空间格局变化的驱动因素, 能够兼顾类型变量与数值变量对服务变化的影响, 分析各因素之间的相互作用, 弥补防风固沙服务空间分布驱动因素分析方面的不足。

因此, 本研究以浑善达克重点生态功能区为例, 通过RWEQ模型模拟1980—2018年防风固沙服务的时空变化趋势, 利用地理探测器开展防风固沙服务空间变化的驱动因素分析, 为浑善达克地区生态建设与恢复以及土地管理提供科学依据。

1 研究区与数据 1.1 研究区概况

浑善达克地区(40.7°N—45.5°N, 111.1°E—118.5°E)地处我国内蒙古高原中段, 总面积16.39万km², 跨越内蒙古自治区与河北省。该区空间范围依据国家重点功能区中的浑善达克沙漠化防治生态功能区范围划定, 包括内蒙古8旗1县(克什克腾旗、阿巴嘎旗、苏尼特左旗、苏尼特右旗、太仆寺旗、镶黄旗、正镶白旗、正蓝旗、多伦县)与河北省6县(丰宁满族自治县、围场满族蒙古族自治县、张北县、康保县、沽源县、尚义县)(图 1)。浑善达克地区地势西南高、东北低, 平均海拔1300 m。该地区属于中纬度干旱、半干旱大陆性季风气候, 年平均气温0—3℃, 大部分地区年降水量在250—400 mm, 年平均相对湿度为60%。春秋季风沙多, 年均风速度为4—5 m/s, 最大风速普遍在24—28 m/s, 年大风日数(>8级)大约60—80 d。草地为主要的土地覆被类型, 自东向西依次为草甸、草原、荒漠草原；土壤类型为栗钙土、棕钙土, 非地带性土壤为风沙土。浑善达克地区整体生态环境脆弱, 是我国北方重要生态屏障之一, 同时也是京津地区主要沙源地, 沙漠化是该地区主要的生态问题。针对土地退化、沙漠化等问题, 2000年以来，开始在浑善达克地区实施京津风沙源治理工程, 2010年浑善达克地区被列为国家沙漠化防治生态功能区。在各类生态工程支持下, 浑善达克地区通过围栏封育、划区轮牧、草地补播等措施实施草地恢复与治理, 沙漠化趋势逐渐减缓^[35]。

1.2 数据来源

本研究使用的数据年份为1980年、1990年、1995年、2000年、2005年、2010年、2015年、2018年。20—20时日累计降水量、日均温、10 m处日均风速等气象日值数据来自中国气象科学数据共享服务网(http://data.cma.cn/)。土壤质地数据来自第二次全国土地调查, 中国科学院南京土壤研究所提供的1:100万土壤数据。雪盖数据来源于寒区旱区科学数据中心(http://bdc.casnw.net)的“中国雪深长时间序列数据集(1979—2018年)”, 空间分辨率为0.25°。2000年之前的归一化植被指数(Normalized Difference Vegetation Index, NDVI)数据来源于GIMMS NDVIg3数据集, 在Google Earth Engine上进行图像计算和提取等处理；2000年之后的NDVI数据(分辨率为1 km)、全国土地利用/覆被数据(空间分辨率为30 m)、GDP空间分布数据(分辨率为1 km)、1:100万地貌类型、1:100万土壤类型、1:100万植被类型数据, 以及30 m高程数据均在中国科学院资源环境科学数据中心(https://www.resdc.cn)下载。路网数据来源于OpenStreetMap(https://www.openstreetmap.org)。本研究所使用的年粮食播种面积、年粮食产量、年肉类产量、年末牲畜数量均来自各旗县及其所属盟、市的相关年份的统计年鉴并求得多年平均值, 人工造林面积数据从《中国林业统计年鉴》中获取。在模型运算过程中, 本研究所用栅格数据均采用Krasovsky_1940_Albers投影, 并统一重采样为30 m分辨率。

2 研究方法 2.1 防风固沙服务物质量模拟

在充分考虑气候、地形、土壤理化性质以及植被覆盖等因素后, 本研究采用修正风蚀方程RWEQ(Revised Wind Equation)模型, 定量模拟浑善达克地区风蚀量的时空动态变化。根据RWEQ模型以及国内实际情况进行参数的调整^{[13, 36—38]}, 计算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中, SL为实际风蚀量(kg m^-2 a^-1)；Q_max为实际风力的最大输沙能力(kg/m)；s表示实际关键地块长度(m)；SR表示潜在风蚀量(kg m^-2 a^-1)；Q_rmax表示潜在风力的最大输沙能力(kg/m)；s_r表示潜在关键地块长度(m)；G为防风固沙的物质量(kg m^-2 a^-1)；z为所计算的下风向距离(m), 本次计算取50 m；WF表示气候因子(kg/m)；EF为土壤可蚀性成分(无量纲)；SCF表示土壤结皮因子(无量纲)；K′表示地表糙度因子(无量纲)；C表示植被因子(无量纲)。

(1) 气候因子WF

气候因子主要考虑了降水、气温、风速、雪盖等自然因素影响下风力对土壤的搬运能力, 其中, 风力因子为气候因子中风速数据的相关计算, 主要计算公式如下：

(8)

(9)

式中, WF为气候因子(kg/m)；Wf为风力因子(m/s³)；g为重力加速度(m/s²)；ρ为空气密度(kg/m³)；SW为土壤湿度因子；SD为雪盖因子(无积雪覆盖天数/研究总天数), 定义雪盖深度＜25.4 mm时为无积雪覆盖；u₁为2 m处起沙风速, 取5 m/s；u₂为由10 m转换为2 m处的月均风速值(m/s)；Nd为各月风速大于5 m/s的天数。

(2) 土壤可蚀性因子EF、土壤结皮因子SCF

土壤可蚀性因子与土壤结皮因子分别表示在一定理化条件下土壤受风蚀影响的程度与土壤结皮抵抗风蚀能力的大小, 二者的表达式分别为：

(10)

(11)

土壤质地数据为国际标准, 采用对数正态分布模型将土壤粒径转换为RWEQ模型需要的美国标准。其中, sa为土壤粗砂含量(%)；si为土壤粉砂含量(%)；cl为土壤粘粒含量(%)；OM为土壤有机质含量(%)；CaCO₃为碳酸钙含量(%), 本研究中取0。

(3) 植被因子C

植被因子表示一定的植被条件对风蚀的抑制程度。

(12)

(13)

式中, SC为植被覆盖度(%), NDVI、NDVI_max、NDVI_min分别为NDVI的实际值、最大值和最小值。

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

地表糙度因子用于模拟地表粗糙程度对风蚀的影响因子。

(14)

式中, α为坡度, 由30 m的DEM数据经过ArcGIS的slope模块计算得到。

2.2 防风固沙服务保有率模拟

RWEQ模型计算的防风固沙量受降水、风场强度等气候变化的影响较大, 不能准确地评估生态系统自身防风固沙的效果。通过计算防风固沙服务保有率能够避免气候因素对服务功能的影响, 进一步分析防风固沙服务能力。防风固沙服务功能保有率计算公式为：

(15)

式中, F表示防风固沙服务功能保有率；G为防风固沙量(kg m^-2 a^-1)；SR为潜在风蚀量(kg m^-2 a^-1)。

2.3 防风固沙服务变化趋势分析

本研究采用最小二乘法线性拟合函数分析浑善达克地区防风固沙服务物质量1980—2018年的变化趋势及其空间分布变化。据1980—2018年重大生态工程与生态建设实施的年份为关键节点, 在栅格尺度上分析1980—2018年、1980—2000年、2000—2010年以及2010—2018年四个时间段浑善达克地区防风固沙服务物质量变化趋势及空间分布特征。

(16)

式中, Trend为防风固沙服务物质量或保有率的变化趋势；n为研究总年数；v_i为第i年防风固沙服务物质量年均值。Trend＞0, 表明防风固沙能力呈增加趋势；Trend＜0, 表明防风固沙能力呈下降趋势。

2.4 防风固沙服务驱动因素分析 2.4.1 地理探测器

地理探测器主要包括4个探测工具：因子探测、风险区探测、生态探测、交互作用探测^[30]。

(1) 因子探测器通过层内方差之和(Within Sum of Squares, SSW)与全区总方差(Total Sum of Squares, SSW)计算各探测因子对防风固沙服务物质量空间分异的解释力, 用q值度量, 计算公式如下：

(17)

式中, h=1, 2, …, L为防风固沙服务物质量或影响因子的分层；N_h和N分别为各层和全区的单元数；σ_h²和σ²分别为各层和全区的防风固沙服务物质量的方差。q值越大, 相应的影响因子对防风固沙服务物质量空间分析的解释力越大。

(2) 风险区探测用于判断两个影响因子的子区域间的防风固沙物质量均值是否有显著差别, 并通过t统计量来检验：

(18)

式中, Y_h为子区域h内防风固沙物质量均值；n_h为子区域h内的样本数量；Var表示方差。

(3) 生态探测通过统计量F评估两个影响因子X1与X2对防风固沙物质量空间分异的影响是否显著：

(19)

式中, N_X1与N_X2分别为影响因子X1与X2的样本数目, SSW_X1和SSW_X2分别表示对应影响因子分层的层内方差之和；L1与L2表示影响因子X1与X2的分层数目。

(4) 交互探测能够判断不同影响因子的共同作用对防风固沙服务物质量的解释力的整体影响(增加、减弱或相互独立)。通过比较影响因子X1与X2对防风固沙物质量各自作用与交互作用的q值(q(X1)、q(X2)、q(X1∩X2))大小, 将交互作用分为5类(表 1)。

2.4.2 基于地理探测器的驱动因素分析

防风固沙服务物质量的变化同时受到气候、土壤、植被类型等自然因素和退牧、禁牧等人类活动的影响。基于RWEQ模型以及防风固沙服务或其他服务驱动因素研究中所涉及的驱动因素的出现频率和研究结果^{[23, 26, 29, 39]}, 本研究选取地貌类型、年降水量、土壤类型等11类自然因素和单位面积国内生产总值(Gross Domestic Product, GDP)、人口密度等9类社会经济因素作为探测因子, 分析浑善达克地区2000、2005、2010、2015、2018年和1980—2018年年均防风固沙服务物质量变化的驱动因素(表 2)。基于浑善达克当地的实际情况以及数据的不同类型, 本研究将浑善达克的地貌类型、土壤类型、植被类型分别分为4类(平原、台地、丘陵、山地)、7类(针叶林、阔叶林、灌丛、荒漠、草原、草丛和草甸、栽培植被)、7类(淋溶土、半淋溶土、钙层土、干旱土、初育土、半水成土、水成土、盐碱土), 数值类变量利用ArcGIS中的自然间断法, 将人口密度、单位面积GDP分为4类, 其余变量分为10类, 分析2000—2018年浑善达克地区单位面积防风固沙量时空变化的驱动因素。

3 结果分析 3.1 防风固沙服务时空动态变化 3.1.1 防风固沙量时空变化

根据浑善达克地区各年单位面积防风固沙量平均值(表 3), 可以发现, 1980—2018年, 浑善达克地区单位面积防风固沙量的变化范围为13.01—22.90 kg m^-2 a^-1, 防风固沙总量为2.18×10¹²—3.81×10¹² kg/a。1980年, 浑善达克全区单位面积防风固沙量和总量均为最高, 随后二者在一定范围内上下波动, 其中, 单位面积防风固沙量在13.01—19.47 kg m^-2 a^-1之间变化, 在2018年小幅度回升。相比于防风固沙总量, 浑善达克地区潜在风蚀量和实际风蚀量有较大幅度的变化。1980年潜在风蚀量最高, 1990年有较大幅度的下降, 在1995年保持稳定, 而后在2000年又有较大幅度下降, 之后一直保持稳定；1980—1995年实际风蚀量基本保持稳定, 在2000年有显著下降, 之后保持稳定。

空间上, 浑善达克地区防风固沙服务总体表现为西北部荒漠、中部沙地以及南部林地普遍高于其他地区的空间分布格局(图 2)。其中, 正蓝旗、正镶白旗、苏尼特左旗的单位面积防风固沙量较高, 数值均在20 kg m^-2 a^-1以上；张北县、沽源县、康保县、尚义县单位防风固沙量较低, 主要在4—10 kg m^-2 a^-1。此外, 单位面积防风固沙量的高值区逐渐从中部偏东地区向西北地区移动。1980—2018年, 浑善达克地区防风固沙服务高值分布面积整体呈现波动减小的趋势。1980年, 高值区域(≥30 kg/m²)面积最大, 为4.63×10⁴ km²；2005年, 面积缩减到最小, 为2.14×10⁴ km²。此外, 2015年单位面积防风固沙量的极高值区域(≥40 kg/m²)面积缩减到最小, 仅为0.63×10⁴ km²；而2018年区域面积又有所回弹, 高值区域面积增长至4.41 km², 其中, 极高值区域面积为1.85×10⁴ km²。

综合1980—2018年浑善达克地区单位面积防风固沙量的结果, 可以发现, 38年来浑善达克地区的风蚀程度有所好转, 潜在风蚀与实际风蚀情况均有不同程度的下降, 防风固沙总量也有小幅下降(表 3)。这可能与内蒙古地区风蚀发生的原动力(即风速)的变化, 以及影响植被固沙作用的气候条件的不稳定有关, 全球气候变暖和内蒙古地区近30年来风速的减小缓解了浑善达克地区风蚀情况^{[6, 19, 26]}。本研究结果还显示, 2000年以来潜在风蚀量、实际风蚀量相比以前有显著下降并保持稳定, 这与2000年以来退耕还林还草和京津风沙源治理工程的实施有较大关系^[35]。

3.1.2 防风固沙保有率时空变化

1980—2018年, 浑善达克地区防风固沙服务保有率自1980年的最低值70.79%波动上升, 至2018年达到最大值94.28%(表 3)。其中, 1980—1995年, 保有率变化幅度较小, 保持在70%左右, 2000年达到89.11%后, 出现一定程度的下降, 2005—2018年, 防风固沙保有率以83.31%为起点, 大致呈现较为稳定的上升态势。

空间上, 浑善达克地区防风固沙服务保有率整体呈现出自东南向西北递减的空间分布特征。东南部的丰宁满族自治县、围场满族蒙古族自治县和克什克腾旗保有率始终高于其他地区, 西北部苏尼特左旗、苏尼特右旗的荒漠及中覆盖度草地地区的保有率则相对较低(图 3)。1980—2018年, 防风固沙服务保有率低值区(保有率＜60%)逐渐向苏尼特右旗北部与中西部沙地边缘地区集中, 且面积显著减小, 其中, 极低值区域(保有率＜40%)面积自1995年开始显著减少。1980年, 低值区面积最大, 以极低值区域为主, 低值区面积为4.83×10⁴ km², 包括极低值面积3.07×10⁴ km²；2000年, 防风固沙服务保有率低值面积显著减小, 而2005年, 面积又扩大至3.06×10⁴ km²。值得注意的是, 极低值的面积没有出现大幅度回升, 且2005年之后的极低值面积变化都不大。2018年, 全区保有率低值区域面积最小, 为0.19×10⁴ km², 其中, 极低值区域面积仅为0.02×10⁴ km²。

整体来看, 浑善达克地区防风固沙服务保有率整体显著提高, 年平均保有率为82.42%。尤其是2000年之后极低值区域面积的显著减少可能与2000年之后实行的京津风沙源治理工程、退耕还林还草等生态恢复工程相关, 说明植被对风沙的抑制作用逐渐显著，生态恢复取得了初步的效果；而2005年, 保有率的短暂性降低, 可能与当年气候相对干旱, 阻碍了植被的生长, 削弱了防风固沙服务功能有关^{[24, 40]}。

3.2 防风固沙服务趋势变化 3.2.1 单位面积防风固沙服务变化趋势

总体而言, 1980—2018年浑善达克地区单位面积防风固沙量以0.43 kg m^-2 a^-1的速度减少。其中, 1980—2000年、2000—2010年, 单位面积防风固沙量分别以1.38 kg m^-2 a^-1和1.04 kg m^-2 a^-1的速度减少, 2010—2018年, 单位面积防风固沙量则以1.77 kg m^-2 a^-1的速度增加(表 4)。

空间格局的变化趋势上来看, 1980—2018年, 浑善达克大部分地区单位面积防风固沙量表现为下降趋势, 增加与显著增加区域集中在苏尼特左旗、苏尼特右旗, 以及正镶白旗、正蓝旗和阿巴嘎旗的部分区域(图 4)。

1980—2000年单位面积防风固沙量变化的整体空间格局与1980—2018年的整体格局大致相同(图 4), 说明2000年以来浑善达克地区单位面积防风固沙量相对稳定。2000—2010年除苏尼特左旗、浑善达克沙地及周边地区为降低或显著降低趋势外, 大部分地区单位面积防风固沙量为增加或显著增加趋势(图 4)。2010—2018年单位面积防风固沙量整体恢复下降的趋势, 增加趋势集中区转移至苏尼特左旗北部、阿巴嘎旗以及浑善达克沙地地区(图 4)。在所有时段内, 苏尼特右旗、苏尼特左旗以及浑善达克沙地地区单位面积防风固沙量均出现明显的增加趋势, 防风固沙能力得到了稳定且连续的提高。

结合表 3可以看出, 1980—2000年潜在风蚀与实际风蚀均呈波动减小趋势, 2000—2010年二者则基本保持稳定。两个时间段单位面积防风固沙量变化趋势呈现明显差异(图 4), 可能由于1980—2000年风速显著下降, 降水年际变化幅度大, 风蚀动力减弱使得大部分地区单位面积防风固沙量为下降趋势；2000年以后, 风速基本保持稳定, 降水增多促进了植被的自然恢复, 使得大部分地区在2000—2010年单位面积防风固沙量为增加趋势(图 5)。

3.2.2 防风固沙保有率变化趋势

1980—2018年浑善达克地区防风固沙保有率以3.46%/a的速度增加(表 4)。其中, 1980—2000年增速较快, 为6.29%/a；2000—2010年转变为减小趋势且速度相对较慢, 为-1.34%/a；2010—2018年保有率恢复增加趋势, 为2.39%/a。

1980—2018年浑善达克地区防风固沙服务保有率总体呈现增加趋势, 其中, 西北荒漠地区及其与中部沙地过渡地区增加趋势显著(图 6)。以2000年、2010年为时间节点, 发现不同时段内, 防风固沙保有率变化趋势明显不同：1980—2000年浑善达克绝大部分地区保有率为增加趋势, 其中, 苏尼特右旗、苏尼特左旗以及正镶白旗为显著增加趋势；2000—2010年以苏尼特右旗、苏尼特左旗、镶黄旗、正镶白旗为主的北部、中部地区呈现降低趋势, 说明该时段内植被对防风固沙服务的影响减弱；2010—2018年浑善达克防风固沙服务保有率又整体恢复增加趋势, 但显著增加区域仅出现在苏尼特右旗中北部的小部分地区(图 6)。

由此可见, 1980—2018年浑善达克地区生态系统防风固沙服务能力有了较为显著的提高, 在2010年之后, 显著增加区域的大幅减少说明浑善达克地区生态系统防风固沙能力已经维持在一个相对稳定的状态。而2000—2010年, 由于风速持续增加导致潜在风蚀量增加, 在单位面积防风固沙量相对稳定的情况下, 该时段内防风固沙服务保有率出现短暂下降。

3.3 防风固沙服务驱动因素分析 3.3.1 单因素作用分析

2000—2018年多年平均的单因素作用结果显示, 对防风固沙量空间变化的影响程度(即解释力q值)较大的探测因子，其影响程度由大到小依次为：土壤类型(75.15%)、年末牲畜数量(25.20%)、年降水量(24.16%)、人工造林面积(21.11%)、年粮食产量(18.96%)、年粮食播种面积(18.04%)、NDVI(16.80%)、植被覆盖度(16.74%)、年肉类产量(16.12%)、年均温(13.98%)、平均风速(11.10%)、高程(10.41%)(图 7)。综合2000、2005、2010、2015、2018年的探测结果, 可以发现, 土壤类型始终是防风固沙服务空间变化的主要驱动因素, 且影响程度稳定在70%左右。同样影响较强且持续时间较久的因子还包括年降水量和平均风速, 解释力分别在24.16%和11.10%左右。

由于自然环境与社会经济条件不同, 各年份防风固沙量的主要影响因子也存在差异：2000年退耕还林还草、京津风沙源治理等生态工程的实施, 年粮食播种面积、年末牲畜数量以及人工造林面积对防风固沙服务空间变化的影响较强。相比于2000年, 2005—2015年年降水量与平均风速对防风固沙服务空间格局变化的影响程度有较为明显的增加, 年降水量的影响相对稳定, 解释力保持在25%左右；平均风速的影响程度则逐年增强。这可能与该时间段内年降水量相对稳定, 风速的持续增强有关。2000—2018年, 由于退牧还草、划区轮牧、草地补播等生态工程的实施, 浑善达克地区年末牲畜数量明显减少并在2010年后保持稳定, 生态系统植被也有所恢复。因此, 年末牲畜数量对防风固沙服务空间变化的影响程度逐年减弱, 而NDVI和植被覆盖度的影响程度则有所增加。

如图 8所示, 单位面积防风固沙量均值随不同探测因子的分类数值的增加而发生不同的变化(探测因子分类间断值见表 5)：单位面积防风固沙量均值随NDVI、植被覆盖度、坡度、起伏度、人口密度、单位面积GDP的增加而波动下降。其中, 植被覆盖度、NDVI在第3分类处对应的单位面积防风固沙量均值最大, 随后为明显的波动下降；而人口密度与单位面积GDP的持续增加, 使单位面积防风固沙量均值呈现不同幅度的持续下降。单位面积防风固沙量均值随平均风速、与城镇的距离、与道路的距离的增加表现为波动上升的趋势。随着人工造林面积、年粮食产量、年粮食播种面积以及年肉类产量增加, 单位面积防风固沙量均值变化较为一致, 均为先波动增加, 在某一分类处达到最大之后, 出现显著下降, 随后保持相对稳定的上下波动。

因此, 在17类数值型驱动因素中, 单位面积防风固沙量均值与平均风速、与城镇的距离、与道路的距离整体表现为正相关关系；与坡度、起伏度、人口密度、单位面积GDP大致表现为负相关关系；与植被覆盖度和NDVI表面为先缓慢增加, 随后显著降低的关系。

在2000—2018年多年平均生态探测的结果中, 年降水量、土壤类型、人口密度与其他所有影响因子之间均存在显著差异, 年均温、NDVI、植被覆盖度、人工造林面积和肉类产量也与其他大部分影响因子之间存在显著性差异, (表 6)。结合因子探测的结果, 可以判断土壤类型、年降水量、人工造林面积是影响浑善达克地区单位面积防风固沙量空间变化的两个最主要因素。

3.3.2 多因素交互作用分析

根据多年平均交互探测器结果(表 7), 发现各因子之间对浑善达克地区防风固沙服务的影响均存在增强型交互作用。除地貌类型、年均温、植被类型、与城镇的距离与其他大部分因子之间为非线性增强关系外, 大部分因子与其他因子之间普遍存在双因子增强关系。土壤类型、年降水量、人工造林面积、牲畜数量与其他因子之间的双因子增强作用对防风固沙服务空间变化的影响程度较大。尤其是土壤类型与人工造林面积、粮食产量、粮食播种面积之间的双因子增强作用最强, 对防风固沙量变化的影响程度(q值)均达到最大值82%。这说明土壤类型也加强了其他因素对防风固沙服务空间分布的影响程度。

而年均温与除地貌与土壤类型外的其他因子之间的非线性增强关系, 在所有因子之间的非线性增强作用中对防风固沙服务空间变化的影响程度最高, 在30%左右。可以看出, 尽管年均温对防风固沙服务空间变化的单因素作用的影响较弱, 但通过改变区域自然环境状态或人类活动, 使得两因子间的交互作用对防风固沙的影响水平高于二者单独作用之和, 从而间接影响了防风固沙服务的空间变化。因此, 年均温对防风固沙量的变化也存在较为重要的间接影响。

4 结论

本研究采用修正风蚀模型模拟了浑善达克地区1980—2018年防风固沙服务及保有率的时空变化, 以重大生态建设实施等年份为时间节点, 分析不同时段内二者的变化趋势, 并通过地理探测器探究年降水量、人口密度等20项自然与人为因子及因子之间的交互作用对该地区防风固沙服务空间格局变化的影响。结果表明：

(1) 1980—2018年, 浑善达克地区风速的下降, 划区轮牧、草地补播等生态恢复政策的实施, 加之适宜植被生长的气候变化, 使得防风固沙服务能力整体得到一定的提高, 单位面积防风固沙量波动减少, 防风固沙服务空间分布差距逐渐缩小；

(2) 土壤类型、年末牲畜数量、年降水量、人工造林面积是造成浑善达克生态功能区防风固沙服务空间差异及分布变化的主要影响因素。分别在初育土、52.84—86.86万头、404.21 mm、5167 hm²时, 各自对应的单位面积防风固沙量均值最大。

(3) 各驱动因素间的交互作用都会放大单因子对浑善达克地区防风固沙服务空间分布的影响。其中, 土壤类型与其他因子的交互作用对防风固沙服务空间变化的影响, 均大于单因子分别对防风固沙服务空间变化的影响, 整体对防风固沙服务空间变化影响程度最大。

(4) 年均温虽然对单位面积防风固沙量空间分布的影响较小, 但当它与除地貌类型与土壤类型的其他因子共同作用时, 对防风固沙服务空间变化影响均高于两个驱动因素各自对防风固沙服务的影响之和, 对单位面积防风固沙量空间分布具有较强的间接影响作用。

5 讨论

在长时间序列下, 浑善达克防风固沙重点生态功能区的防风固沙能力得到了一定程度的提高, 风蚀程度有所缓解, 单位面积防风固沙量高值区从中部、东部沙地向西北荒漠地区转移。特别是在2000年之后, 防风固沙能力与单位面积防风固沙量均稳定在相对较低的区间范围内。结合地理探测器的分析结果, 可以得出, 1980—2018年浑善达克生态功能区的上述变化, 主要与该段时间内土壤类型, 年均风速、年降水量、人工造林面积以及牲畜数量的时空格局变化有关。1980—2000年, 年均风速、年降水量等因素年际变化不稳定, 京津风沙源治理工程尚未开始, 是导致浑善达克地区防风固沙服务能力整体较弱, 风蚀程度较高的主要原因。2000年之后, 围栏丰育、划区轮牧、草地补播等生态保护政策的实施, 使得人工造林面积、牲畜数量在时空分布上较为稳定, 加之年均风速、年降水量等气候条件促进植被恢复, 使得浑善达克生态功能区防风固沙服务能力整体提高, 区域间防风固沙服务能力的差异也逐渐缩小。此外, 在各驱动因素的交互作用中, 土壤类型的时空分布决定了区域内地理条件的差异, 从而放大其他因子对防风固沙服务空间分布的影响, 使二者对防风固沙服务分布的影响大于单因子影响但小于二者各自的影响之和。年均温能够影响除地貌类型和土壤类型外其他因子的时空分布, 当年均温与其他因子共同作用于防风固沙服务空间分布时, 其影响程度高于两个单因子独立的影响程度之和, 具有较强的间接作用。

但是, 两种以上因子之间的交互作用对防风固沙服务能力及时空变化的影响如何, 周边地区对浑善达克生态功能区防风固沙服务是否产生影响, 还需要进一步地深入研究, 为浑善达克地区生态建设与生态恢复提供科学依据。