水文保留曲面下草原露天矿内排土场地貌近自然重塑

2022-08-18夏嘉南李根生卞正富雷少刚金慧霞

煤炭学报 2022年7期

夏嘉南，李根生，卞正富，雷少刚，金慧霞,3

(1.中国矿业大学公共管理学院，江苏徐州 221116；2.中国矿业大学矿山生态修复教育部工程研究中心，江苏徐州 221116；3.温州职业技术学院，浙江温州 325035)

露天开采作为我国草原矿区的重要开采方式，其内排土场占地面积随开采工作面的不断推进而增大。在此过程中，传统地貌重塑方法往往将原始地貌转变为“台阶式”或“锥状堆积式”的人工堆叠地貌，使原有的地貌及水文特征消失，导致矿区重塑地貌与周边自然地貌割裂。虽然在后续复垦设计中，内排土场通过输水管网布设以满足复垦区重建植被的生长需求，但未能实现矿区内外地貌水系的“无痕”拼接，从而在高强度降雨条件下，易引发边界处严重表土水蚀等问题。相对人工规则地貌，自然稳定地貌作为长期演变的结果，具有较好的适宜性与系统性。因此，可尝试在采排复过程中保留自然地貌及其水文特征以破解上述困境。

当前，GeoFluv模型作为主要的近自然地貌重塑方法，被广泛应用于采损区地貌重塑，且塑形效果在美国墨西哥州和怀俄明州、西班牙、澳大利亚等矿区均得到了良好的检验，并在国内得到了初步的运用。然而，受模型框架影响，重塑区水文沟道位置需人为主观设定，虽后续通过自然稳定沟道特征值加以修正，但整体上其地貌求解结果主观性较高。NbS(Nature-based Solutions)理论谋求人与生态系统互惠，与“保护优先，自然恢复为主”的中国生态文明建设理念相契合。在上述研究及理论的指导下，团队内诸多学者均提出师法自然稳定地貌的内排土场地貌重塑方法。其中，李恒等以内蒙古胜利矿为例，通过将周边自然稳定地貌剖面形态应用至排土场边坡，通过GeoWep模拟论证仿自然形态边坡较传统“台阶式”边坡更具抗水蚀能力。陈航等以内蒙古胜利矿为例，通过提取周边自然稳定地貌特征参数，在不考虑采复周期的情景下采用GeoFlue软件构建研究区内排土场近自然重塑地貌，并基于CLiDE模拟验证其表土抗水蚀能力优于传统内排土场地貌重塑结果。然而，在上述研究中，内排土场采排复一体下复填区的可用土方量会随采复周期变化。同时，相较于周边自然稳定地貌统计结果，区域原始自然稳定地貌的特殊性(受区域气候及相邻区的长期作用影响)也未得到充分考虑。在前置研究中，笔者以新疆黑山露天矿为例，提出基于SAS(Surface Adjustment Surface)的内排土场地貌近自然重塑模型，强调通过局部调整的方式以尽可能保留原始地貌起伏特征，同时避免模型构建中水文沟道的人为主观设计，且CLiDE水蚀模拟表明，NNDL_SAS(Near-natural Design Landform with Surface Adjustment Surface)相较于TDL(Traditional Design Landform)在10 a当地气候条件下可减少约39.07%的表土水蚀量。然而，后续论证表明，上述模型求解形态局部调整可能会造成水系断流、淤积等不可控问题。因此，需要一种既贴合周边自然地貌，又保留原始地貌及水文特征的内排土场地貌重塑方法，以削弱前置研究中地貌设计结果水系断流、淤积等问题。

为此，笔者在前置研究的基础上，以新疆黑山露天矿为例，提出一种基于HPCS(Hydrological Preserved Curved Surface)的内排土场近自然地貌重塑模型，求解并以NNDL_HPCS(Near-natural Design Landform with Hydrological Preserved Curved Surface)为实验组，NNDL_SAS,TDL和PNL(Pristine Natural Landform)为对照组，结合GIS空间分析软件与CLiDE演化模型，模拟并讨论NNDL_HPCS在地貌融合、土方调配和土壤抗水蚀3方面的改良效果。为相同及相似地区NbS理论下矿山近自然生态修复提供方法支持。

