基于自然边坡模型的草原煤矿排土场坡形重塑

李 恒，雷少刚，黄云鑫，程 伟，邓 彪，张周爱

(1.中国矿业大学 矿山生态修复教育部工程研究中心，江苏 徐州 221116; 2.神华北电胜利能源有限公司，内蒙古 锡林浩特 026015; 3.神华宝日希勒能源有限公司，内蒙古 呼伦贝尔 021500)

采矿活动会对生态环境造成极大的破坏，主要体现在改变了矿区原有地表形态和自然景观，从而引发一系列的土地与环境问题，因此将受采矿扰动的区域一定程度地恢复到扰动前的生态服务功能，使之能与周围的自然环境相协调是十分重要的[1-2]。在土地复垦与生态重建过程中，地形重塑、土壤重构和植被重建被认为是3个十分重要的环节，其中地形重塑为基础[3]。但以往的土地复垦方案更重视土壤重构[4-5]和植被重建[6]方面，而地形重塑受制于理论和技术层面的问题，往往会被研究者所轻视。在这种情况下，排土场往往被重塑为简单的阶梯式地形，未能与当地自然环境相融合，不仅养护成本较高，且缺乏长期稳定性，从而产生较为严重的水土流失和景观破碎现象。因此寻找合理的地形重塑理论依据和方法显得十分重要。

近自然地形重塑作为一种仿造当地自然形态进行地形重塑模型构建的理论，在矿区复垦活动中被国内外诸多学者应用。国外学者和研究单位在结合水文地质学[7-8]、自然地理学[9]、系统动力学等相关理论和原则的基础上，以塑造不会轻易变动，且与周围未损毁土地相适应的近自然地形为目的[10-11]，从流域角度开发出了GeoFluv[12]，GeoWEPP[13-15]等模型进行废弃矿区的地形重塑。国内学者也逐渐意识到近自然地形重塑的优势，并应用河流动力学[16-17]和流域地貌学[18-20]的相关理论和分析重新构建项目区的流域格局。现有近自然地形重塑的相关研究主要着眼于以水文为基础的损毁地形整体构建[21]，在排土场等工地形的坡面形态方面研究较少。

近自然地形重塑的主要理论依据是“师法自然”，即学习受当地气候、土壤、植被、水文条件影响下的自然地表形态，并应用其重塑矿区损毁地形。笔者以内蒙古胜利矿区为研究区，为设计出适宜胜利矿区自然条件的排土场边坡模型，将研究区内未受采矿扰动的自然地表形态的关键特征参数化并进行统计分析，分析其特征参数之间的相关性，在此基础上构建出与当地环境相适应，土壤损失量较小的近自然边坡模型。

1 研究区概况

研究区胜利矿区(图1)位于内蒙古自治区锡林浩特市西北部，锡林河以西的缓坡丘陵干草原与河谷冲积、湖积平原的过渡带,西部属缓坡丘陵,东部为河谷、湖积平原。气候类型属中温带半干旱、干旱大陆性季风气候[22]。胜利矿区排土场共5层台阶，每层台阶高15 m，边坡角为33°，排土场整体坡角为17°。由于排土场简单的采用阶梯式方法构建，长期受到雨水的冲刷侵蚀，土壤受侵蚀强度远大于周边自然环境，由此导致排土场边坡已经出现沟蚀、垮塌等现象。

图1 研究区区位Fig.1 Location of the study area

2 研究方法

2.1 DEM数据获取

数据获取时间为4月，进行边坡特征参数提取的参照区位于胜利矿区内的自然区域，其中陡坡提取区域位于胜利矿区西北部，缓坡提取区域位于胜利矿区西部。在研究区未受采矿扰动的自然区域内，使用大疆精灵3无人机获取高精度的地面数字高程模型(DEM)数据。数据范围为3 000 m×3 000 m，像元大小为0.068 m×0.068 m。

2.2 边坡特征参数提取及边坡模型构建

在Arcmap10.2中对DEM数据进行填洼后进行坡向分析和获取等高线两部分处理，进行这两部分处理的原因是坡向分析后的图像可较为清晰的显示出分水岭，而等高线分布可以显示出参照区域的坡度。以此为基准获取自然地形部分的边坡剖面线，通过边坡剖面线提取可知胜利矿区自然地形部分的主要边坡类型为反S型边坡，共提取出52组反S型边坡剖面线，其中包括28组陡坡部分的边坡剖面线及24组缓坡部分的边坡剖面线。

