鄂尔多斯市防风固沙功能时空变化及驱动因素分析

马 苏，刘军会，康玉麟，刘 洋，郑颖娟，周甲男

中国环境科学研究院环境信息研究所，北京 100012

土壤风蚀是干旱与半干旱地区广泛存在的严重环境问题之一[1]. 土壤风蚀导致的土地沙化和土地退化降低了土壤的生产力，削弱了土壤保持服务功能[2-3]，威胁着区域经济发展、粮食安全和人类福祉[4-5]. 根据2015年发布的《中国荒漠化和沙化状况公报》，截至2014年，我国仍有荒漠化土地361.16×104km2和沙化土地172.12×104km2. 为遏制土地沙化与环境恶化的趋势，国家相继实施了一系列重大的生态保护与建设工程，如“三北”防护林、黄土高原中上游退耕还林还草、京津风沙源治理等，使得我国北方地区土壤风蚀得到了有效缓解[6-7]，但是生态系统的恢复是漫长复杂的过程，及时监测、评估工程区生态状况及功能变化尤为重要[8]. 而现有研究多侧重于植被覆盖状况的监测，对区域防风固沙功能的动态变化及其空间差异的研究相对较少.

防风固沙是生态系统通过其结构与过程减少由于风蚀所导致的土壤侵蚀的作用[9]，这种保持土壤、抑制风蚀过程的功能即为防风固沙功能[10-11]. 风蚀模型是评估防风固沙功能的主要技术手段，并已经发展演化出通用风蚀方程(WEQ)[12]、德克萨斯侵蚀分析模型(TEAM)[13]、Bocharov模型[14]、风蚀预报系统(WEPS)[15]和修正的风蚀模型(RWEQ)[16-17]等多种模型. 其中，RWEQ是目前应用最广泛的土壤风蚀模型[18-23]，该模型充分考虑气候变化、植被覆盖状况、土壤质地、地表粗糙度等因素，通过参数的修正和公式调整可以应用到我国北方土壤风蚀状况的评估[24-32].

鄂尔多斯市地处我国北方防沙屏障带的南侧，境内分布有库布齐沙漠和毛乌素沙地，土壤风蚀严重，生态环境脆弱，是京津地区主要风沙源之一，其防风固沙能力不仅影响当地生态屏障建设，而且关系到环京津地区的环境安全. 因此该研究基于长时间序列卫星遥感影像以及气象、植被、土壤等数据资料，采用修正的风蚀模型(RWEQ)与GIS技术，研究2000−2018年鄂尔多斯市防风固沙功能的时空变化规律，并从气候变化和人类活动角度探讨驱动影响机理，以期为鄂尔多斯市科学实施生态修复和筑牢生态安全屏障提供科学参考.

1 研究区概况

鄂 尔 多 斯 市 位于106°42′40″E~111°27′20″E、37°35′24″N~40°51′40″N之 间，总 面 积 约8.68×104km2(见图1). 地处黄河“几字弯”腹地，三面黄河环抱，平均海拔1 000~1 500 m，北部为黄河冲积平原，东部为丘陵沟壑水土流失区和砒砂岩裸露区，西北部为库布齐沙漠，东南部为毛乌素沙地，沙化土地约占全市总面积的48%. 气候类型为温带大陆性干旱半干旱气候，年降水量170~350 mm，年蒸发量高达2 000~3 000 mm；多大风，年均风速为3.6 m/s，最大风速可达22 m/s. 植被以典型草原、荒漠草原和草原化荒漠为主. 土壤类型主要包括栗钙土、棕钙土和风沙土等.

图 1 鄂尔多斯市的地理位置示意Fig.1 The location of Ordos

2 数据来源及处理

遥感数据：归一化植被指数(NDVI)为MOD13Q1产品，该数据已经经过几何精纠正、辐射校正、大气校正等预处理. 利用MRT工具进行投影和格式转换批处理，然后采用最大值合成法(maximum value composite)获得NDVI数据，并基于像元二分法生产2000−2018年的植被覆盖度.