1 研究区概况与数据来源

1.1 研究区概况

研究区位于新疆黑山露天矿矿权范围内，地理坐标：87.43°E—87.49°E，43.23°N—43.23°N，周边地形以山间谷地为主，属大陆半干旱及高寒气候；区内表土以第四系为主，野生植被类型为典型草原植被。采区东西长约4 874 m、南北长约1 000 m；采坑预设边界为812 m×1 000 m的矩形，内排作业由东至西(图1)。

图1 研究区采前原始DEM和真彩色卫星影像Fig.1 DEM of the natural surface of the study area before acquisition and RGB satellite image of the study area

依据研究区煤炭资源储量核实报告，区域可采煤层为6,7，8,9，11,12-1,12-2,13-1和13-2号。其中13-2号煤层为全区可露采、其余为局部可露采(图2)。内排土场台阶高度均设置为10 m。其中西帮为开采帮，帮坡角为20°；南北两侧边帮帮坡角均为13°；东帮为内排帮，帮坡角为8°。

1.2 数据来源

研究区采前原始DEM(Digital Elevation Model)、煤层顶板DEM和煤层底板DEM均来源于前置研究，空间分辨率均为30 m×30 m；内排土场开采和复填子区空间分布数据与前置研究相同，均基于内排土场采复子区位置识别方法获取；复填子区内的可用土方量数据通过可用土方量计算式求解获取；复填区HPCS数据通过HPCS预构建模型获取；复填子区NNDL_HPCS通过Matlab软件基于土方平衡下区域坡度最缓原则定向筛选不同控制点系数下预构建的HPCS获得；复填子区NNDL_SAS来源于前置研究成果中“近自然设计DEM”；研究区PNL数据来源于前置研究成果中“自然原始地貌”；研究区内排土TDL与前置研究成果中的“传统设计DEM”一致(30°斜坡加平台式“阶梯状”布局)；地貌融合效果中，各DEM相对于PNL的高程方差值通过ArcGIS空间叠加分析方法，逐栅格累计统计DEM与PNL间栅格高程差的平方值获取；地貌融合效果中，自然水文保留率指标通过ArcGIS坡向模块和空间叠加分析方法，提取并对比研究区NNDL_HPCS,NNDL_SAS和TDL三者图像栅格坡向与PNL空间分布一致性获取；土方调配中，NNDL_HPCS,NNDL_SAS和TDL的土方平均运距均由Matlab基于每立方土体空间直线位移距离统计获取；最后运用CLiDE演化模型和Matlab软件，分别模拟并统计研究区NNDL_HPCS,NNDL_SAS,TDL和PNL四者10 a的土壤水蚀量，以评估NNDL_HPCS的抗水蚀能力。

图2 研究区各煤层厚度及最底层煤层底板标高空间分布Fig.2 Spatial distribution of the coal seam thickness and bottom plate elevation of the lowest coal seam in the study area

2 研究方法

2.1 基于HPCS的内排土场近自然地貌重塑模型

2.1.1 采复子区空间位置识别

依据露天矿采复周期，获取开采子区m和复填子区f的水平投影空间位置(图3)。

图3 矿区采复周期n中采前和复填后竖直剖面示意((a)，(b))以及采复子区网格单元示意(c)Fig.3 Schematic diagram of the vertical section before and after mining in the mining area under the open-pit mining-refill cycle n， schematic diagram of the grid unit of the mining-refill sub-area

如图3所示，依据露天矿实际采复周期及其时间先后顺序，沿开采方向将露天矿区划分为个子区。其中，开采子区和复填子区分别对应区域～m和区域～f。在露天矿内排构建的采复过程中，采坑每前进一个开采子区m对应修复后侧一个复填子区f(图3(a)，(b))，直至结束内排过程，矿坑由区域不断迭代前进至采矿终止边界。

2.1.2 可用土方量计算

基于输入数据栅格空间分辨率，将开采及复填子区均细分为row行col列的网格单元组(图3(c))。其中，网格单元大小与栅格分辨率一致。在采复周期中，开采子区m可用于复填子区f的土方量计算式为