反S型边坡以拐点为分界处分为上侧的凸起部分和下侧的凹陷部分。反S型边坡与其他类型边坡有部分共同的边坡特征参数，包括坡高、坡长、凸面曲率和凹面曲率，但反S型边坡与其他类型边坡的不同之处在于其凸起部分和凹陷部分存在拐点，因此这两部分分别在水平面和铅垂面上有不同的占比。为此笔者提出了凸面水平占比和凸面竖直占比两个新的边坡特征参数，通过这两个参数可以确定拐点在边坡上的位置，进而可以确定反S型边坡的具体形状。

根据剖面线数据计算反S型边坡的特征参数，包括坡高、坡长、凸面曲率、凹面曲率、凸面水平占比、凸面垂直占比。其中，坡高指边坡坡面在铅垂面上的投影距离，坡长指边坡坡面在水平面上的投影距离，凸/凹面曲率指边坡凸起/凹陷部分形成的圆弧的曲率，凸面水平占比指边坡上侧凸起部分在水平面上的投影长度占坡长的百分比，凸面竖直占比指边坡凸起部分在铅垂面上的投影长度占坡高的比例，根据这些边坡特征参数可确定反S型边坡的坡形。分析边坡特征参数之间的相关性，利用统计学方法构建胜利矿区自然边坡的边坡模型。

2.3 边坡土壤损失量对比

利用土壤水蚀预报模型WEPP对构建的自然边坡模型和胜利矿区原排土场边坡模型进行土壤水蚀分析。WEPP模型在进行模拟时需要气象、土壤及植被覆盖情况等基础数据。其中气象数据以研究区2016年的气象数据为基础，使用气候发生器CLIGEN生成模型所需要的气象数据文件。土壤数据包括沙粒含量、黏粒含量、有机质含量、石砾含量、阳离子交换量及土壤容重，将这些数据输入WEPP中即可生成土壤数据文件。由于本次研究主要考虑反S型边坡模型与原排土场边坡模型之间土壤损失量的差异，为控制变量先假设边坡坡面的植被覆盖程度为0，即全部为裸地的情况。在确定了气象、土壤、植被覆盖情况等基础数据后，通过对比不同坡高、不同模拟年份条件下反S型边坡模型和原排土场模型的土壤损失量和土壤流失曲线图，以验证构建的近自然边坡模型的适用性。

3 研究结果

3.1 边坡特征参数相关性分析

为得到胜利矿区自然边坡特征参数之间的拟合公式，首先需要分析这些边坡特征参数之间的相关性。在利用SPSS进行相关性分析之前，由于选取的边坡特征参数有可能存在偏离一般规律的数值，因此首先要对偏离值进行去除。在去除了偏离值之后，根据坡高和坡长的关系制作散点图(图2)，根据图像可看出，胜利矿区缓坡部分与陡坡部分有着不同的规律，因此需要对两部分分别进行分析。

图2 坡高对坡长的散点Fig.2 Scatter plot and fitting curves of slope height to slope length

陡坡部分的边坡剖面线坡度在22%～46%，对陡坡部分的边坡特征参数进行相关性分析的结果发现，坡高与坡长、凸面水平占比与凸面竖直占比在0.01水平上显著相关;坡长与凸面曲率、坡长与凹面曲率在0.05水平上显著相关。

缓坡部分的边坡剖面线坡度在6%～10%，对缓坡部分的边坡特征参数进行相关性分析的结果发现，坡高与坡长、凸面曲率与凹面曲率、凸面水平占比与凹面竖直占比分别在0.01水平上显著相关。

3.2 边坡模型构建

3.2.1陡坡部分边坡模型构建

根据陡坡部分边坡特征参数之间的相关性分析，选取坡高对坡长、坡长对凸面曲率和凹面曲率、凸面水平占比对凸面竖直占比这4组相关特征参数进行曲线拟合，通过这种拟合得到的公式组，可以得到一个以坡高为自变量的近自然边坡模型。在进行曲线拟合之前，为检验模型的准确性，将提取两组数据作为验证数据。在曲线拟合的过程中(图3)，通过分析R2值，参数的T检验和方差的F检验后，分别得到了4组相关特征参数的拟合式:

s=e3.829+0.024d

(1)

qt=e2.530-0.007s

(2)

qa=e2.039-0.011s

(3)

zt=e-2.219+2.673pt

(4)

式中，d为坡高;s为坡长;qt为凸面曲率;qa为凹面曲率;pt为凸面水平占比;zt为凸面竖直占比。

图3 陡坡部分边坡特征参数拟合Fig.3 Slope characteristic parameter fitting diagram of steep part

