露天矿内排土场近自然地貌重塑研究——以新疆黑山露天矿为例

夏嘉南，李根生，卞正富，雷少刚，3，宫传刚，李 恒

（1.中国矿业大学 矿山生态修复教育部工程研究中心，江苏 徐州 221116；2.中国矿业大学 公共管理学院，江苏 徐州 221116；3. 山东省煤田地质局第二勘探队 山东省采煤塌陷地和采空区治理工程研究中心，山东 济宁 272000）

0 引 言

内排土场作为露天矿区重要的组成部分，其占地面积随采复工作面推进不断增大[1]。在内排过程中，传统地貌重塑方法往往将原始微起伏地貌景观转变为平台和边坡相间的“梯田式”人工堆叠景观，使原始地表要素及水系特征消失[2]。虽然在后续复垦设计中，内排土场内布设输水管网以满足复垦植被生长需求，但整体上未能实现其与周边地貌融合，易引发土地退化等问题[3-4]。因此，需尝试在采复过程中保留原始自然地貌特征以破解上述困境。

近年来，不少专家学者就内排土场规划及内排方案开展了研究。例如，王东等[5]以内排土场空间利用最大化为目标，基于采场工作线及边坡角、内排土场边坡角、基底倾角对排土线的影响规律，构建了采场工作帮、内排土场、端帮和基底的空间几何模型；栗嘉彬[6]以技术最佳适用性为目标，根据物料用途，构建了露天矿表土最优采集周期、最佳覆土厚度等求解模型；OGGERI等[7]以最小复填成本为目标，根据场地、运输等诸因素，提出了多情境下露天矿采复处置方法。然而上述研究中，未充分考虑重塑地貌结果与周边自然地貌起伏贴合情况，导致内排土场往往被重塑为平台加斜坡的“简单梯田式”形态，人工地貌形态明显，易引发严重的景观破碎化问题[8-10]。相较于人工规则地貌，自然地貌作为长期演变的结果，具有高度的区域适宜性与稳定性，可作为修复区地貌重塑的参考对象[3,10-11]。且国外实践证明[12-13]，基于自然的地貌重塑能更好地适应当地水气条件，在有效提高区域抗水蚀能力的同时仅需要较少的维护费用[14]。因此，迫切需要一种贴合周边自然地貌的内排土场重塑方法，用于解决内排土场设计结果与周边自然地貌间贴合的问题[15]。

基于此，笔者以新疆黑山露天矿为例，提出一种采排复一体下基于调整曲面的内排土场近自然地貌重塑模型，并结合地理空间分析技术与MATLAB软件，构建内排土场近自然地貌重塑结果。并通过CLiDE景观演化模型，模拟了研究区10 a间原始地貌、传统设计地貌及近自然设计地貌三者的土壤水蚀情况，对比评价了近自然设计地貌抗水蚀效果，为相同及相似地区内排土场地貌形态设计提供借鉴。

1 研究区概况与数据来源

1.1 研究区概况

研究区位于新疆托克逊县西北约90 km处，地理坐标：87°26′04″～87°29′40″E，43°13′48″～43°14′01″N，周边地貌以山间谷地为主；年均降水量约160.3 mm，年均气温约1.6 ℃，冬寒夏凉，属大陆干旱及高寒气候。区内表土以第四系为主，野生植被类型为典型草原植被。

依据地质勘探资料，选取吉布拉克煤矿探矿权范围以南为研究区（图1）。区域东西长约4 874 m、南北长约1 000 m。其中，预设采坑边界东西长约812 m、南北长约1 000 m，采矿作业由东至西；区域可采煤层分别为6、7、8、9、11、12-1、12-2、13-1和13-2号煤层，其中13-2号煤层为全区可露采、其余为局部可露采。

图1 研究区示意Fig.1 Schematic of the study area

1.2 数据来源

研究区采前自然原始数字高程数据（DEM）来源于地理空间数据云，拍摄时间2011年12月，条带号87，行编号43，空间分辨率30 m×30 m；研究区煤层底板、顶板标高及水平坐标数据来源于《新疆托克逊县黑山矿区黑山露天煤矿资源储量核实报告》。此外，基于上述报告通过ArcGIS克里金插值模块获取煤层底板及顶板空间分布数据；内排土场采、复子区空间分布数据基于内排土场采复子区位置识别方法获取；复填子区内的可用土方量数据通过可用土方量计算式求解获取；复填区地表调整曲面数据通过地表调整曲面预构建模型获取；复填子区地表近自然DEM通过MATLAB软件基于动态土方量一致、且区域斜率较缓的原则定向筛选不同B样条控制点位置下的地表调整曲面获得；研究区自然原始DEM数据通过ArcGIS空间叠加提取获得；研究区内排土场传统设计DEM数据基于土方量动态平衡下“30°斜坡加平台”的设计思路获取；最后运用CLiDE演化模型与MATLAB软件，分别模拟并统计研究区原始DEM、传统设计DEM及近自然设计DEM三者10 a的土壤水蚀量，以评估近自然整形结果的抗侵蚀能力。