(1)

式中，(m)为开采子区m可用于复填子区f的土方量，其值是开采子区m总岩土体积(最底层底板以上所有物质)与其可采矿层总物质体积之差；为土体修正系数，是煤层受剥离损失、夹矸、岩土膨胀等多因素综合影响下实际土方量与剥离土方量之比；，分别为采复子区内网格单元的行、列序号，定义域分别为1～row与1～col；(m)为开采子区网格单元m(,)所在区域的采前地面平均标高；(m)为开采子区网格单元m(,)所在区域最底层矿层的底板平均标高；(m)为开采子区网格单元m(,)所在区域的可采矿层平均厚度；(m)为开采子区网格单元m(,)的水平投影面积，是行列分辨率的乘积，其余参数解释同上。

2.1.3 复填区HPCS预构建

自然稳定地貌作为区域长时序演化的结果，其不仅形态稳定，且地表水文特征与上下流域间高度匹配。为构建与周边自然地貌相衔接的重塑地貌，依据NbS理论，应对其水文特征进行保留，以实现重塑地貌与周边自然地貌相融合。HPCS模型是在上述理论基础上，通过对水文沟道和坡顶进行识别，令其中水文沟道和边界处高程不变，以保证地貌重塑前后内外水系相通，汇水面关系不变；通过调整区域内各坡顶高程，以实现土方量平衡；并通过目标筛选，以获取最优地貌设计形态的地表高程调整模型。在数值上，其是重塑地貌DEM较采前自然地貌DEM高程变化的空间分布数据。

在HPCS求解过程中，首先需构建一批形态各异的HPCS预构建模型，以为后续土方控制和坡度缓和优化提供可供操作的基底数据。故基于其模型定义与NbS理论，为尽可能保留复填子区f原始水文特征，通过ArcGIS水文分析模块与Matlab软件，创建图4所示HPCS预构建模型。

图4 HPCS预构建示意Fig.4 Schematic diagram of pre-construction of HPCS

如图4所示，对原有自然稳定地貌进行概化，以减少NNDL_HPCS求解结果与PNL间的汇水面分布差异。其中，依据概化要素于原有地貌水文所起作用，划分为分割汇水面的地形突出点(坡顶点)和汇水面内水流集中的沟道线(水系或者沟道)。其中，水文沟道网格提取结果通过ArcGIS水文分析模块获取，其最小汇水面积参数通过均值变点法(OTSU)确定。原始地貌的地形突出点则基于式(2)提取：

(2)

其中，(,)为复填子区网格单元(，)的取值，凸点取值为1，其余为0；area为以网格单元(，)为中点的*高程搜索窗口，area(2，2)为窗口中心网格单元(，)所对应的高程值，其余参数解释同上。当搜索网格单元中心点单元area(2，2)高程为其窗口内最大值，定义该中心网格为凸点。通过Matlab遍历填挖后复填子区PNL，识别凸点网格单元所在位置。在此基础上，将复填子区凸点网格识别结果与水文沟道网格提取结果进行空间叠加，构建式(3)判别式对HPCS网格单元进行预赋值：

(3)

其中，为HPCS网格单元预赋值结果;(,)为其网格单元(,)的取值，凸点网格赋值为w，其余为0；w为凸点网格(，)所对应的无量纲系数，定义域为[0,1]，其余参数解释同上。在此基础上，基于式(4)所示模糊算法，预构建HPCS：

(4)

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

同SAS中的土方动态平衡方法，通过竖直伸缩变化使预构建的HPCS符合土方平衡要求，其计算式为

(5)

′(,)=(,)

(6)

′(f)=′(,)+(f)

(7)

式中，为伸缩变化系数，其值是复填子区f剥离的总岩土(包含矿层)体积减去由开采子区m获取的膨胀后的可用土方量之后，与预构建的HPCS总体积比值的相反数；′(,)为预构建HPCS伸缩变换后网格单元(,)的值，是复填子区f最终设计结果(NNDL_HPCS)与PNL间高程值的变化；′(f)为NNDL_HPCS中复填子区网格单元f(，)所对应的地表标高，其值等于该处的自然原始标高与伸缩变换后的HPCS值之和。