其中由于凸面水平占比和凸面竖直占比与其他边坡特征参数没有相关性，而最终的边坡模型只能有坡高一个自变量，为获取特征规律应尽量采取与自然边坡特征参数相近的数值，因此求取实测数据凸面水平占比的均值为46.82%，根据式(4)计算凸面竖直占比的参考值为38.01%。由此构建出陡坡部分的自然边坡模型。将验证数据的坡高代入模型中，所得到的坡长、凸面曲率和凹面曲率数据与实测数据的相差幅度均低于10%，可说明模型的准确性。

3.2.2缓坡部分边坡模型构建

根据缓坡部分边坡特征参数的相关性分析，分别选取坡高对坡长、凸面曲率对凹面曲率、凸面水平占比对凹面竖直占比这3组相关特征参数进行曲线拟合(图4)。同样提取出两组数据作为验证，用上述方法拟合出了3组特征参数的拟合式:

s=e3.782+0.103d

(5)

qa=e0.24-0.443qt

(6)

zt=e-1.187+0.736pt

(7)

同样对实测数据的凸面曲率和凸面水平占比求取均值，并按照式(6)，(7)计算凹面曲率和凸面竖直占比的数值。经计算的凸面曲率参考值为1.41×10-3，凹面曲率为0.69×10-3，凸面水平占比为52.39%，凸面竖直占比为45.26%。

3.3 反S型边坡模型与原排土场边坡模型对比

考虑到实际情况，由于在进行排土场边坡重塑时不可能按照缓坡的坡度进行放坡处理，而陡坡部分边坡模型的坡度与原排土场边坡相近，为此在对比分析时使用陡坡部分边坡模型与原排土场边坡模型进行分析。在进行水蚀分析的时候，为确保分析结果的准确性，应当对除坡型因素外的其他变量进行控制，主要包括边坡高度，土壤类型，气候数据，植被覆盖以及模拟年份。为了使分析结果更加符合原排土场的实际情况，在进行水蚀分析模拟时，应当以一个或者数个台阶的高度作为基准考虑。在综合考虑边坡设计需要的基础上，笔者最终选取边坡高度分别为15，30，45和60 m来构建不同的边坡模型。而在模拟年份这一变量中，由于边坡形状或许会因为常年受到侵蚀而发生改变，导致土壤损失量发生变化，因此需要划分不同的模拟年份。本次分析将以1，10和50 a分别进行水蚀分析。分析结果见表1。

图4 缓坡部分边坡特征参数拟合Fig.4 Slope characteristic parameter fitting diagram of gentle part

表1 不同边坡的特征参数及年均土壤损失量
Table 1 Characteristic parameters annual soil loss of each slope

参数原排土场边坡反S型边坡坡高/m1530456015304560坡长/m56.86113.72170.58227.4465.9694.54135.50194.22边坡特征参数凸面曲率/10-3————7.916.484.863.22凹面曲率/10-3————3.722.721.730.91凸面水平占比/%————46.8246.8246.8246.82凸面竖直占比/%————38.0138.0138.0138.011 a23.9348.5072.1495.3917.0031.8642.2750.64土壤损失量/(kg·m-2)10 a32.0765.4696.40127.9021.6741.1854.0559.9150 a35.0071.78106.30141.2023.3544.3258.7065.65

由表1可知，在模拟坡高为15 m时，反S型边坡1，10，50 a的土壤损失量分别为原排土场边坡的71.04%,67.57%,66.71%；在模拟坡高为30 m时，反S型边坡1,10,50 a的土壤损失量分别为原排土场边坡的65.69%,62.91%,61.74%；在模拟坡高为45 m时，反S型边坡1,10,50 a的土壤损失量分别为原排土场边坡的58.59%,56.07%,55.22%；在模拟坡高为60 m时，反S型边坡1,10,50 a的土壤损失量分别为原排土场边坡的53.09%,46.84%,46.49%。

在对各个坡型的年均土壤损失量进行了对比分析后，为了对不同坡型土壤流失的分布情况进行对比，还需要根据两种坡型的土壤流失曲线图进行分析，图5分别为模拟年份为1 a，边坡高度为30 m的条件下两种坡型的土壤流失曲线图。其中红色的曲线表示坡面剖面线，绿线表示坡面上相应点的侵蚀情况，灰色的面积图形在Y轴上的数值为实际上的土壤流失量或土壤沉积量。

图5 原排土场边坡和反S型边坡的土壤流失曲线Fig.5 Soil loss curves of the slope of the original dump and anti-s slope

