胜利一号露天煤矿北排土场土壤物理性质空间分布研究

王党朝 刘慧芳 肖 武,3 邹玉珠 王 鑫

(1. 神华北电胜利能源有限公司，内蒙古自治区锡林浩特市，026000；2. 中国矿业大学(北京)土地复垦与生态重建研究所，北京市海淀区，100083；3.浙江大学公共管理学院，浙江省杭州市，310058)

土壤是不均一且具有高度空间变异性的时空连续体。露天煤矿排土场主要是由剥离土层中的废弃物及煤层上覆的岩土层、煤粉和表土堆积构成。排土场的堆积侵占了大片土地，其结构的特殊性极易造成侵蚀和水土流失，在大风条件下，极易产生扬尘污染大气环境，因此采取因地制宜的方法，对排土场进行科学地复垦治理具有必要的现实意义。露天煤矿排土场受排土秩序和机械碾压的影响，土壤结构和层次等理化性质发生了较大变化，其质量也随之发生变化，土壤质量的高低是土地复垦成功与否的关键所在。

内蒙古胜利煤田矿区位于我国生态脆弱的干旱、半干旱地区，该地区植被覆盖率较低，水土流失和土地荒漠化现象严重。大型露天煤矿开采对原本脆弱的东部草原生态区造成了破坏，加剧了当地草原生态系统的进一步退化，已引起众多学者的关注。土壤物理性质主要包括土壤容重、孔隙度、土壤水力性质、土壤团聚体等，土壤的孔隙数量大小和容重的改变在较大程度上影响着土壤中水、肥、气、热等肥力因素的变化和供应状况。近年来，随着国家对露天矿区复垦工作的重视，陆续有学者对排土场土壤的空间变异性展开了研究。郭凌俐等人以传统统计学和地统计学相结合的方法，选取黄土区露天煤矿排土场复垦初期0～40 cm的土壤研究其土壤颗粒组成的空间变异规律，并运用普通克里格法进行最优无偏线性插值，制作了不同粒径土壤颗粒含量的空间分布图；瞿春艳等人利用地统计学方法分析了海州露天煤矿排土场土壤pH值、全氮、全磷、全钾等营养元素的空间变异性，分析得出这些指标均具有强烈的空间相关性，其空间变异性主要是由土壤母质、地形、气候等非人为因素引起的；吕刚和黄龙等人采用传统统计学、地统计学和GIS技术相结合的方法，以海州露天矿排土场边坡为研究对象，研究排土场土石混合边坡土壤抗冲性的空间变异特征及砾石对其影响；赵红梅等人对大柳塔采煤塌陷区土壤含水量的空间变异特征进行分析后认为，地表塌陷使得表层土壤结构变疏松，土壤粒度变粗，使其持水能力下降，土壤含水量比正常区域显著降低。

由此可见，在不同尺度、不同土地利用类型不同植被覆盖度的条件下，土壤的空间变异性均存在且各有差异。对于内蒙古胜利煤田矿区而言，由于其本身脆弱的生态环境，加上人为扰动，研究该地区排土场土壤的物理性状对于指导该矿区植被重建和生态恢复具有必要性。土壤特性空间变异性的研究方法己经由最初的经典统计学方法发展到了时序分析方法、地统计学方法、随机模拟、分形和分数维方法及其应用等研究方法，目前运用较多的是地统计学方法。胜利煤田胜利一号露天矿北排土场是该矿区最早的外排土场之一，将其作为研究对象，运用描述性统计法和地统计法分析土壤容重、孔隙度和含水量的空间分布，利用ArcGIS10.3绘制土壤物理性质空间分布图，并在此基础上分析原因，可为因地制宜地复垦治理整个胜利煤田的排土场提供理论基础，也可更好地指导土壤重构工作的顺利实施。

1 材料与方法

1.1 研究区概况

