矿区排矸场地形重塑中自然地形要素的计算与分析——以宁夏宁东羊场湾矿区废弃地为例

蔡懿轩， 杨 刚†， 张成梁， 顾清敏

(1.北京林业大学信息学院，100083，北京； 2.国家林业草原林业智能信息处理工程技术研究中心，100083，北京；3.轻工业环境保护研究所，100089，北京； 4.国家能源集团宁夏煤业有限责任公司羊场湾矿区，751410，宁夏灵武)

矿区废弃地的生态修复中，地形重塑是第一步、也是最核心的部分。地形构建不合理将导致矿区严重的水土流失，使后期的植被重建受阻，从而导致极大的维持和养护成本。近年来，近自然地形重塑的理念和方法逐渐被提出并应用。大自然具有自我更新和再生能力，其自身原始的地形特征为生态重建提供最重要的参考依据。近自然地形重塑就是以矿区周边未被扰动的自然地形为参照来进行地形设计。因此，充分了解当地未被扰动的自然区域的地形、水文等特征已成为近自然地形重塑的先决条件。而在此过程中，存在一系列有待探讨的问题：1)参照区域的选择问题，即选择哪块儿区域作为近自然地形重塑的参照区域；2)参数的选择问题，即选择哪些地形特征参数进行计算和分析；3)参数的计算。

近年来，国内外许多学者对近自然地形重塑方法进行探讨[1-3]：Terrence等[4]针对矿区废弃地的坡面和沟道治理，以集水区为单元讨论地形重塑方法。杨翠霞等[2]探讨露天开采矿区废弃地近自然地形重塑的理论体系和技术方法，并结合Geofluv模型重建周口店采石场废弃地。Bugosh等[5]基于流域地貌学理论提出地形重塑模型GeoFluv。但这些工作重点放在地形的流域水力侵蚀作用上，且大多仅考虑流域尺度的参数计算。其对于参照区域选择、参数选择等问题尚未有系统工作，还未形成成熟的解决方案。本研究以宁夏宁东羊场湾矿区一区排矸场为实例，讨论近自然地形重塑中的参照区域确定、特征参数选取及参数计算问题，以给出可行解决方案，其分析结果可为该矿区排矸场近自然地形重塑提供参照目标，并为该地区地理分析、环境评价及生态环境治理等提供基础数据。

1 研究区概况

羊场湾矿区位于宁夏回族自治区灵武市宁东镇境内(图1a)，为我国宁东能源化工基地的主力矿井之一。该地区位于毛乌素沙漠的边缘地带，属温带干旱气候区，典型的大陆性气候，多年平均降雨量212 mm，年平均蒸发量2 862.2 mm，研究区全年盛行西北风，年平均风速3.1 m/s，全年>17 m/s的大风超过50 d。研究区地势较为平缓，属低缓剥蚀残丘地貌。研究区内土壤主要以淡灰钙土为主，部分地区有明显土壤沙质化现象，地表抗侵蚀能力弱；常年大风，风力对地表土壤的搬运作用和塑形作用明显。

羊场湾矿区平均年排矸量可达38万t。排矸场位于矿区东南处，主要依据地形自上而下倾倒的方式进行排矸，目前已形成约10 m高差的陡坡，并继续向东侧推进。矸石山坡度达40°，边坡易发生侵蚀现象；在极端气象条件下，易发生山体滑坡、塌方。随着矸石堆放的增多，对当地土地资源、大气环境、水资源等有强烈扰动。

2 研究方法

2.1 数据准备

本研究采用资源三号卫星拍摄5 m分辨率DEM数据为基础数据(图1b)。从资源三号卫星的一景数据中选取矿区周边区域87.33 km2的数据作为试验区，数据获取时间为2019年8月，地图数据范围为： E 106°32′～106°41′、N 37°57′～38°03′。图中待重塑区划定范围为羊场湾矿区一区排矸场区域。

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

2.2 参照区域的选取

在近自然地形重塑时，需根据矿区待重塑区域的情况选择与之有类似特征、对其重塑有参照价值的未干扰区域作为参照目标。