根据原排土场边坡模型的土壤流失曲线图(图5(a))可知，原排土场边坡的土壤流失从第1坡面开始，而在第2平台上有土壤颗粒的沉积，由此可看出平台设置对边坡保土能力有一定的影响。但在第2坡面位置，土壤损失量比第1坡面位置提升了一倍左右，而最终在第2坡面底部土壤损失量达到最大，最大分离量为116 kg/m2。

根据反S型边坡模型的土壤流失曲线图(图5(b))中可知，边坡的土壤损失量从坡顶开始增长，在到达边坡长度50 m左右的时候在47 kg/m2处上下波动并保持到坡底。最大分离量为48.4 kg/m2。而在边坡长度为51.5 m处的位置土壤损失量有一个骤减，之后继续增加。由于此位置接近反S型边坡凸起部分和凹陷部分的拐点位置，因此推测出现此峰值的原因是此位置剖面曲率趋于平缓，地表径流流速减慢，导致地表径流中携带的土壤颗粒得以沉降，在此处部分堆积使得土壤损失量减少。

4 讨 论

4.1 反S型边坡与其他类型边坡的对比

由于自然地形的复杂性，在胜利矿区未受采矿扰动的自然区域内对边坡剖面线进行提取时，获得的坡型并不只有反S型边坡一种。陡坡部分在提取过程中总共获得60条边坡剖面线(图1)，其中反S型边坡28条，凸型边坡21条，S型边坡6条，凹型边坡4条，直线型边坡1条，反S型边坡和凸型边坡分别占边坡提取总数的46.67%和35%。缓坡部分在提取过程中总共获得55条边坡剖面线，其中反S型边坡24条，凸型边坡15条，S型边坡3条，凹型边坡6条，直线型边坡7条，反S型边坡和凸型边坡分别占边坡提取总数的43.63%和27.28%。由于后3种边坡在本次提取中获得的数量过少，可以认为其并非为胜利矿区自然边坡的主要坡型。而在对凸型边坡进行特征参数提取和边坡模型构建并进行水蚀分析后，发现该边坡模型的土壤损失量高于反S型边坡，因此采用反S型边坡模型对排土场边坡进行设计。

在对陡坡部分反S型边坡进行模型构建时，统计实测数据的凸面水平占比并求取均值为46.82%，根据凸面水平占比和凸面竖直占比的拟合式(4)，计算得出凸面竖直占比为38.01%，即凹面部分的竖直占比为61.99%，十分接近61.80%。这表明陡坡部分的反S型边坡模型，其凸面向凹面过渡的拐点在坡高方向上的投影位置在黄金分割点左右。据此可以推测，自然界的反S型边坡，其凹面部分和凸面部分在竖直方向上的拐点位于黄金分割点附近，蕴含着丰富的美学价值。

4.2 反S型边坡的保土能力和长期稳定性

保土能力是评价边坡模型优劣的重要依据。边坡的保土能力越强，土壤损失量越少，其稳定性就越高。随着坡高的增加，无论是原排土场边坡模型还是反S型边坡模型的年均土壤损失量都有一定程度的提升。原排土场边坡模型在坡高分别为30,45,60 m时的1 a年均土壤损失量分别比15 m时的1 a年均土壤损失量增加了102.67%,201.46%,298.62%，而反S型边坡模型同等条件下的增加百分比分别为87.41%,148.65%,197.88%，相较于原排土场边坡模型分别减少了15.26%,52.82%,100.74%，并且原排土场边坡模型的年均土壤损失量增幅差异不大，而反S型边坡模型的年均土壤损失量增幅呈现递减趋势。这表明随着坡高的增加，反S型边坡模型的保土能力相较于原排土场边坡模型越来越强。

另一方面，随着模拟年份的增加，两种边坡模型的年均土壤损失量都有一定程度的提升，这可能是因为边坡在经历数年的侵蚀导致边坡形状发生改变后，其保土能力有所下降，因此边坡模型的长期稳定性也是一个重要的评价依据。随着模拟年份的增加，边坡模型在不同坡高条件下的年均土壤损失量增加百分比相差不多，如原排土场模型在坡高为15,30,45,60 m的10 a年均土壤损失量相较于1 a年均土壤损失量增加百分比分别为34.02%,34.97%,33.63%,34.08%，故可以用均值代替。原排土场边坡模型各坡高条件下的10 a年均土壤损失量相较于1 a年均土壤损失量提升了34.17%左右，50 a年均土壤损失量相较于1 a年均土壤损失量提升了47.41%左右;而反S型边坡模型分别提升了28.2%和38.44%左右。这表明反S型边坡模型在长时间尺度下的保土能力优于原排土场边坡模型。