在此基础上，为筛选出区域表土抗水蚀能力较强(坡度平缓)且与周边地貌融合的NNDL_HPCS，构建相邻网格单元间坡角为对象的区域坡度评分标准。在坡度评分标准设置中，为避免区域不同坡度组合间的分值混淆，同时减少软件运算中非必要的内存占用，以1°为步长，令相邻步长间高坡度分值是低坡度分值乘上相邻网格数(即8个)的值再加1。其计算式为

(8)

其中，()为坡度夹角(°)所对应的坡度评分值，其值域∈；为正整数，定义域在2°～90°。在区域坡度评分标准构建的基础上，构建区域坡度评分模型，其计算式为

(9)

其中，(f)为复填子区f的区域坡度总评分值；p()为子区f内网格单元(,)沿某一朝向的坡度评分值;∑p为子单元(,)与周边8个网格单元间坡角所对应评分值的总和。为确保最终筛选对象与周边已知地貌无缝拼接，除复填子区f与矿坑邻接处网格单元外，其余网格单元都需统计与周边8个网格单元的坡度评分值。

将设为w的最小变化步长，通过上述HPCS预构建模型，由0至1遍历获取一批不同w组合下的预构建HPCS，进而获取一组形态各异的复填子区′。在此基础上，以区域坡度评分最小为目标，定向筛选复填子区′求解结果，最终获取原有自然形态表面下，土方量动态平衡且坡度平缓的复填子区f内NNDL_HPCS。

2.2 水文保留率

坡向是地貌的重要表征参数，对地表径流起着直接导向作用，同时受区域光照、温度、雨量、风速等因子的综合影响，间接决定区域植被类型及其生长状况。基于灰色系统理论，以自然稳定地貌坡向空间特征为塑形导向。将水文保留率定义为重塑地貌与采前自然地貌坡向空间分布上一致性，即2者相同坡向区域的面积占总区域面积的比例。通过ArcGIS坡向组件，提取PNL,NNDL_HPCS,NNDL_SAS和TDL的坡向空间栅格数据。以PNL坡向空间栅格数据为基准，通过Matlab软件逐栅格进行空间叠加分析，计算NNDL_HPCS,NNDL_SAS,TDL三者与PNL的区域坡向一致性，其计算式为

(10)

其中，为地貌相对于PNL的水文保留率，其定义域在[0,1]，若2者坡向空间一致性越高(即相同位置2者间的坡向夹角越小)，其取值越大；row，col分别为地貌坡向栅格数据的总行、总列数；为地貌的坡向空间栅格数据，其中(,)为其第行第列的栅格的坡向值，其取值大小符合罗马盘规则；为PNL的坡向空间栅格数据，其中(,)是其第行第列栅格的坡向值，其取值大小也符合罗马盘规则。在ArcGIS坡向模块中，平台坡向赋值为-1，考虑到平台坡向的万向性(可视作任意方向)，默认将平台处的坡向认定为一致。

2.3 土方直线运距

基于露天采排复一体化理念，令开采工作面土方由上层台阶逐层向下剥离，后侧复填工作面由煤层剥离后底板处逐层至上复填。其中，复填区当且仅当相邻低处台阶的可用土方量(考虑膨胀系数)完全满足复填所需时，才进行相邻高处台阶的复填作业。由于实际采排复一体化工作中，土方运移受运输通道等诸因素影响，其难以精确计算，故依据上述思路，按照研究区边帮设计参数，以10 m为分段间隔，通过Matlab建立采复前后土方的空间对应索引，并依据其每m土体中心点的位移，计算其土方平均运距，以此简要评价NNDL_HPCS,NNDL_SAS和TDL的土方调配特征。其计算式为

(11)

其中，为土方平均运距；为运移土方按时间的排序，取值为以1为步长从0至(m)的递进值；(m,)，′(m,)分别为开采子区m中第运移的立方土体运移前后中心点所对应的三维坐标，其余参数解释同上。通过统计各采复周期中各立方土体运移前后位移所对应的模长均值，简要评价重塑地貌土方运移特征。