一般而言，流域侵蚀作用对地貌形成影响显著，与矿区的地质安全息息相关。研究区属西北干旱荒漠区，年平均降雨量很低，但其降雨非常集中，偶发性暴雨对地面冲蚀作用非常强烈，流域侵蚀作用依然是影响其地表形态的最重要因素。因此，在进行参照区域的选取和分析时，依据流域的邻近性、相似性等因素，并以流域为选取单位进行参照区域的确定；同时兼顾气候、土壤、岩性等其他自然因素。

具体而言，首先通过水文分析得到相关区域的流域信息，在此基础上依据如下原则进行参照区域选取：

1) 邻近性原则。选取待重塑区的邻近流域为参照区域；

2) 近似性原则。选择与待重塑区在水文、气候、土壤、植被和岩性等相似的流域区域；

3) 稳定性原则。选择处于内应力和外应力基本平衡的稳定状态的参照区域。

在此基础上，还需对所选取流域的相关特征进行筛查，剔除明显偏离整体样本数据的“异常”流域。其中采用箱线图方法[6]进行异常流域的剔除。

2.3 参数的选择

与地形特征相关的参数有许多[7-9]，但需重点选取对近自然地形重塑有参考意义的参数进行计算和分析。本研究主要参考杨翠霞等[10]给出的近自然地形重塑指标及美国GeofluvTM模型[11]中的地形计算参数。杨翠霞等[10]综合前人工作，选取流域面积、流域圆度、沟道密度等7项指标参数来刻画流域地貌的多方面形态特征，这些参数与矿区地形重塑设计密切相关，因此选取这些参数作为自然流域地貌特征的指标。

上述流域地貌特征重点在于流域宏观形态的把握，而在进行近自然沟道设计时，本研究重点参考了美国GeoFluvTM模型中有关沟道形态和局部地表形态的若干参数。GeoFluvTM模型是Bugosh等[5]开发的用于废弃矿区地形重建的数学模型。该模型可根据给定的沟道参数和地形参数来模拟生成近自然地形。

流域地貌特征是以流域为单元进行地形特征的描述，其测量尺度一般为数km2到数百km2。但地表更小尺度内的局部形状特征同样影响自然地形的稳定性，对近自然地形重塑具有重要参照作用。为此，选取坡度、坡向、局部地形起伏度、局部地表粗糙度等4项参数来刻画局部地表形态特征。其中坡度、坡向为坡面特征；局部地形起伏度、局部地表粗糙度分别刻画地形在竖直方向、水平方向上的变化特征。

综上所述，拟将测量参数归为3类，共14项参数：1)流域地貌形态：沟道级别、流域面积、主沟道长度、沟道密度、流域圆度、流域高程差、平均坡度；2)沟道形态：蜿蜒度、沟道比降、分水岭到沟道源头最小距离；3)局部地表形态：坡度、坡向、局部地形起伏度、局部地表粗糙度。

2.4 参数的计算

拟计算的地形参数涉及地形多种类型的形态特征，需采用不同的计算方法，且有计算上的先后次序。为此我们制定系统的计算路线(图2)：首先进行水文分析并确定参照区域；其次分2个分支进行参数计算，左侧分支进行流域地貌特征和沟道形态特征的计算，右侧分支对局部地表形态特征进行计算。

图2 参数计算的技术路线Fig.2 Technical route of parameter calculation

在“汇流量阈值的确定”步骤中首先拟合出沟道密度与汇流量阈值间的关系曲线，其次采用均值变点分析法[12]找到拐点来确定河网汇流量阈值。均值变点法是目前最佳分析窗口计算中的主流方法，相比人工作图法等方法，其分析效果更符合实际。在3.3.3节中局部地形起伏度计算中，采用均值变点法来确定最佳分析窗口尺寸。

“分水岭到沟道源头最小距离”这项参数来自于GeoFluv模型[11]。当降水发生后，坡面上部(靠近分水岭的部分)不会立刻形成侵蚀的沟道，而是需要经过一定距离的流量汇集，对坡面形成足够的侵蚀力形成沟道。参数“分水岭到沟道源头最小距离”就是为反映其现象。该参数与降水量、降水强度、地表粗糙度、土壤的抗冲性和渗透性都有关系，是地形重塑中近自然沟道设计的重要参考值[11]。首先基于地形位置指数[13]来识别地形类别，由此确定分水岭区域；其次针对每个子流域确定入水口；之后利用ArcGIS测量分水岭到沟道源头最小距离值。