胜利煤田位于内蒙古自治区锡林浩特市北郊3 km，地理坐标为东经115°30′～116°26′，北纬43°57′～44°14′，地处内蒙古高原的中部，大兴安岭西延的北坡，属于NE-SW向高原、丘陵地形，地形标高970～1326 m，该地区内地貌形态由构造剥蚀地形、剥蚀堆积地形、侵蚀堆积地形和熔岩台地等4个地貌单元组成。本区属半干旱草原气候，春秋两季风大干旱，夏季炎热降雨少，冬季寒冷风雪多，年温差较大。每年6～8月为雨季，年平均降水量为294.74 mm，年平均蒸发量为1794.64 mm；冻结期为10月初至12月上旬，解冻期为翌年3月末至4月中旬，最大冻土深度为2.89 m。胜利一号露天矿在排土场绿化复垦上大胆尝试、不断摸索，经过几年的努力逐渐摸索出一套适合锡林浩特地区气候条件和排土场自然环境的绿化复垦“六步法”，即“一排、二覆、三沙障、四种、五灌、六养护”。北排土场于2006-2009年治理，占地面积为107万m2，排土场高度为60 m，4个台阶，绿化面积为101万m2，其中坡面35.5万m2、平盘面积为65.5万m2。主要有沙打旺、苜蓿、披碱草、柠条、羊柴、沙蒿和冰草等植物，均为抗旱耐旱的先锋物种。研究区位置如图1所示。

1.2 研究方法

1.2.1 采样布局

土壤采样点位布设一般遵循全面性、代表性、客观性、可行性和连续性的原则，布点方法则有简单随机法、网格法、分层法、地统计法和方差四叉树法。根据北排土场的特点，分别从边坡和平盘采集土壤样品。土壤样品的采集时间为2017年6月，采样月份天气干燥，没有出现大于5 mm的降水天气，采用网格法结合3S技术进行采样点布设，实际采样过程中对个别样点位置进行了调整，最终共采集95个样点土壤样品。所有样品采自0～20 cm土壤表层，利用环刀采集土样后放入铝盒，装入自封袋，采样同时记录样点编号、取样坐标等信息。

图1 研究区位置图

1.2.2 土壤物理指标测定

土壤样品带回实验室，烘干称重测定土壤含水量、容重和孔隙度。土壤含水量采用烘干称重法获得，容重采用环刀法获得，比重采用排水称重法获得，孔隙度采用土壤容重和比重计算获得。每个指标重复测定3次，取平均值。

1.3 数据处理

1.3.1 数据分析方法

(1)半变异函数。半变异函数(semi-variograms)也称为半方差函数，是地统计学中用以研究土壤空间变异性的关键函数之一，能够定量地反映土壤性质在不同距离上观察值之间的变化。这是一个连续变化的函数，反映了土壤性质的空间变异具有连续性。在一维条件下，当空间点x在一维x轴上变化时，区域变量Z(x)在点x和(x+h)处的值Z(x)与Z(x+h)差的方差的一半，定义为区域化变量Z(x)在x轴方向上的变异函数，记为y(h)，见式(1)：

(1)

式中：y(h)——半方差函数，h为样本间距又称为步长或位差；

N(h)——间距为矢量h的所有观测点对个数；

Z(xi)、Z(xi+h)——Z(x)在xi和xi+h位置上的实测值。

用GS+软件进行半方差函数拟合，其中C0称为块金值，它表示h值很小时，两点间取值的变化，反映了区域化变量Z(x)内部随机性的可能程度；A为变程，它的大小反映了研究对象中某一区域化变量的变化程度，也反映了影响范围；基台值由拱高C和块金值C0之和组成，反映某区域化变量在研究范围内变异的强度，通常使用块金值与基台值之比C0/(C0+C)或拱高与基台值之比C/(C0+C)，反映随机部分或自相关分占总空间异质性的大小。根据Cambardella等人的研究结果，当C0/(C0+C)<25%时，土壤物理性质有很强的空间依赖程度；当25%≤C0/(C0+C)≤75%时，土壤物理性质有中等程度的空间依赖程度；当C0/(C0+C)>75%时，相应的土壤物理性质有很弱的空间依赖。决定系数表示变异函数与模型拟合的程度，决定系数越大，模型拟合度越高。

(2)模型验证。通过平均误差ME、标准化平均误差MSE、均方根误差RMSE、平均标准差ASE和标准化均方根误差RMSSE这5个指标作为插值精度指标。如果预测误差是无偏的，则ME、MSE应该接近于零，如果RMSE接近于ASE，则可以得到正确的变异性，说明预测值的变异性与测量值相同；如果RMSE大于ASE，说明在预测中低估了变异性；如果RMSE小于ASE说明在预测中低估了变异性；如果RMSE小于ASE，说明在预测中会把变异性夸大。如果预测标准差是正确的，RMSSE应该接近于1。

(3)克里金插值法。本研究采用普通克里格法进行空间插值，即在半方差函数y(h)的基础上，用若干个采样点的实测Z(xi)来推求未知点的Z0(x0)的一种方法，通过采用线性加权法来表达未测点x0的预测值，见式(2)：

(2)