2 研究方法

2.1 内排土场近自然地貌重塑模型构建

2.1.1 采复子区空间位置识别

依据露天矿采复周期，获取图2所示各周期下开采子区amn和复填子区afn水平投影空间位置。

图2 露天矿采复周期下开采及复填子区空间位置提取示意Fig.2 Schematic of extraction of spatial location of the mining sub-area and the refilling sub-area under the open-pit mining-refill cycle

如图2所示，依据露天矿实际采复周期，沿开采方向将露天矿区按采复周期时间先后顺序划分为n个子区。其中，开采子区和复填子区分别对应区域amn和区域afn。在露天矿内排构建的采复过程中，采坑每前进一个开采子区amn对应修复后侧一个复填子区afn，直至内排过程结束，矿坑由区域Ab移至区域Aa。

2.1.2 可用土方量计算

为便于后续土方量计算，将开采及复填子区均细分为rowc行colc列的网格单元组（图2）。在采复周期n中，开采子区amn可用于对应复填子区afn的土方量，计算式如下：

（1）

式中:V（amn）为开采子区amn可用于对应复填子区afn的土方量，其值是开采子区afn最下层可采矿层底板以上的总物质体积减去区域内所有可采矿层总物质体积；i,j分别为开采及复填子区内网格单元的行、列序号，取取值分别在1～rowc和1～colc之间；H（amnij）为开采子区amn网格单元（i,j）所在区域的采前地面标高；hd（amnij）为开采子区amn网格单元（i,j）所在区域的最底层矿层底板标高；ht（amnij）为开采子区amn网格单元（i,j）所在区域的可采矿层总厚度；S（amnij）为开采子区amn子单元（i,j）的水平投影面积，其余参数解释同上。

2.1.3 复填区地表调整曲面预构建

依据区域原有自然稳定地貌，构建图3所示调整曲面模型，以尽可能保留复填子区afn原有地貌起伏特征。

图3 地表调整曲面预构建示意Fig.3 Schematic of pre-built adjustment surface

如图3所示，将反映复填子区最终地貌设计形态较其原始自然地貌空间高程变化的曲面定义为地表调整曲面。为保证地表调整曲面能使复填子区表面与周边自然地貌“无缝”融合，选取边界处高程变化小，中间变化大，且与区域自然地貌相似的连续放缓曲面。参照“船体”构建思路，采用连续性较优的B样条曲线[16]为“龙骨”，三角函数为侧方“支架”，构建复填子区afn所对应的地表调整曲面。图中虚线箭头为开采方向；x,y,z三轴正方向分别对应开采方向、水平垂直开采方向及竖直方向；fb和ft分别对应B样条曲线函数和三角函数。其中，B样条曲线fb位于复填子区afn水平投影中心线位置，采用控制点构建，其计算式如下：

fb（j）=Bezier（v,tb,ts,te）

（2）

（3）

其中:fb（j）为复填子区afn第j列的B样条曲线值，定义域在0～1；Bezier为竖直平面上的2维B样条曲线函数，形态由起始点0,y1、终止点1,yNC和中间NC-2个控制点共同控制；NC为Bezier函数中起点和终点在内所有控制点的个数，取值为>2的整数；v为NC行2列的矩阵，由上至下行数据分别对应起点、控制点1至控制点NC-2、终点的坐标；tb为起点所对应网格单元的列序号；ts为Bezier函数取值间隔，是终点与起点所对应网格单元的列序号之差，取值为大于0的正数；te是终点所对应网格单元的列序号；x1～xNC-2分别为控制点1至控制点NC-2在中心线上以起点为原点的水平投影距离，其中x1≤x2≤…≤xNC-2；y1～yNC分别是起点、控制点1至控制点NC-2、终点的预设高程值，定义域皆在0～1之间，其余参数解释同上。在确定复填子区B样条曲线函数形态的基础上，通过三角函数确定地表调整曲面，其计算式如下：

（4）

式中:ft（i,j）为复填子区afn地表调整曲面网格单元（i,j）所对应的预设高程值，取值在0～1；π为圆周率，其余参数解释同上。