土地覆被数据：使用2000年的Landsat-TM遥感影像数据和2018年的Landsat8遥感影像数据，结合高分遥感影像，采用CART决策树分类分别得到2000年和2018年的土地覆被图.

气象数据：来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn)以及鄂尔多斯市气象局8个气象站点的数据. 数据内容包括鄂尔多斯市各气象站点的编号、经纬度和海拔，以及每个气象站点在相应分析时间尺度内的降水量、风速、气温、日照时数等；采用Kriging对长时间序列离散点的气象数据进行插值处理，以得到气象要素的栅格图像.

土壤数据：来源于中国科学院南京土壤研究所提供的1:1 000 000土壤图及所附的土壤属性表，包括有机质、砂粒和黏粒含量等属性.

其他数据：雪盖数据采用中国科学院旱区寒区科学数据中心(http://westdc.westgis.ac.cn)的中国雪深长时间序列数据集计算得到. 数字高程(DEM)来自中国科学院资源环境科学数据中心(http://www.resdc.cn)，空间分辨率为90 m.

所有数据在ArcGIS中进行运算时投影坐标系统均为Krasovsky_1940_Albers，数据均重采样为250 m的空间精度.

3 研究方法

3.1 RWEQ模型

RWEQ模型是一种以较高时空分辨率对区域土壤风蚀状况进行长时间序列估算，从而有效预测防风固沙量的模型[33-36]. 因此，该研究采用RWEQ模型测算防风固沙量，进而评估鄂尔多斯市防风固沙功能及其动态变化，用单位面积防风固沙量的多少代表防风固沙服务功能的强弱. 防风固沙量(SRA)为潜在土壤风蚀量( SLp) 与实际土壤风蚀量(SLr)的差值，具体计算公式如下：

式中：SLp为潜在风蚀量，kg/m2；Z为下风向距离；sr为潜在关键地块长度，m；Qpmax为潜在风力的最大输沙能力，kg/m；sp为实际关键地块长度，m；WF为气象因子，kg/m；EF和SCF分别为土壤可蚀性因子和土壤结皮因子；K和C分别为土壤粗糙因子与植被因子；SLr为实际风蚀量，kg/m2；Qrmax为实际风力的最大输沙能力，kg/m；SRA表示单位面积防风固沙量，kg/m2.

a) 气象因子(WF). 气象因子表征了在考虑降雨、温度、日照及雪盖等因素下，风力对土壤颗粒的搬运能力.

式中：g为重力加速度，m/s2，取值9.8 m/s2；ρ为空气密度，kg/m3，该文取1.205 kg/m3；SW为土壤湿度因子；SD为雪盖因子(无积雪覆盖天数与研究总天数之比，雪盖深度<25.4 mm为无积雪覆盖)，结合鄂尔多斯市的地理位置以及雪盖数据的处理结果，该文取1；Ui为2 m处的风速，m/s；Ut为2 m处的临界风速，笔者所在课题组于2021年8−9月在伊金霍洛旗选取22个土壤剖面，采用TDR100时域反射仪测度土壤水分含量，发现土壤含水量均在1%左右，而刘玉章等[37]的风洞试验研究表明，起沙风速一般随沙中含水量的增加而增加，在含水量为0.8%左右时，起沙风速为4.5~5.0 m/s，所以鄂尔多斯市的临界风速取5 m/s；Nd为观测风速大于临界风速的天数，d；N为一次试验中观察的总次数，该研究中进行了1年累计，所以设定N为365.

b) 土壤可蚀性因子(EF)与土壤结皮因子(SCF).土壤可蚀性因子与土壤结皮因子的计算公式分别如式(9)(10)所示.

式中：sa为土壤粗砂含量(0.2~2 mm)，%；si为土壤粉砂含量，%；cl为土壤黏粒含量，%；[OM]为土壤有机质含量，%；[CaCO3]为碳酸钙含量，%.