3 结果与分析

3.1 流域分析及参照区域的确定

3.1.1 汇流量阈值分析及河网提取 以2000栅格为间隔，提取不同汇流量阈值下的河网，并计算其沟道密度，得到沟道密度与汇流量阈值间的关系曲线(图3)。

图3 沟道密度与汇流量阈值拟合关系Fig.3 Fitting relationship between gully density and catchment flow threshold

从图3中看出，随着汇流量阈值的增大，沟道密度逐渐递减趋于平缓。其拟合曲线的相关性系数R趋于1，符合幂函数关系。根据均值变点分析法[12]，计算原始样本的统计量S与样本分段后的统计量Si，S与Si之间最大差值对应的点即为变点。观察S与Si的差值变化曲线(图4)，当阈值为4 000时差值达到最大，因此确定河网提取的最佳汇流量阈值为4 000。

利用ArcGIS得到阈值4 000河网提取结果(图5)。其中河网中河道总长度为228.04 km，河道数目为438条。

图5 河网提取结果(汇流量阈值为4 000)，粉红色框区域则是3.2节分析后所选定的参照区域Fig.5 River network extraction results (catchment flow threshold is 4 000), the pink framed area is the reference area selected after the analysis in Section 3.2

3.1.2 流域分析及参照流域的确定 根据图5给出的待重塑区及其周边的沟道状况，对沟道进行分级[14]。图中排矸场西侧存在部分人工区域，因此主要从排矸场东侧的未扰动的自然流域中选择参照流域进行分析。

根据邻近性原则，重点考虑与待重塑区邻近并相交的2个流域区域：待重塑区北部南向北的流域(图6中紫色流域1)与待重塑区南部从东向西的流域(图6中绿色流域2)。待重塑区北部所覆盖的Ⅰ级沟道由南向北汇合成1个Ⅱ级沟道，并与右侧的另一个由东南向西北流的支流汇入流域1的主沟道(图中黄色虚线沟道)。而待重塑区南部的Ⅰ级沟道由北向南汇合成1个Ⅱ级沟道，并汇入到流域2主沟道(图中的蓝绿色虚线沟道)。待重塑区所覆盖的都是Ⅰ级沟道，并与流域1和流域2相联通。由于位置邻近，且这2个自然流域与待重塑区在气候、土壤、植被和岩性等方面都具有很大的相似性，依据相似性原则，选择2个流域内含Ⅰ级沟道的流域作为参照流域(图6)。

图6 分级水文网与Ⅰ级子流域划分示意图Fig.6 Schematic diagram of classification of hierarchical hydrological network and level Ⅰ sub-watershed

图6为初步确定的29个Ⅰ级子流域区域。其中子流域1、3和4与待重塑区域相交，因此剔除不作为参照流域。将除上述3个子流域外的26个Ⅰ级流域作为初步选定的参照流域，26个子流域的6项参数计算结果见表1。

这26个子流域中，有些子流域与其他流域形态特征差别很大，对数据统计造成干扰；为此通过箱型图方法[6]对子流域6项特征参数进行统计，剔除“异常”流域，最终筛选出21个子流域作为参照流域。在此基础上，将这21个子流域区域所属的2条主流域(流域1和流域2)所辖区域作为参照区域(图5中的粉红色框线区域)进行后续局部地表形态特征的计算和分析。

3.2 流域地貌形态特征及沟道形态特征

采用SPSS软件对表1中21个参照流域的6项流域地貌特征参数进行S-W检验，其显著性>0.05。表明参数分布都较为集中，符合正态分布，可作为地形重塑参考目标。21个流域的参数均值为：平均子流域面积0.25 km2；平均主沟道长度0.54 km；平均沟道密度2.29；平均流域圆度0.27；流域高程差32.33 m；平均流域坡度4.61 °。

表1 初步选定的26个Ⅰ级子流域地貌形态特征Tab.1 Morphological characteristics of the 26 level Ⅰ sub-watersheds initially selected

进一步地，对参照流域内的23条Ⅰ级沟道(6、24号子流域内分别有2条 Ⅰ 级沟道)沟道形态特征参数进行测量(表2)。

