排土场植被覆盖度及地形变化对流域水土流失的影响

2022-06-06雷少刚

煤炭科学技术 2022年4期

邓彪，陈航，李恒，雷少刚

(1.神华北电胜利能源有限公司，内蒙古锡林浩特 026015；2.中国矿业大学矿山生态修复教育部工程研究中心，江苏徐州 221116；3.浙江大学公共管理学院，浙江杭州 310058)

0 引言

21世纪以前，我国采煤方式以井工为主，进入本世纪后露天煤矿发展迅速，产量已居世界第二，我国已成为露天采煤大国和强国，与此同时，伴随而来的一系列环境问题日益突出[1]。由地下岩土堆积而成的外排土场破坏了原地表植被，并使得矿区水土流失、土地荒漠化现象加剧，生态环境严重恶化[2-3]。矿区作为一种人为作用诱发产生的特殊侵蚀区域，其严重的水土流失问题越来越受到国家、地方政府以及广大工作者的高度重视，已成为制约中国北方典型草原区实施可持续发展战略的一个热点问题[4-5]。影响水土流失的因子主要包括气候、地形、土壤以及地表覆被等因素[6]。区域气候因素人为无法改变；占地面积和经济条件的限制也使得排土场在短期内难以降低边坡角度，无法从结构上优化稳定性；虽然排土场边坡和平台覆有地表土，理化性质的差异也使得不同土壤之间难以紧密贴合，在雨水渗透的条件下，表层土体易发生滑动；因此在人为因素下，只能通过改变地表覆被来缓解坡面水土流失[7-8]。但由于我国露天矿多处于黄土高原和内蒙古高原等干旱、半干旱的生态脆弱区，水资源匮乏，盲目重建植被浪费物力财力，如何通过排土场植被重建缓解水土流失、实现效益最大化逐渐成为研究热点。

关于地表覆被与水土流失的关系，各领域学者分别从斑块、坡面以及流域尺度上做了相关研究[9-10]。在斑块尺度上，张思毅等[11]通过室内模拟，研究不同覆盖度的植被在不同雨强下的减沙效应及水动力学机理，结果表明坡面的达西-韦斯巴赫阻力系数随着植被覆盖度呈指数函数增加，植被覆盖度为88%时的土壤临界剪切力约是裸地的2.5倍。甘卓亭等[12]通过人工模拟降雨，定量研究不同生长阶段的植被对坡面侵蚀产沙的贡献率，试验表明不同生长阶段植被的减流效果依赖于不同植被器官作用，并且不同生长阶段的作用效果也不相同。坡面尺度上，陈立顶等[13]利用同位素示踪法针对不同坡面类型(坡长、坡度、坡形)评估土壤侵蚀效应，研究发现坡面形态对土壤侵蚀具有重要影响，并且不同植被空间配置对土壤侵蚀具有明显差异。史倩华等[14]测定不同年限和植被配置的排土场侵蚀状态指标，通过结果对比表明合理的植被配置在改善土壤性质、控制侵蚀状态等方面均有较好效果。流域尺度上，聂小军等[15]利用同位素示踪法分析了未采区、自恢复沉陷区与植被修复区的土壤侵蚀与养分特征，研究了矿区土壤侵蚀与养分的演变规律。索安宁等[16]基于流域降水与水土流失结果的统计，分析得出水土流失变化主要由地表形态引起的，降水量变化的贡献很小，农田水利建设可能是导致水土流失变化的主要原因。通过以上研究发现，前人研究多集中在斑块与坡面等中小尺度，研究流域大尺度多从特定气候和地貌覆被条件入手，依赖遥感监测手段，从模型模拟方面开展。但总体来讲，大尺度水土流失预测仍然是个难点，需要进一步深入探讨。

笔者以矿群为研究尺度，在地形、土壤、气候和土地利用数据确定的条件下，分析了影响锡林浩特矿群所在流域水土流失量增加的主要原因，在此基础上分析了地表覆被和地形变化对水土流失模拟结果的影响，通过设置排土场不同植被覆盖度，分区域模拟矿群所在流域水土流失，分析变化规律，并结合相关数学模型计算缓解水土流失的最佳植被覆盖度范围，以此来指导矿区植被重建工作，实现效益最大化。

1 研究区概况