2.4 抗水蚀能力评价

CLiDE模型是基于CAESAR模型开发的景观地貌演化模型。相比于传统的侵蚀预测模型，在获取侵蚀和沉积的时序数据的同时，还能获取最终发育的地貌形态。并可模拟地表径流-地貌-地下水3者间的交互作用和补给关系，极其适用于地貌修复效果评价。为评价内排土场设计结果的抗水蚀能力，将研究区NNDL_HPCS设为实验组，NNDL_SAS,TDL和PNL设为对照组，套用前置研究设定参数(气象数据来源于中国气象网)构建CLiDE模型，模拟4者10 a的土壤侵蚀过程。其中，依据研究区典型草原植被类型，将土地利用类型统一设置为第5类(适用于半野生草地、粗糙的草地、中性草地等)，并基于区域植被以1 a生为主，将其植被恢复时间统一设置为1 a。并基于Matlab软件统计各组地貌演变前后土方变化量，以评价其抗水蚀效果。其中为排除矿坑最终位置对土壤水蚀模拟结果的干扰，在PNL对应位置处设置相同大小的矿坑。其计算式为

第一，专业人才匮乏。当下养老服务存在的主要困境是养老专业人才队伍的匮乏，导致服务专业水平不足，影响服务质量和服务成效。在养老服务提供方的人才队伍中，主要包括志愿者、居委会工作人员，以及非专业出身的其他领域转行而从事养老服务的人员等，在养老服务的能力和水平等方面都有待考量，这对于前期智慧养老服务模式的建设是一个亟待突破的限制性因素。更值得注意的是，智慧养老的运营模式不仅需要的是专业的养老服务人才，更是需要懂大数据技术、大数据工具的数据专业领域的养老服务复合型人才。

(12)

式中，为土壤水蚀总量；(f)，(f)分别为复填子区网格单元f(,)中地貌演化前后的地表标高均值；(f)为复填子区网格单元f(,)的水平投影面积，为行列分辨率的乘积。

3 结果与分析

3.1 数据准备结果及分析

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

考虑到研究区可采煤层非惟一(9层)，且其厚度与底板标高存在差异，故依据三维地质资料，以13-2号煤层的底板标高为研究区煤层底板标高，并采用空间叠加方法叠加各煤层厚度数据，以获取研究区煤层顶板标高等效数据，实现本研究地貌重塑模型于多矿层复杂情景下的应用。

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

依据2.1.1节采复子区空间位置识别方法，将研究区划分为45个一一对应的采复子区，与已有研究保持一致。其中，采复子区网格行(row)列(col)数均为32行3列，其水平投影为东西长90 m，南北长960 m的矩形。通过2.1节内排土场近自然地貌重塑模型，将w的最小变化步长取0.01，土体修正系数取1.1，求解研究区NNDL_HPCS(图5(b))。

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

3.2.1 地貌融合结果与分析

在视觉效果上，相较于TDL(图5(d))“阶梯状”地貌空间形态，NNDL_HPCS(图5(b))和NNDL_SAS(图5(c))两者与PNL(图5(a))较为相似，存在复杂的自然地貌纹理特征；在表土高程变化上，受土方平衡约束，NNDL_HPCS,NNDL_SAS,TDL三者较PNL高程变化幅度分别为-112～171,-205.90～50.03,-160.21～73.60 m，变化区间大小分别为NNDL_HPCS(283 m) > NNDL_SAS(250.93 m)>TDL(233.81 m)；在空间形变特征上，受“开采子区地表原始标高”、“复填子区煤层底板标高”、“复填子区煤层顶板等效标高”、“膨胀系数”4个要素的直接影响，本研究NNDL_HPCS与已有成果NNDL_SAS均较PNL发生空间上的不规则形变。例如，在东西方向上，复填子区～NNDL均较PNL有所抬升，而在其他区域均较原始地貌有所下降。与前两者相似，TDL较PNL也具有明显的空间形变特征(图5(d)，6(d))，然而与前两者不同，受斜坡设计影响，TDL与周边原始地貌间“分割”明显，具“斜坡-平台-斜坡”3段式的地貌形态特征(图5(d)，6(d))。