c) 植被覆盖因子(C). 植被覆盖度是影响风蚀的关键因子，其覆盖程度直接影响近地表风速及土壤粗糙度.

式中：SC为植被覆盖度，%；NDVImin为裸地对应NDVI值，取累计频率为2%的值；NDVImax为NDVI最大值，取累计频率为98%的值.

d) 地表粗糙因子(K). 地表粗糙因子的计算公式如式(13)所示.

式中，slope为地形坡度，根据DEM利用ArcGIS软件中的Slope工具实现.

3.2 Theil-Sen趋势分析

该方法在栅格尺度上对防风固沙量的变化趋势进行模拟，能够反映研究时段内防风固沙量变化幅度的空间分布，其计算公式：

式中：SLOPE为回归方程的斜率；m为数据集时间段的年数；Fj为某像元第j年的防风固沙功能，t/km2.若SLOPE为正，表示随着时间变化，防风固沙功能呈上升趋势；反之，当SLOPE值为负数时，表示随着时间变化，防风固沙功能呈现下降趋势；SLOPE的数值大小反映防风固沙功能变化的程度.

3.3 Pearson相关分析

Pearson相关系数能够表征2个变量之间的相关程度，该研究以像元为计算单元，分析降水对防风固沙功能的影响，其计算公式为

式中：x为像元的防风固沙功能，t/km2；y为降水量，mm.

4 结果与讨论

4.1 防风固沙功能现状特征

为了更好地表达防风固沙功能的空间差异，依据研究区防风固沙功能值域及自然断裂法，设置防风固沙功能在0~0.7×104t/km2之间为低值区，在0.7×104~2.9×104t/km2之间为较低区，在2.9×104~5.0×104t/km2之间为一般区，在5.0×104~7.8×104t/km2之间为较高区，在7.8×104t/km2以上为高值区.

评估结果表明，2018年鄂尔多斯市防风固沙总量为95.07×108t，单位面积防风固沙量为10.95×104t/km2. 空间上，防风固沙功能存在明显地域差异(见图2)，高值区域主要分布在鄂尔多斯市东部和北部黄河沿岸，包括准格尔旗、达拉特旗东部和北部、东胜区、康巴什区、伊金霍洛旗东部等地区，面积占比为25.14%；较高区域零散分布在鄂尔多斯市中部和南部，包括鄂托克旗西北部、鄂托克前旗大部、乌审旗南部、达拉特旗南部等区域，面积占比为33.03%；一般区域集中分布在杭锦旗、鄂托克旗和乌审旗的部分区域，面积占比为15.58%；较低区和低值区主要分布在库布齐沙漠和毛乌素沙地，面积占比分别为15.27%、10.97%.

图 2 2018年鄂尔多斯市防风固沙功能的空间分布Fig.2 Spatial pattern of wind-breaking and sand-fixing function in the Ordos in 2018

4.2 防风固沙功能时空变化特征

从时间变化(见图3)上看，2000−2018年，鄂尔多斯市防风固沙功能呈显著增加趋势. 单位面积防风固沙量由1.62×104t/km2增至10.95×104t/ km2，以平均每年26.96%的变化率波动增加. 从不同阶段分析，2000−2009年，单位面积防风固沙量变化相对平缓，年均变化率为31.68%，期间防风固沙量增加了1.19×104t/km2；2009−2018年增长幅度较大，年均变化率为63.03%，期间单位面积防风固沙量增加了8.14×104t/km2.

图 3 2000−2018年鄂尔多斯市防风固沙功能的变化Fig.3 Annual changes of wind-breaking and sandfixing function in Ordos from 2000 to 2018