本研究区为内蒙古自治区锡林郭勒盟锡林浩特市周边胜利矿区群(图1)，包括胜利一号露天矿、西二矿、西三矿，胜利东二号露天煤矿以及乌兰图嘎锗矿。截至2016年底各矿区共有9座外排土场，占地面积约27.89 km2。研究区地处温带丛生禾草典型腹地，除河滩、丘间洼地和盐化湖盆地外均为典型草原，地理坐标为43.90°～44.23°N，115.41°～116.44°E，海拔970～1 212 m。该区属于半干旱草原气候，极端最高气温38.3 ℃，极端最低气温-42.4 ℃，平均气温1.7 ℃；年最大降雨量481.0 mm，年最小降雨量146.7 mm，年平均降水量294.74 mm，年均蒸发量1 794.64 mm。春季多风，风速2.1～8.4 m/s,年均风速3.5 m/s,瞬时最大风速36.6 m/s，最大冻土深度2.89 m，年无霜期122 d。研究区土壤类型主要由栗钙土、草甸栗钙土、草甸土等组成，部分区域由于草场退化而形成沙化、砾石化栗钙土，土壤有机质含量低，土壤肥力差。研究区植被以芨芨草、针茅、羊草、猪毛菜为主。矿区排土场复垦始于2006年，在堆积岩土上覆以采坑剥离表土，再配以适生灌草和浇灌设施，以保证植被生长所需的基本条件。绿化植被主要有沙打旺、紫花苜蓿、柠条、油蒿、披碱草等。但由于表土堆积时人工作业的不均匀性，使得排土场土壤异质性差别较大，同时由于矿区地处生态脆弱区，排土场水土流失较为严重。

图1 研究区区位

2 水土流失对植被覆盖度和地形变氏响应的研究方法

2.1 GeoWEEP模型

通过总结国内外土壤侵蚀模型的特点，选用GeoWEPP模型研究矿区排土场水土流失的变化规律。土壤侵蚀研究经历了由经验模型、基于物理机制的土壤侵蚀模型，向具有灵活可视化地理用户界面的新一代物理模型发展的过程，模型的发展为流域水土保持措施优化提供了更为有效的评估工具[17]。经验模型结构形式往往比较简单，可在宏观上对土壤侵蚀进行评价，但难以精细动态模拟土壤侵蚀的时空变化过程，且缺乏对侵蚀过程及其机理的深入研究，而采用物理过程模型探索土壤侵蚀是当前的研究热点[18]。GeoWEPP模型由WEPP模型改进而来，是具有基于物理过程和连续模拟优势的分布式水文模型，该模型将GIS与WEPP模型结合起来，既可以反映自然界中下垫面的复杂性、降水和侵蚀产沙的时空不均匀性，也可以描述侵蚀的动态变化和泥沙输移、沉降过程的时空变化(图2)；不仅考虑研究区域的地形地貌，气候，土壤与土地利用方式等条件，同时提高了模型的模拟效率[19]。

图2 GeoWEPP模型技术框架

利用GeoWEPP模型，通过输入开采前后的地形数据。及研究区的土壤、气候及土地利用数据，在不同植被覆盖度下对研究区进行为期10 a的水土流失模拟。

2.2 数据准备

2.2.1 地形数据

地形数据是模拟水土流失的基础数据，采前地形数据来源于DLR的数字高程数据(DEM)，由美国太空总署NASA和国防部国家测绘局NIMA2000年联合测量的SRTM数据。采后地形数据则根据2016年底在矿区所测高程点经克里金插值得到，其中因开采活动而形成的排土场和采坑等地形变化区域，根据各矿区的生产报告和开采计划，结合RTK测量仪进行实地高程测量，将其处理成符合现状的地形数据并镶嵌到初始高程数据上作为现状地形数据。考虑研究区数据量较大，因此需要将地形数据重采样到30 m×30 m，并在进行排土场不同植被覆盖度下水土流失模拟时分别对胜利东二号矿区流域和西北各矿区流域进行模拟，以确保模拟结果的相对准确性。

2.2.2 土壤数据

根据内蒙古土壤类型划分确定研究区土壤类型，包括栗钙土、潮土、沼泽土和草甸土4种，其中沼泽土和草甸土为主要土壤类型。通过野外采集样本分析获取土壤的粒径组成、有机质含量等信息，按照土壤类型统计相关指标，计算每种类型各指标的平均值，建立土壤数据库。研究区主要土壤参数见表1。