如图6所示，以研究区西边界为起点，以PNL为基准，由西向东逐列分别统计NNDL_HPCS,NNDL_SAS,TDL相对于PNL的高程方差。整体上，NNDL_HPCS,NNDL_SAS和TDL三者相对PNL的高程总方差分别为288 942.8,582 248.8和465 452.8 m。其中，NNDL_HPCS相对于PNL的高程方差最小，相较NNDL_SAS与TDL分别低50.37%与37.92%，在绝对高程上更接近于PNL，形态保留效果最优。随距研究区西边界距离变化，在0～810 m段，NNDL_HPCS,NNDL_SAS,TDL,PNL四者均位于矿坑位置，列相对高程方差值均为0；在810～1 800 m段和4 590～4 890 m段，NNDL_HPCS列相对高程方差远高于NNDL_SAS和TDL，其均值分别为后2者均值的11.83倍与6.99倍，和17.16倍与6.43倍；在1 800～4 590 m段，NNDL_HPCS列相对高程方差值则远低于NNDL_SAS和TDL，其均值分别为后两者均值的0.13倍与0.18倍。由地质资料可知，在810～1 800 m段和4 590～4 890 m段，随着前侧开采子区煤层底板走低，可采煤层厚度变化较少，导致后侧复填可用土方量不断增加，NNDL_HPCS,NNDL_SAS和TDL对应复填子区均形成“凸包地形”；在1 800～4 590 m段，随着前侧开采子区煤层底板走低，可采煤层厚度增大，导致后侧复填可用土方量稳中有减，NNDL_HPCS,NNDL_SAS和TDL对应复填子区均形成“凹陷地形”(图6)。即在“凸包地形”段，NNDL_HPCS与PNL栅格绝对高程重合效果上远弱于NNDL_SAS和TDL，在“凹陷地形”段则远优于2者。

图6 研究区PNL,NNDL_HPCS,NNDL_SAS,TDL三维示意(a)～(d)，及其复填子区物质等效厚度随距研究区西边界距离的折线统计(e)和3者相对于PNL的列高程方差随距研究区西边界距离的折线统计(f)Fig.6 Schematic diagram of the PNL，NNDL_HPCS，NNDL_SAS,TDL(a)～(d)，the broken line statistics chart of the equivalent thickness of the material in the refill subregion with the distance from the western boundary(e) and they relative to the PNL column elevation variance with the distance from the western boundary in the study area(f)

图7 研究区NNDL_HPCS,NNDL_SAS,TDL水文保留率随距研究区西边界距离的折线统计Fig.7 Hydrological retention rate of NNDL_HPCS，NNDL_SAS， and TDL in the study area is broken-line statistics with the distance from the western boundary of the study area

如图7所示，以研究区西边界为起点，以PNL为基准，由西向东逐列分别统计NNDL_HPCS,NNDL_SAS,TDL相对于PNL的水文保留率。整体上，NNDL_HPCS,NNDL_SAS和TDL三者水文保留率分别为73.53%,70.68%和65.00%(不包括矿坑部分)。其中，NNDL_HPCS水文保留率最高，相较NNDL_SAS与TDL分别高2.85%与8.53%，在边坡朝向上更接近于PNL，水文保留效果相较最优。随距研究区西边界距离变化上，与高程方差类似，在0～810 m段，NNDL_HPCS,NNDL_SAS和TDL三者均位于矿坑位置，水文保留率均为1；在810～1 800 m段和4 590～4 890 m段，NNDL_HPCS列水文保留率(59.19%和60.60%)低于NNDL_SAS和TDL，其均值分别较后2者均值低23.27%与9.31%，和11.84%与5.23%；在1 800～4 590 m段，NNDL_HPCS列水文保留率(79.44%)高于NNDL_SAS和TDL，其均值分别较后者均值高13.75%和16.52%。与上述高程方差类似，在“凸包地形”段，NNDL_HPCS与PNL栅格水文保留率弱于NNDL_SAS和TDL，在“凹陷地形”段则优于两者(图6)。