2.1.4 曲面土方量控制及坡度缓和优化

为保证采复周期中可用土方量动态平衡，通过竖直伸缩变化使地表调整曲面符合土方平衡要求，其计算式如下：

（5）

（6）

（7）

在此基础上，为筛选出区域整体坡度较缓的复填子区地表高程模型，以提高其表土抗水蚀能力，构建以网格单元与周边网格单元间坡角为对象的区域坡度评分标准，其计算式如下：

（8）

其中：fp（a）为坡度夹角a所对应的坡度评分值，其值为大于0的整数；au为正整数，其取值在2°～90°，其余参数解释同上。在区域坡度评分标准构建的基础上，构建区域坡度评分模型，其计算式如下：

（9）

式中：p（afn）为复填子区afn的区域坡度总评分；fpij（a）为子区afn内网格单元（i,j）沿某一朝向的坡度评分值；∑fpij（a）为子单元（i,j）与周边8个网格单元间坡角所对应评分值的总和。

依据上述评分模型，通过MATLAB软件，调整2.1.3节中B样条曲线起点、控制点1至控制点NC-2、终点的坐标位置，遍历构建一组形态各异的复填子区地表高程模型。并以区域坡度总评分最小为目标，定向筛选复填子区Ht′求解结果，最终获取原有自然形态表面下，土方量动态平衡且坡度平缓的复填子区afn地表高程近自然设计形态。

2.2 设计结果抗水蚀能力评价

CLiDE模型是基于CAESAR模型开发的景观地貌演化模型[17]。与传统的侵蚀预测模型相比，其不仅能获取侵蚀和沉积的时序数据，还能获取最终发育的地貌形态。同时嵌入了Lisflood-FP流体力学相关模块与地下水分布式水文模块，可以模拟地表径流-地貌-地下水三者间的交互作用和补给关系[18]，故极其适用于地貌修复效果评价[10]。为评价内排土场设计结果的抗水蚀能力，将研究区内排土场近自然设计结果为试验组，原始自然地貌及传统复填方法下内排土场设计结果为对照组，结合实地调查情况构建CLiDE模型，模拟其10 a的土壤侵蚀过程。并基于MATLAB软件统计各组地貌演变前后土方变化量，以评价其抗水蚀效果。

3 结果与分析

3.1 数据准备结果及分析

研究区可采煤层共计9层，依据储量报告，采用ArcGIS克里金插值模块获取研究区13-2号煤层底板标高平面分布图及其各煤层厚度平面分布图（图4）。为避免数据空间分辨率差异造成的干扰，将影像空间分辨率统一设置为30 m×30 m，数据形式均为162列32行的栅格影像。

图4 13-2号煤层底板标高平面分布与研究区各煤层厚度平面分布Fig.4 Plane distribution map of floor elevation of No. 13-2 coal seam and plane distribution map of the thickness of each coal seam

3.2 近自然设计地貌构建结果及分析

依据2.1.1节采复子区空间位置识别方法，在30 m×30 m栅格影像数据下，按采复周期将研究区划分为45个一一对应的采复子区。其中，采复子区网格行数为32，列数为3，即复填子区和开采子区均设计为水平投影东西长90 m，南北长960 m的水平矩形区域，以使采复工作前后采坑的水平投影大小不变（图5）。

图5 研究区开采子区及复填子区空间位置示意Fig.5 Schematic of the spatial location of the mining sub-area and the refill sub-area

在此基础上，依据2.1节所示内排土场地表近自然设计方法，通过MATLAB软件以0.01为控制点数值最小变化间距，膨胀系数k取均值1.1，求解研究区内排土场地貌近自然设计DEM结果（图6）。其中，考虑到矿区可采煤层非惟一（9层），且厚度及底板标高存在差异，依据研究区三维地质资料，以13-2号煤层的底板标高为研究区煤层底板标高，并采用空间叠加方法叠加各煤层厚度数据，获取研究区煤层顶板标高等效数据，以实现本研究内排土场近自然地貌重塑模型于多矿层复杂情景下的应用。

图6 研究区内排土场DEM图Fig.6 DEM of the internal dump in the study area