从空间变化(见图4)上看，2000−2018年鄂尔多斯市有54.89%的区域防风固沙功能增加. 其中，防风固沙功能明显增加区主要分布在准格尔旗和乌审旗的东南部等区域，其面积占鄂尔多斯市总面积的18.41%；防风固沙功能一般增加区主要分布在杭锦旗的库布齐沙漠北部、乌审旗的毛乌素沙地、鄂托克前旗、达拉特旗、伊金霍洛旗、东胜区、康巴什区等区域，面积约占鄂尔多斯市总面积的36.48%. 此外，相比2000年，2018年鄂尔多斯市有15.33%的区域防风固沙功能下降，主要分布在杭锦旗北部的库布齐沙漠和鄂托克旗西部. 其余29.76%的防风固沙功能没有发生明显变化.

图 4 鄂尔多斯市2000—2018年防风固沙功能的空间分布Fig.4 The spatial distribution pattern of wind-breaking and sand-fixation function in Ordos from 2000 to 2018

4.3 防风固沙功能驱动因素分析

4.3.1降雨量对防风固沙功能的影响

利用R语言编程进行降雨量和防风固沙功能的显著性检验，结果表明年降雨量与防风固沙功能存在显著相关(P<0.01，r=0.69). 2000−2018年鄂尔多斯市的降水量在200.58~430.96 mm之间，多年平均降雨量为317.3 mm，呈现波动增加的趋势，且与防风固沙能力变化趋势呈现较好的一致性. 从防风固沙能力变化较大的年份分析，2012年降雨量较2011年大幅增加，2012年防风固沙功能也显著增加；而2017年降雨量下降，防风固沙功能也小幅下降(见图5).

图 5 2000—2018年鄂尔多斯市防风固沙功能和降雨量的变化Fig.5 Annual changes of rainfall and WBSF function in Ordos from 2000 to 2018

图 6 鄂尔多斯市降雨量与防风固沙功能相关系数的空间分布Fig.6 The Spatial pattern of correlation coefficients between rainfall and WBSF function in Ordos

空间上，78.59%的区域相关系数为0.4~0.9(见图6)，主要分布在降雨量≥200 mm区域，该区域为防风固沙功能一般区及以上，而相关系数为负值的占比较小，主要分布在库布齐沙漠，该区域为防风固沙功能低值区. 说明该区域降雨量对防风固沙功能的变化产生了积极影响.

4.3.2土地利用对防风固沙功能的影响

根据鄂尔多斯市土地利用类型转移矩阵(见表1)，2000−2018年，鄂尔多斯土地利用类型变化的主要特点是城市扩张、林地和草地增加以及沙地、耕地与湿地资源减少. 林地、草地、建设用地表现出不同程度的增加，其增幅分别为11.03%、1.45%、163.18%，其中林地和草地面积的增加主要来源于沙地，沙地整体呈现出净转出的特征，建设用地的增加主要来源于草地. 耕地和湿地总面积减少，其减幅分别为2.09%和6.71%，其中耕地主要转换为建设用地以及退耕还草.

表 1 鄂尔多斯市2000−2018年土地利用类型转移矩阵Table 1 Transferring matrix of land use type in Ordos from 2000 to 2018 km2

通过对2000−2018年土地利用类型变化和防风固沙功能变化的叠加(见图7)，可统计分析不同土地利用类型之间的转化对防风固沙功能的影响. 结果表明，鄂尔多斯市土地利用类型发生变化的面积占研究区总面积的25.62%，而在土地利用类型发生变化的区域，其防风固沙量的变化占研究区总防风固沙量变化的29.86%，说明人类活动主导下的土地利用类型变化对防风固沙功能的影响不容忽视. 影响防风固沙功能的用地类型主要在林地、草地与耕地和沙地之间进行. 鄂尔多斯市近20年实施了大面积植树造林、京津风沙源治理工程以及禁牧、休牧等措施，2018年沙地面积较2000年缩减了1 378.56 km2，减幅为6.51%，沙地主要转变为草地、林地和耕地，这部分区域防风固沙功能增加了35.00%. 草地和林地增加区域防风固沙功能分别增加了20.38%和8.78%. 而耕地转为林地，防风固沙功能降低了2.47%，耕地向林地转换大多是人为造林的结果，新造林地在生长初期植被覆盖度较低，导致短期内地表植被的保护作用下降，固沙功能轻度降低. 综合来说，与生态环境质量提高相关的土地利用类型变化对区域防风固沙功能起到了一定的增强作用.