式中：Z0(x0)——未知点x0处的线形、无偏、最优估计值；

n——实测值个数λi为计算Z0(x0)时实测值Z(xi)的权重。

1.3.2 数据处理

在本研究中，土壤物理指标的描述性统计分析、K-S检验和相关分析过程均在SPSS19.0软件中完成。土壤中物理指标半变异函数模型的计算和选取采用GS+9.0软件进行，计算完成后将求取的模型参数输入到ArcGIS 10.3软件中，采用Geostatistical Analyst模块的Ordinary Kriging进行分析，最终得到土壤物理指标的空间分布图。

2 结果分析

2.1 物理指标描述性统计分析

为了掌握北排土场土壤容重、孔隙度和土壤含水量的整体变异情况，首先利用SPSS20.0对所选指标进行描述性统计分析，土壤物理性质描述性统计结果见表1。

由表1可以看出，排土场土壤的容重均值高达1.51 g/cm3，土壤紧实度远高于自然演变的土壤，这可能与排土场的形成过程受到机械压实有关。根据变异系数划分，当CV<10%时，变现为弱变异性；当10%100%时，表现为强变异性。容重的变异系数为10.28%，属于中等变异。

表1 土壤物理性质描述性统计表

孔隙度与土壤含水量是反映土壤水分物理特性的基本指标，能较好的反映土壤的疏松性、结构性以及持水性。排土场土壤质地多为砂质土，土壤孔隙度的变化范围从21.95%到51.65%，平均孔隙度为35.80%，土壤较紧实，疏松性差。变异系数为15.74%，属中等变异。

土壤含水量的变化范围为0.82%～13.30%，均值为4.04%，变异系数较大，为68.89%，属中等变异。研究区土壤含水量属于较低水平，复垦后土壤持水性能差，有待改良。

计算变异函数和进行Kriging插值时，要求所选土壤物理性质的变量数据均服从正态分布，对不符合正态分布的变量数据进行对数转换后再进行检验。从偏度值、峰度值以及K-S检验的P值综合考虑，土壤容重、孔隙度均符合正态分布，土壤含水量经对数转换后符合正态分布。

2.2 物理指标空间变异特征

传统的描述性统计能从整体上反映排土场复垦土壤质量，但不能定量化描述某些特定部位的土壤变化特征，因此选取合适的半变异函数模型对北排土场复垦后的土壤物理性质的空间变异性进行分析，北排土场土壤物理性质的变异函数理论模型及有关参数见表2。

表2 北排土场土壤物理性质的变异函数理论模型及有关参数

根据决定系数(R2)最大和残差平方和(RSS)最小的原则，容重和含水量选用指数模型拟合，孔隙度选用球状模型拟合。变程又叫做空间自相关距离，在空间自相关距离内土壤具有某种程度上的空间相关性。从表1中土壤理化性质的空间自相关距离结果分析，该小尺度研究区域内各指标之间的空间自相关距离变化较小，含水量的空间自相关距离最小为56.68 m，表明土壤含水量的空间异质性程度最大，影响含水量的因素较复杂。容重的空间自相关距离最大，为120 m。

结构性因素和随机性因素共同作用影响土壤理化性质的空间变异性。在变异函数模型中，块金值C0表示随机性因素(土地利用类型、人类干扰、动物活动)引起的变异，C表示由结构性因素(土壤母质、地形、气候)引起的变异，基台值C0+C是随机变异和结构变异之和。C0/C0+C表示空间变异性程度，如果该值较高，说明空间变异性主要是由随机性因素引起的，反之则由结构性因素引起的。比值小于25%，表明空间相关性很强，在25%～75%之间，表明具有中等空间相关性，大于75%说明空间相关性很弱；土壤含水量、孔隙度的C0/C+C0小于25%，表明这两种性质具有很强的空间自相关性。虽然排土场的复垦工程经历了很大程度的人工干扰，但其变异性主要还是由于土壤母质等结构性因素引起的，究其原因，可能是复垦所选用的土壤非原先表土，而由上覆岩层土壤组成，经实验测定，这些土壤质地偏砂，沙砾含量占比较大，持水性能较差。而孔隙度和含水量正是反映土壤结构性、持水性的重要指标。容重的C0/C+C0值为39.44%，介于25%～75%之间，属于中等空间自相关性，这说明变异是由随机因素和结构因素共同引起的，反映出复垦措施在某种程度上会影响土壤理化性质的空间变异性。

2.3 物理性质的空间分布格局及原因