表1 研究区土壤参数

2.2.3 气候数据

由国家气象网站下载锡林浩特市2011—2016年的气候数据，统计月平均降雨量，月最高和最低气温，降雨天数等数据，在多年气象资料统计参数的基础上，经气候发生器CLIGEN生成模型所需要的气候数据。

2.2.4 土地利用数据

研究区地类主要包括湿地、裸地建设用地、草地和排土场植被，分别由所用地形数据对应的2000年和2016年夏季Landsat影像分类生成。其中自然草原区植被类型为芨芨草、针茅、羊草，排土场植被类型为沙打旺、紫花苜蓿、柠条。将自然草原区和排土场植被分别对应植被数据库中的不同植被类型，并将草地和排土场植被分别按照植被覆盖度从10%～90%每隔10%设置9个等级，据此进行水土流失模拟效果的对比分析。

2.3 非线性均值变点分析

为获取研究区植被覆盖度的最佳阈值，采用非线性均值变点法进行不同植被覆盖度下水土流失模拟结果的分析。均值变点法是一种最常见最直观的数学算法，它以平均值为分析对象，通过计算整个样本数据的离差平方和SS与分段样本的统计量Si之差来确定变点，变点的存在会使S和Si的差距增大，该方法对恰有一个变点的检验最为有效[20]，多个变点时有可能因为均值的多次升降反而抵消了S和Si之间的差距。均值变点法的基本原理如下：

1)令i=2，3，…，n，对每个i将样本分为2段：X1，X2，…，Xi-1和Xi，Xi+1，…，Xn，计算每段样本的算术平均值Xi1和Xi2及统计量：

(1)

(2)

(3)

3)计算期望值E：

E=SS-Si,i=2,3,…,n

(4)

3 结果与分析

3.1 地形格局变化因素分析及水土流失模拟

为分析研究区地形变化的主要原因，分别对研究区采前地形与现状地形的高程和坡度进行统计，统计结果如图3所示。由图3可看出，现状地形下采坑和排土场的存在使得研究区高程分布变广，排土场斜坡-平台的堆积形态使平台高程点占比有明显突起，突起间隔与排土场平台间隔相关；采前状态下研究区整体坡度呈西高东低状态，地势起伏变化不明显，而矿坑挖损和排土场堆积增加了相应区域的整体坡度，相比于采前地形，露天矿业的发展使当地地貌发生明显变化，这或许是导致流域水土流失增加的主要原因。

图3 采前和现状条件下的高程坡度变化

为分析植被类型变化与地形变化对流域水土流失结果的影响区别，本研究通过控制变量法，利用三组模拟试验进行了地形格局变化下流域水土流失效果的对比。第1组在植被类型为自然草原区植被的条件下，对采前地形进行为期30 a的水土流失模拟；第2组在地形不变的前提下，对采前地形的植被类型进行替换，即将排土场在采前地形中对应区域的自然草原区植被类型变为排土场植被类型进行模拟，研究排土场植被类型变化对流域水土流失结果的影响；第3组在植被类型不变的前提下，将采前地形替换为现状地形，研究排土场地形变化对流域水土流失结果的影响，模拟结果见表2。

表2 地形格局变化模拟结果

将植被类型变化和地形变化模拟结果分别与采前模拟做差值，如图4所示，图4中NAN指其他情况。

图4 地形格局变化水土流失模拟结果

从水土流失差值图中可以发现，在地形不变的情况下，仅仅改变排土场在采前地形中对应区域的自然草原区植被类型对区域水土流失效果变化的影响很小，流失总量上看植被类型变化后的研究区水土流失量仅低于采前地形水土流失量0.2%，而沉积值略高，这说明排土场重建植被相较于自然草原区植被的水土保持效果更好一些；而地形发生改变时模拟的水土流失量高于采前地形水土流失量4%，这说明地形变化带来的水土流失影响效果更大。纵向来看，排土场重建植被以沙打旺、柠条为主，一方面借助植被根系稳固边坡结构，另一方面通过改善土壤性质缓解水土流失，但其缓解效果明显弱于地形变化带来的水土流失增量，这充分说明，仅改变排土场植被类型对地形变化带来的水土流失影响的缓解程度并不高，需要结合对排土场的地貌重塑以控制水土流失。