图 7 2000—2018年鄂尔多斯市土地利用类型和防风固沙功能的空间变化Fig.7 The spatial change of land use and WBSF function in Ordos from 2000 to 2018

4.4 RWEQ模型适用性讨论

该研究评估了2000−2018年鄂尔多斯市防风固沙功能变化及长时间序列的空间差异，发现鄂尔多斯市防风固沙功能明显提升，这与张照营[38]和张燕婷[39]得到的内蒙古自治区防沙屏障带和北方防沙带的防风固沙能力明显增强的结论均一致. 该研究植被覆盖度的变化与邵艳莹等[40-42]的研究结果一致. 该研究计算得出的防风固沙能力稍高于江凌等[43]估算的2000−2010年内蒙古自治区的年均防风固沙能力(4 880 t/km2)，原因在于评估时段和区域不同. 另外，笔者采用小时数据作为气象因子计算的数据来源，时间分辨率更高.

RWEQ模型是以物理过程理论为基础的经验模型[44]，其计算参数具有较强的地域性，具体应用时需进行必要的参数修正. 但迟文峰等[45]利用内蒙古高原137Cs示踪技术监测成果验证RWEQ模型反演结果，发现二者具有较好的拟合性(R2=0.83,P<0.01)，从侧面说明RWEQ模型在我国北方区域具有可行性. 笔者所在课题组通过在伊金霍洛旗实地采样，测量22个土壤剖面的土壤含水量来修正当地起沙风速，其他参数在综合考虑与借鉴前人研究成果的基础上，采用了我国北方沙化地区修正后的关键参数和计算公式，并且为满足模型对数据空间和时间分辨率的要求，对输入参数进行了插值、拟合、渐进等优化，提高了模型的精确度. 在今后的研究中，可进一步对模型参数如土地粗糙度以及土壤结皮和可蚀性因子的土壤粒径转换等问题进行本地化调整，构建具有更广泛适用性的土壤风蚀模型，以便更准确地反映区域防风固沙功能的变化状况，进而引导防风固沙生态治理工程有针对性地投入.

5 结论

a) 2018年鄂尔多斯市防风固沙总量为95.07×108t，单位面积防风固沙量为10.95×104t/km2；防风固沙功能明显存在区域差异，整体呈现由东向西递减的趋势，高值区面积占比为25.14%，低值区面积占比为10.97%.

b) 2000−2018年，鄂尔多斯市防风固沙功能以平均每年26.96%的变化率波动增加. 54.89%的区域防风固沙功能增强，主要分布在准格尔旗和乌审旗的东南部；29.76%的区域防风固沙功能稳定，主要位于杭锦旗、鄂托克旗、伊金霍洛旗的部分区域；而防风固沙功能降低区域为15.33%，主要位于杭锦旗北部的库布齐沙漠和鄂托克旗西部地区.

c) 在驱动因素中，鄂尔多斯防风固沙功能变化与降雨量呈显著相关(P<0.01，R2=0.69)，多年平均降雨量为317.3 mm，呈现波动增加的趋势，且与防风固沙能力的变化趋势呈现较好的一致性，78.59%的区域二者相关系数为0.4~0.9. 土地利用类型变化中，沙地面积减幅为6.51%，其防风固沙功能增加了35.00%；草地和林地增加区域防风固沙功能分别增加了20.38%和8.78%，退耕还林还草、沙地改善等有利于生态环境质量提高的土地利用类型转化，对改善区域防风固沙功能有积极作用.

d) 生态保护和恢复工程建设应该综合考虑气候变化和人类活动的复合影响，发挥气候变化带来的正面效应，结合鄂尔多斯市的自然地理环境特点，注重生态工程布局的区域适宜性，对生态系统进行科学的保护与修复，持续发挥北方防沙带生态安全屏障作用.