4.3 近自然边坡模型的应用和局限

地面坡度越小，地表越平缓，其受侵蚀的影响效果也就越低[23]，从这种角度而言使用缓坡部分的边坡模型貌似是最佳选择。但考虑到实际情况，现有排土场的坡角为17°，若把坡度放缓到缓坡部分(大多<10%)，那将占用坡底部分极多的土地，与实际情况不符。而陡坡部分的反S型边坡模型，其坡角与原排土场的坡角相近，但其保土能力和长期稳定性的提升效果十分显著，故使用陡坡部分的反S型边坡模型对胜利矿区排土场进行设计是可行的。

本文构建边坡模型的方法不仅适用于排土场边坡，还可用于自然地形和设计地形的衔接区域。长久以来，各近自然重塑研究都着眼于以流域为主的整体地形构建[24]，而构建形成的地形往往与自然地形不能很好的进行衔接。本文使用的边坡模型构建方法为解决设计地形与自然地形的衔接问题提供了一种思路，即根据研究区周边自然地形的边坡特征参数，采用统计学方法构建出自然地形和设计地形的过渡区，使得设计地形能较好的与自然地形进行衔接，与当地的流域地貌相协调。

另一方面，本文构建边坡模型的方法也可用于矿区复垦规划中的排土场设计。根据开采前获取的埋深、采深、采厚等地质信息估算外排土场的土方量，计算松散系数和下沉系数后设计外排土场高度，利用本文所构建的近自然边坡模型对外排土场边坡进行优化设计。在对内排土场进行设计时，由于采掘区域的煤层埋深和采深一般具有空间异质性，可以利用本文构建的近自然边坡模型将内排土场重塑为近自然地形而非传统复垦模式中的平地，并对内、外排土场与周围地貌的衔接区域进行优化，将具有一定高度差的衔接区域设计成近自然的反S型边坡，使其能与周围的自然环境相协调。

本文对边坡模型保土能力的评价是基于WEPP土壤水蚀预报模型进行的。该模型的优点是能在短时间内模拟不同时间尺度下的土壤侵蚀过程，相较于野外不同时段实测来进行对比的方法而言，简化了研究过程并降低了研究难度。但应用WEPP模型进行评价的方法存在一定的局限性。一方面，WEPP模型是基于美国本土环境开发的水蚀模型，因此在很多文献中提到过WEPP模型的模拟过程中有部分经验参数，虽然通过气候、土壤等数据能减弱这些经验参数的影响，但还是有可能造成模拟结果上的误差。另一方面，WEPP模型是基于水力侵蚀进行的预报模型，对以水力侵蚀为主的胜利矿区而言适用性较好，但对于风力侵蚀、冻融侵蚀影响效果远大于水力侵蚀的区域，WEPP模型的预测结果可能与实际结果有偏差。

5 结 论

(1)胜利矿区自然区域内边坡类型主要为反S型边坡。通过对边坡形态进行分析，提出了反S型边坡不同于其他类型边坡的两个边坡特征参数，即凸面水平占比和凸面竖直占比。分别对陡坡部分和缓坡部分的反S型边坡的边坡特征参数进行相关性分析发现，陡坡部分边坡特征参数的坡高与坡长、凸面水平占比与凸面竖直占比、坡长与凸面曲率、坡长与凹面曲率显著相关;缓坡部分边坡特征参数的坡高与坡长、凸面曲率与凹面曲率、凸面水平占比与凹面水平占比显著相关。

(2)在对边坡特征参数进行拟合后分别得出陡坡部分和缓坡部分的反S型边坡模型，并发现陡坡部分的反S型边坡模型，其凸面和凹面在竖直方向上的占比符合黄金比例。考虑到实际情况，确定了以陡坡部分反S型边坡模型与原排土场边坡模型进行对比分析。

(3)在WEPP模型中输入胜利矿区的气候、土壤、植被覆盖等数据，对不同坡高、不同模拟年份条件下反S型边坡模型和原排土场边坡模型进行了水蚀分析，得到了两种边坡模型的土壤损失量及土壤流失曲线图。结果表明随着坡高和模拟年份的增加，反S型边坡模型的保土能力和长期稳定性明显优于原排土场边坡模型，验证了反S型边坡模型在排土场边坡设计中的适用性。

(4)本文基于自然边坡特征参数构建近自然边坡模型的方法不仅适用于排土场坡形设计，也为解决矿区土地复垦活动中设计地形与自然地形的衔接问题提供了一种思路，可作为其他地区进行地形重塑的理论依据。