通过分析水土流失图、沉积图(图5)与采坑和排土场的位置关系可以看出，土壤侵蚀主要发生在采坑边缘以及排土场边坡处，排土场平台侵蚀效果不明显，河网形成过程中对采坑的填洼使得采坑形成侵蚀沉积区，排土场与自然地貌的衔接处也有大量沉积产生。由此可知，排土场边坡为研究区水土流失的主要源地。因此，单一方面改变矿区植被类型不能有效缓解地形变化后的水土流失效果，应该在地形重塑的基础上进行植被恢复，从而达到研究区生态环境优化的效益最大化。

图5 地形格局变化土壤沉积模拟结果

3.2 不同植被覆盖度下排土场的水土流失模拟

3.2.1 模拟结果分析

由上述分析可知，地形变化前后研究区水土流失变化较大的区域主要分布在排土场周围。为分析排土场不同植被覆盖度情况下流域水土流失的模拟效果，本次研究利用GeoWEPP模型，以50 a为模拟时长，通过输入土壤、气候和土地利用数据，在10%的植被覆盖度间隔下，分别对胜利东二号矿区和西南各矿区所在流域进行了9组排土场水土流失模拟试验，代表现状地形条件下排土场所在流域在植被覆盖度处于10%～90%条件下的水土流失状况。模拟结果如图6、图7所示。

图6 胜利东二号矿区流域不同植被覆盖度下水土流失模拟结果

为了保证模拟结果的准确性，该模型在数据选择上通过专业部门采集和实地采样分析得到。但由于特定流域水土流失精确数据难以获取，模型模拟的结果缺少实测数据校正，因此通过对模拟结果进行无量纲处理，研究排土场植被覆盖度对流域水土流失的影响，确保数据间对比的合理性。无量纲处理通过以下公式进行：

Pi=(Ai-Amin)/(Amax-Amin)

(2)

其中，Ai为第i个模拟情况下的水土流失模拟结果；Amin为最小模拟结果；Amax为最大模拟结果。分别统计2个流域的水土流失量值，及经归一化处理的水土流失模拟结果，结果见表3(其中由于90%植被覆盖度下的归一化值为0，为获取拟合曲线将其舍去)。

表3 排土场不同植被覆盖度下研究区的水土流失模拟结果

以排土场植被覆盖度为自变量，水土流失归一化值为因变量分别进行拟合，水土流失拟合结果如图8所示。

由拟合曲线可知：随着排土场植被覆盖度的增加，对应水土流失归一化值逐渐减小，两者拟合指数函数效果较好，决定系数R2均高于92%。变化曲线说明随着排土场植被覆盖度的提高，流域内水土流失情况有一定程度上的缓解，并且植被覆盖度越高，效果越明显。考虑到现实状态下的植被覆盖度不会达到试验水平，并且拟合函数逐渐变缓也说明随着植被覆盖度的提高，投入物力财力带来的边际生态效益到达某一变点后开始下降，符合边际效用递减规律。因此需要确定合适的植被覆盖度作为排土场植被重建的目标，以实现植被重建效益最大化。由于拟合函数的特殊性，直观手段不能取得最佳阈值，需要通过均值变点法计算获取。

3.2.2 均值变点法求取植被覆盖度最佳范围

由图8可知，水土流失归一化拟合曲线呈幂函数分布，满足非线性均值变点法分析应用要求，利用式(1)—式(4)分别对胜利东二号矿区流域、西北各矿区流域的水土流失归一化值的均值差E进行计算，据此绘制均值变点统计分析图(图9)。

图8 水土流失量归一化拟合结果

从图9中可以看出，水土流失归一化值的均值差E随植被覆盖度增加呈先升高后下降的变化，这表示在变点值之后，随着排土场植被覆盖度的增加，投入所带来的边际生态效益逐渐减少，符合边际效用递减规律。对胜利东二矿区流域而言，变点出现在植被覆盖度40%处；而对西北各矿区流域和矿群总体统计结果而言，变点出现在植被覆盖度50%处。考虑到实际植被重建施工过程中的合理性和可操作性，认为锡林浩特矿群排土场的植被覆盖度最佳为40%～50%。