根据河流分类标准[15]，当蜿蜒度S=1.0时，为顺直河段；S< 1.2时为低度蜿蜒，S=[1.2，1.4]为中度蜿蜒，S>1.4为高度蜿蜒。从表2中得出目标流域内沟道的蜿蜒度范围为[1.030，1.316]，其中57%的沟道蜿蜒度<1.2；35%的沟道蜿蜒度在[1.2，1.3]之间，说明大部分沟道属于中低度蜿蜒沟道。参照流域的沟道比降在[0.001 0，0.045 1]之间，其中52%的沟道比降<0.02；30%的沟道比降在[0.020 0，0.030 0]之间；除11号沟道，沟道比降都<0.040 0。

表2 Ⅰ级沟道形态参数Tab.2 Level Ⅰ gully shape parameters

对这3项参数数据进行S-W检验，3项参数显著性> 0.05，表明数据分布集中，服从正态分布。此外，对3项参数进行箱线图分析，发现无异常值。因此，可将此3项参数的均值作为近自然地形重塑设计的参考目标。

3.3 局部地表形态特征

3.3.1 坡度 利用ArcGIS对参照区域进行坡度分析，根据我国水力侵蚀强度分级指标中面蚀分级标准，将坡度分为6个等级生成面蚀分级标准下的坡度分级图(图7)。从图中可知，自然地形中大部分区域为≤5 °的不易发生面蚀的平地和微坡面，面积占72.74%；而>15 °的陡坡区域仅占1.42%。

3.3.2 坡向 坡向利用ArcGIS对参照区域进行坡向分析，生成坡向分级占比图(图7b)。根据图中的占比数据统计可知，自然地形的主要坡向为北、东、西、东北、西北，共占区域面积的70.41%。若按照阴坡、阳坡、半阴坡、半阳坡对坡向重分类，比例分别为30.18%、19.34%、24.45%和26.04%；阴坡面积明显大于阳坡。

图7 局部地表形态特征参数分布图(括号中的百分数为该等级区域占整体面积的比例)Fig.7 Distribution map of local surface morphological characteristic parameters (the percentage in parentheses after the level is the percentage of the area of the level in the overall area)

研究区的主风向为西北风，而参照区域中偏北、偏西坡向所占比例较大。从图中可看出许多子沟道的分水岭呈东北- 西南走向，因而形成许多垂直于西北风向的坡地。这种垂直于主力风向的山脊地形对防风固沙有一定作用，这对于地形重塑有一定参考意义。

3.3.3 局部地形起伏度 局部地形起伏度指局部窗口内地形高程最大值与最小值之差。为确定地形起伏度最佳分析窗口，依次计算网格2×2、3×3……10×10下的平均地形起伏度，通过分析地形起伏度与网格面积间的关系曲线，采用均值变点分析法确定7×7网格为最佳分析窗口，在该网格窗口下的地形起伏度情况见图7c。从图中可知，比例最大的为(1,2] m起伏度的区域(35.01%)；≤3 m的地形区域占72.1%；而>5 m的区域仅占7%。说明整体地形起伏度较小，地势平坦。

3.3.4 局部地表粗糙度 局部地表粗糙度可反映地形在水平方向上的复杂度，是衡量地表侵蚀程度的重要量化指标之一。在ArcGIS中计算得到分级地表粗糙度(图7d)。从图中可见，自然地形96.63%区域的地表粗糙度在[1.000,1.019]之间，粗糙度较低，表明侵蚀程度较小。

4 结论

1)流域侵蚀作用依然是影响西北干旱荒漠区地表形态的最重要因素，因此在参照区域选择上，提出以流域为选取单位，依据流域的邻近性、相似性，并兼顾气候、土壤、岩性等其他自然因素的策略。

2)从流域地貌形态、沟道形态和局部地表形态3方面出发，选取14项参数作为刻画近自然地形特征的因子。这3方面的参数综合近自然地形重塑对不同尺度(流域尺度、局部地形小尺度)以及不同地形构建单元(流域、沟道、坡面)的参数需求，以满足近自然地形重塑的设计之需。

3)基于研究区流域分析，选取21个I级子流域作为参照流域。统计得到这21个参照子流域的流域地貌特征参数值及沟道特征参数值的分布比较集中，显著性> 0.05，表明这些参数值有代表性，可作为地形重塑的参照目标。