如图6所示，视觉效果上，研究区内排土场近自然设计DEM与原始自然地貌及纹理特征较为相似，且高程值均在+2 000.26～+2 959.00 m；表土高程上，受土方平衡约束，近自然设计DEM较自然原始DEM存在差异，变化幅度在-205.90～50.03 m。一方面，受研究区“开采子区地表原始标高”、“复填子区煤层底板标高”和“复填子区煤层顶板等效标高”的直接影响，近自然设计结果较原始地貌发生不规则形变。例如，在东西方向上，复填子区af1和af14近自然设计DEM较原始地貌有所抬升，而在其他区域较原始地貌有所下降。另一方面，受近自然地貌重塑模型中地表调整曲面模型影响（2.1.3节，地表调整曲面突出部分集中在区域中部），使得近自然设计DEM较原始地貌间形变也集中在区域中部。

与前者相似，相较于区域自然原始DEM，传统设计结果也具有中部形变明显的空间立体特征（图6c、图7c），然而与前者不同，受斜坡设计影响，内排土场传统设计DEM与周边原始地貌间地貌“分割”明显。例如，在同一位置由北至南分别做原始DEM、近自然设计DEM及传统设计DEM三者的竖直剖面，在相同可用土方量控制下，相较于前两者地貌“波动式”形态（图7a、图7b），传统设计DEM具有明显的“斜坡-平台-斜坡”三段式的地貌形态特征（图7c）。

在土方运移上，依据采复子区土方对应关系，运用MATLAB计算获取近自然设计及传统设计下土方平均运距。其中，近自然设计DEM与传统设计DEM两者土方平均运距分别为822.20 m/m3与821.71 m/m3，前者较后者仅高出约0.06%，两者土方平均运距极为接近（相差0.26 m/m3）。

3.3 设计结果抗水蚀能力评价

地貌稳定性主要受气候、地质、植被等诸要素影响。其组成物料强度直接影响边坡的稳定性，是其抗侵蚀能力的主要指标。依据地质调查报告，确定研究区内排土场地层结构及岩土物理力学参数，并结合中国天气网获取研究区2007—2017年气候参数。通过CLiDE演化模型，分别模拟3.2节所示10 a间内排土场自然原始DEM、传统设计DEM及近自然设计DEM的演变过程。其中，为保证模拟结果具有可比性，将研究区内排土场传统设计DEM的采复子区空间位置与3.2节近自然设计DEM设定为一致，其最终设计结果如图7c所示。同时，为避免模型边缘效应影响，将研究区周边1 000 m范围共同纳入运算（图7）。最后通过MATLAB叠加及对比分析，统计研究区内排土场自然原始DEM、传统设计DEM及近自然设计DEM下10 a的土壤水蚀总量（表1）。

如表1所示，整体上，本研究内排土场10 a土壤水蚀量在4.098 3×105～23.681 5×105m3。其中，自然原始DEM土壤抗水蚀能力最强，模拟情境下其土壤10 a水蚀总量仅占样本平均水平（14.069 5×105m3）的29.13%。在人工地貌设计下，本研究近自然设计DEM相较于传统设计DEM可减少约39.07%的土壤水蚀量，其土壤水蚀总量逼近样本平均水平。

为进一步分析排土场各种地貌设计方式下对区域整体土壤水蚀的影响，将各DEM演化模拟结果与其初始DEM求差，获取研究区内排土场表土高度的空间变化影像（图8），其中正数为该区域表土高度增高高度，负数为降低高度。

图7 研究区内排土场自然原始、传统设计及近自然设计下DEM影像及横向剖面示意Fig.7 DEM image and cross-sectional schematic of the natural primitive, traditional design and near-natural design of the internal dump

表1 研究区内排土场自然原始、传统设计及近自然设计下区域土壤10 a总水侵蚀量

图8 研究区内排土场演化前后自然原始、传统设计及近自然设计下DEM表层高度变化Fig.8 The surface height changes of the natural original DEM, the traditional design DEM, and the near-natural design DEM before and after the evolution of the dump

如图8所示，整体上，各样本土壤在10 a的演变模拟中均发生了不同程度的运移，使区域原有DEM发生不同程度的改变。在表土高度变化上，各DEM样本高程变化区间均在-0.62～2.79 m之间，其变化幅度由大到小依次为：传统设计DEM（3.39 m）、近自然设计DEM（3.35 m）和自然原始DEM（3.06 m）。排序结果与各样本土壤10 a水蚀总量排序一致。在空间分布上，各样本的土壤堆积及流失区域均有所不同。其中，传统设计DEM表土流失及堆积区与研究区边界高度相关，且主要沿边界（图8研究区边界）呈向内规则分布，视觉上与周边自然地貌间具有明显割裂，景观及水文融合性较低。与前者不同，近自然设计DEM表土堆积及流失区空间分布与自然原始DEM较为相似，两者均呈现出高度复杂的“离散化”空间分布。在周边自然地貌融合上，两者与周边自然地貌间无明显“痕迹”间隔，表土流失及堆积区域空间分布与研究区边界无明显相关，景观及水文衔接性较好。并且，以自然原始DEM表土流失区（演化后高程下降区域）和堆积区（演化后高程上升区域）空间位置为参照的空间叠加分析表明，相较传统设计DEM（空间重叠率：54.31%），近自然设计DEM土壤流失及堆积区具有更高的空间重叠率（57.27%）。