在选取了合适的拟合模型的基础上，运用ArcGIS10.3的Geostatistical Analyst模块，选择Ordinary Kriging，采用交叉验证的方法对插值精度进行验证，并得到容重、孔隙度和土壤含水量的空间分布图。交叉验证结果见表3。

表3 交叉验证结果

由表3可以看出，容重、孔隙度和含水量的ME和MSE都接近于0，RMSSE都接近1，这表明3种模型分别模拟预测对应的物理性质的精度均较高，插值图的精度符合要求。从RMSE和ASE这2个评价指标来看，所选模型均高估了其对应的物理指标的变异性。

容重、孔隙度和土壤含水量的空间分布图如图2所示，三者的空间分布图之间存在明显的对应关系。

图2 北排土场土壤物理性质空间分布图

由图2可以看出，北排土场的中部地区容重值较大、孔隙度小且含水量小。这可能与复垦工序有关，中部地区是整个北排土场的最高处，且中间有一条主路，行驶车辆的反复碾压使得土壤紧实度增加、容重增大、孔隙度减小且水分含量降低。

由图2(a)中可以看出，北排土场的东北方向、正南方向容重值较高，西、西南、东、东南方向容重值偏低，其余部位为中等水平；图2(b)中所示的孔隙度空间分布图则表现出相反的趋势；图2(c)可以看出，北、西北、西、西南方向的含水量要整体高于东、东南、南方向，这可能是由于阳面坡与阴面坡相比较水分蒸发量大所致，与此对应阳坡植被比阴坡植被长势要好，可能使得阳坡植物对土壤的改良作用略优于阴坡植物，使得容重降低，孔隙度增大。由此可见，图中所反映的这3个物理性质之间关系的整体趋势与相关分析的结果一致，容重与孔隙度、含水量显著负相关，孔隙度与含水量呈显著正相关。

3 结果分析

通过研究胜利一号露天矿北排土场的土壤发现，排土场土壤的容重均值高达1.51 g/cm3，土壤紧实度远高于自然演变的土壤，这与排土场在形成过程中受到机械压实有关；土壤孔隙度的变化范围为21.95%～51.65%，平均孔隙度为35.80%，土壤较紧实，疏松性差；土壤含水量较低，保水持水性能较差。造成上述情况的原因，一方面是北排土场复垦年限较短，恢复期9 a左右，土壤内部结构恢复程度较低；另一方面，内蒙古地处干旱半干旱区，蒸发量较大，降雨量较小，且采样时间为6月初，还未开始进入雨季，所以土壤含水量属于较低水平。土壤孔隙度和含水量都具有很强的空间自相关性，受土壤母质等结构性因素影响较大，主要原因是经历采-排-复工序，原有土体结构被破坏，排土场复垦所用的土壤夹杂矸石、砾石等，结构复杂。土壤容重表现出了中等空间变异性，结合空间分布图和相关性分析可以看出，土壤容重与孔隙度、含水量呈负相关，土壤孔隙度和含水量呈正相关。

复垦过程中的机械反复碾压导致排土场容重增大，紧实度过大的土壤孔隙度必然变小，不利于水分的入渗又导致含水量较小。因此，结合排土场特点，宜采取物理措施和植物措施相结合的方法，通过深耕并配以合理的灌排设施，以调节土壤含水量，通过种植抗旱耐旱的先锋作物，向土壤中施肥增加有机质，以改良土壤结构。

4 结论

(1)北排土场的土壤容重、孔隙度和含水量均属于中等变异。与原状土壤相比，土壤容重均值偏高，孔隙度、含水量均值偏低。因此，矿区需要采取翻耕等疏松土壤的物理措施并结合植物措施进行土壤改良，可栽种一些耐旱型先锋植物，使土壤物理结构有利于植被的生长。

(2)通过半方差函数分析，土壤孔隙度和含水量的C0/C0+C值分别为8.19%和3.31%，具有较强的空间自相关性，说明这二者的空间变异主要是由土壤母质等结构性因素引起的；容重的C0/C0+C的值为39.45%，属于中等空间自相关性，这说明变异是由随机因素和结构因素共同引起的，反映出复垦措施在某种程度上会影响土壤理化性质的空间变异性。

(3)这3项物理指标均显示出了明显的空间分布格局，并且其分布格局在一定程度上具有对应性，这种对应性与相关性分析的结果是一致的。表明一种物理性质发生改变的时候，必然会导致其他的物理性质发生改变。所以，针对矿区特定的土壤状况，制定针对性的改良措施，是实现该地区植被恢复与生态重建的重要途径。