在土方运移上，依据采复子区土方对应关系，运用Matlab计算获取NNDL_HPCS,NNDL_SAS及TDL三者的土方平均运距(PNL不存在采挖，故不进行计算)。计算结果表明，3者土方平均运距NNDL_HPCS(863.23 m/m)>NNDL_SAS(822.20 m/m)>TDL(821.71 m/m)。其中，NNDL_SAS土方平均运距较TDL仅增加约0.06%(0.49 m/m)，NNDL_HPCS则较前者约增加5.05%(41.03 m/m)，其施工成本势必较前2者有所增加。

3.2.3 土壤抗水蚀能力求解结果与分析

地貌稳定性主要受气候、地质、植被等诸多要素的影响。其组成物料强度直接影响边坡的稳定性，是其抗侵蚀能力的主要指标。根据地质调查报告，确定研究区内排土场结构及岩土物理力学参数，并结合中国天气网(http://www.weather.com.cn/)获取研究区2007—2017年逐日气象参数。通过CLiDE演化模型，分别模拟NNDL_HPCS,NNDL_SAS,TDL,PNL四者10 a的土壤水蚀过程。其中，为保证模拟结果的可比性，使NNDL_HPCS,NNDL_SAS和TDL三者的采复子区空间位置保持一致，最终输入范围如图6所示。同时，为避免边缘效应对模型模拟结果的影响，将研究区周边1 000 m范围共同纳入运算(图6)。最终通过Matlab软件进行叠加对比分析，统计研究区内排土场NNDL_HPCS,NNDL_SAS,TDL和PNL四者10 a的土壤水蚀总量(表1)。

表1 研究区NNDL_HPCS,NNDL_SAS,TDL,PNL 土壤10 a总水侵蚀量Table 1 10-year total soil water erosion scale of the NNDL_HPCS，NNDL_SAS，TDLand PNL of the dump site

由表1可知，整体上，本研究内排土场10 a土壤水蚀总量在4.098 3×10～23.681 5×10m。其中，PNL土壤抗水蚀能力最强，模拟情景下其土壤10 a水蚀总量仅占4者平均水平(12.175 6×10m)的33.66%。在人工地貌设计下，本研究NNDL_HPCS和已有研究NNDL_SAS相较TDL分别可减少约39.07%和72.58%的土壤水侵总量，内排土场土壤抗水蚀能力提升明显。并且，相较于已有研究中NNDL_SAS土壤抗水蚀效果，本研究NNDL_HPCS可进一步减少约54.99%的土壤水蚀总量。

如图8，9所示，以研究区内排土场西边界为起点，以PNL10 a土壤水蚀为基准(将其列统计值定义为对应列的0值)，由西向东逐列分别统计NNDL_HPCS,NNDL_SAS,TDL相对于PNL的10 a列土壤水蚀量。与水文保留率和绝对高程方差的统计结果不同，在810～1 800 m段，10 a土壤水蚀量中NNDL_SAS(-6.78×10m)

图8 NNDL_HPCS,NNDL_SAS,TDL相对于PNL的土壤水蚀量随距研究区西边界距离变化Fig.8 Soil water intrusion of NNDL_HPCS，NNDL_SAS，and TDL relative to PNL varies with the distance from the west boundary of the study area

图9 NNDL_HPCS,NNDL_SAS,TDL和PNL的土壤水蚀量随距研究区西边界距离变化直方图Fig.9 Histogram of soil water erosion of NNDL_HPCS， NNDL_SAS，TDL and PNL with the distance from the western boundary of the study area

即整体上，研究区PNL在全段均发生较低水平的土壤水蚀；NNDL_SAS土壤水蚀主要发生在初始段(4 590～4 890 m段)与中段(1 800～4 590 m段)，土壤堆积主要发生在终段(810～1 800 m段)；NNDL_HPCS土壤水蚀主要发生在起始段与终段，土壤堆积主要发生在中段；TDL则在全段均发生较高水平的土壤水蚀。

4 讨论

煤炭作为我国重要的能源矿产，其露采比例随社会经济的发展不断增加。内排土场作为煤炭露采区的重要组成部分，其地貌形态不仅直接影响矿区人员设施及生产安全，也直接和间接影响区域生态环境系统稳定。受传统复填方式影响，内排土场往往表现为“梯田式”人工堆垫地貌，从而使原地表要素及水文特征丧失，并在集中降雨时容易引发严重的表土水蚀等问题。诸多实践与研究均表明，NbS理念下基于自然稳定特征的排土场地貌重塑具有重要的理论与实践意义，其在有效提高区域表土稳定性的同时所需维护费用较低。与上述结论相同，团队前置研究均表明，排土场近自然地貌的表土稳定性优于传统“梯田式”地貌。且前置研究中，基于SAS的采排复一体下内排土场NNDL较TDL 10 a水蚀量可减少约39.07%，初步验证了NNDL_SAS在理论上具备可行性。在此基础上，本研究进一步提出了基于HPCS的内排土场地貌近自然重塑方法。与前者类似，本文所提方法综合考虑了采排复一体化下土方平衡约束、区域坡度平滑控制与边界地貌融合3方面约束，在地貌融合效果上不仅优于NNDL_SAS，还进一步提高了研究区的抗水蚀能力(10 a土壤水蚀量分别较NNDL_SAS和TDL减少约55%和72.58%)。然而与前者不同的是，本研究基于HPCS的内排土场地貌近自然重塑方法不再仅着眼于地貌重塑结果与周边自然地貌间的衔接，而是致力于将复填可用土方再分配于原有的各个水文细小单元之中，在灰色系统理论支撑下通过保留或微调整此类要素以尽可能维护区域原有生态系统。

在露天矿采排复一体化过程中，已有坡面重塑方法能够较好地应用于传统内排土场复填区地貌重塑，并在缩短土方运距、减小占地面积、提升技术适用性等方向得到了充分发展。然而在地貌形态上，上述内排土场多呈现出“棱角分明”的“梯田式”空间布局，存在边坡表土易发侵蚀，景观明显破碎化等问题。虽然后续部分学者提出在原有施工工艺的基础上通过“弱化”排土场斜坡与平台间“分界”的方式提高区域表土稳定性及其融合效果，但实际上仍是“去除棱角”的“台阶状”地貌，未有效解决上述问题。并且，受采复周期、可采煤层厚度等诸要素的影响，复填区可用土方随采复周期动态变化，各周期间内排形态存在差异。本研究和前置研究中，笔者均采用采复子区一一对应的方式简化采复过程中复填可用土方的动态变化问题，并相较于TDL在土方运距小幅增加(0.06%～5.05%)的情况下大幅减少(39.07%～72.58%)表土水蚀量和获得更高的地貌融合效果。然而值得注意的是，相较于传统规则状内排土场地貌，NNDL形态较为复杂。在实际施工作业中，为避免近自然重塑地貌不规则起伏所带来的安全风险，应在排土场次高层平盘上运用卡车运输加单斗挖掘的形式进行近自然地貌重塑，其余台阶均应严格遵照传统内排土场地貌设计规范；同时应充分运用RTK等定位技术，确保排土场近自然重塑地貌特征要素(例如2.1.3节内沟道、坡定点)位置符合设计，通过“点线构面”的思路指导地貌近自然重塑。

受地貌发育影响，自然原始地貌并不总是稳定的(例如黄土高原冲刷地貌)。在特定条件下，传统复垦后矿区或将改良原有地貌。由于本研究HPCS模型以PNL为学习对象，因此在使用该模型前首先需保证PNL地貌稳定(例如本研究PNL 10 a的土壤侵蚀量最小，仅为TDL的28.4%)。并且，本研究各样本10 a土壤水蚀模拟量随距西边界距离变化的统计结果表明，在局部地区基于PNL的NNDL构建结果表土抗蚀能力甚至会劣于传统方式(图9中，1 800～4 590 m段，NNDL_SAS劣于TDL；4 590～4 890 m段，NNDL_HPCS劣于TDL)。因此，需充分论证NNDL,TDL和PNL三者的稳定性，从而确定有利于因地制宜的内排土场地貌重塑方式。