4 讨 论

内排土场作为露天矿山的重要组成部分，其占地面积随工作面的推进不断增大。然而，受传统复填方式影响，内排土场区域原始微起伏地貌往往被替换为斜坡加平台的“梯田式”人工堆垫地貌，使原有地表要素及水文特征丧失，并在集中降雨条件下，易引发表土水蚀等问题[2-4]。为此，专家学者提出依据临近自然稳定地貌特征以实现排土场地貌重塑[3,8,10]。例如，有学者以胜利一号矿排土场边坡为例，提出以仿自然“反S形”边坡替代传统斜坡，减少区域表土侵蚀[10]；同样，有学者以胜利一号矿为例，提出依据内排土场周边自然地貌参数，基于区域整体土方平衡，构建内排土场近自然地貌结果，实现景观内外融合，且相较传统设计地貌具备更好的表土稳定性[3]。与上述研究相似，本研究同样以自然地貌为学习对象，且内排土场近自然设计地貌较传统设计地貌具有更强的抗水蚀能力（10年土壤水蚀减少约39.07%）。然而不同的是，上述研究强调以区域地貌统计特征为依据，本研究则强调以原始自然稳定地貌为依据，即在排土场塑形过程中尽可能保留原始自然稳定地貌特征，以实现矿区内外地貌相连。

在露天矿采复过程中，受土体膨胀系数、可采煤层标高、采复周期及采前地表标高等诸要素影响，使得各采复周期间后方复填可用土方量存在差异；且由于复填区形态各异，造成内排土场复填地貌存在走高、走低等诸多情况。在近自然设计中，已有坡面近自然重塑方法虽然能够较好地应用于传统内排土场复填区地貌设计，然而其实际形态仍以去除边角的“梯田式”地貌为主[8,10]。后者提出的内排土场近自然塑形方法虽具有良好的景观融合效果，并考虑到开采前后土方平衡，然而设计中未考虑采复周期影响[3]，易造成采复周期中复填区实际可用土方量无法或超出原有设计表面等问题。受采复周期影响，本研究近自然设计地貌也呈现出部分区域较自然原始表面抬高、沉降等高程形变特征。然而与已有方法不同的是，本研究内排土场近自然地貌重塑方法通过采复子区一一对应的方式解决采复过程中土方动态变化的问题，并相较于传统设计地貌，近自然设计地貌与周边自然地貌无明显分界，模拟构建上可实现采后地形与周边自然地形的融合。

值得注意的是，在自然地貌中，区域原始地貌并不总是稳定的[19]。在一定条件下，传统复垦后矿区或将更有利于区域生态恢复[20]。由于本研究近自然地貌重塑方法以采前原始地貌为学习对象，故在构建内排土场地貌形态前，需充分论证原有自然地貌是否稳定（例如，本研究自然原始地貌10年水蚀量最小，仅为传统设计地貌的28.40%），以选择有利于区域生态恢复的内排土场设计方法。同时，受地下煤层剥离影响，本研究设计结果无法在形态上与原始自然地貌保持一致，并受采后松散土体空间不规则压实沉降[21]，设计结果可能远离预期。故在后续研究中，可进一步从力学视角，优化本研究基于调整曲面的内排土场近自然地貌重塑模型。

5 结 论

1）研究区内排土场近自然设计结果表明，提出的近自然地貌重塑方法可在土方动态平衡的前提下，依据采复子区采前地表、矿层顶板及底板数字高程数据，获取可用地表标高设计结果；在土方调配上，近自然设计地貌土方平均运距为822.20 m/m3，较传统设计地貌高0.06%，两者土方平均运距极为接近；且在视觉效果上，与周边自然景观“无痕”拼接，对区域开采前后地貌改变较少。

2）研究区10 a土壤水蚀模拟结果表明，本研究内排土场近自然整形结果较传统“梯田式”减少约39.07%的土壤流失量，表土抗水蚀能力较优；在土壤演化运移中，本研究近自然设计地貌与区域原始地貌在土壤堆积及流失上空间重叠率达57.27%，较传统设计结果高2.96%，且视觉上呈“离散化”分布，与周边自然地貌衔接较好。