晋北露天煤矿区重构土壤容重差异研究： 探地雷达识别与反演

徐韩笑，曹银贵,2①，罗古拜，白中科,2

〔1.中国地质大学(北京)土地科学技术学院，北京 100083；2.自然资源部土地整治重点实验室，北京 100035〕

露天开采是国际上主要产煤国家资源开发的重要方法，其产量占比一般在50%以上，有的国家甚至超过90%[1]。煤炭资源的大规模开采，一方面满足经济建设需要，另一方面给局部地区带来众多生态环境和社会问题[2-4]。露天矿区在开采过程中，产生大量剥离物，开采形成的排土场是露天矿区最主要高生态风险区，极易导致人为水土流失，破坏矿区生态系统，进而引发各种环境地质灾害。因此，矿区排土场复垦对露天矿区可持续发展有不可替代的意义[5-7]。在矿区土地复垦过程中，土壤物理性质对矿区土地复垦效果具有显著影响，对植被重建和生态恢复起决定性作用[8-11]。土壤容重是土壤重要物理性质，是衡量土壤紧实度和土壤肥力的重要指标之一，对研究土壤持水保肥能力具有不可或缺的意义[12]。

复垦过程中，矿区排土场土壤容重因机械反复碾压极易增大[13-18]，土壤紧实度过大，孔隙度变小，进而影响植被生长,影响矿区复垦进程。同时，重构土壤典型物理性质空间分异特征明显，而传统土壤采样方法的采样深度和采样数量受复垦地障碍限制，给科学揭示土壤空间分异特征及形成机制带来一定困难，并且传统土壤采样方法对土壤破坏性大，操作繁琐，不适于大范围使用[19-21]。探地雷达是近年发展起来的一种探测地下目标的无损探测技术，与其他常规地下探测方法相比，具有探测速度快、探测过程连续、分辨率高、操作方便灵活等优点[22-23]。目前，探地雷达在国内外土壤探测方面的应用逐渐展开，而在矿区尤其是大型露天矿区重构土壤探测的应用较少，特别对矿区重构土壤容重的研究较少。已有研究主要针对探地雷达无损探测技术的理论层面[23-24]，通常采用实验室设计实验构造出具有单一变量的样本。有研究[25-26]表明土壤体积含水率对介质的介电常数影响较大，只有在土壤体积含水率低于0.5%条件下，介电常数与介质容重呈显著正相关。但在实际应用过程中，土壤体积含水率低于0.5%的情况极少，而针对其他土壤含水条件下土壤容重与介电常数相关关系的研究较少[27]，导致探地雷达在探测土壤容重方面应用推广的理论依据不够充分。

因此，以中煤平朔安太堡露天矿区南排土场、西排土场、西排土场扩大区和内排土场为研究区域，以复垦土地重构土壤容重数据为基础，通过对比分析法、显著性分析法和拟合分析法对黄土露天矿区复垦地重构土壤容重变化规律进行研究，分析黄土露天矿区重构土壤容重总体特征和深度特征，以及不同重构土壤容重在探地雷达图像上的差异，验证探地雷达技术应用在重构土壤容重研究中的可行性，建立基于探地雷达图像的重构土壤容重-介电常数拟合模型，旨在为利用探地雷达图像进行快速无损探测重构土壤容重及精准复垦提供有力支撑。

1 材料与方法

1.1 研究区概况

平朔矿区位于山西省朔州市，地处北纬39°23′～39°37′、东经112°10′～113°30′，矿区总面积为380 km2，属于温带半干旱大陆季风性气候区，地带性植被类型属于干草原植被(图1)。研究区是中煤平朔安太堡露天矿区的排土场，海拔高度为1 360～1 465 m，地带性土壤主要为黄棉土、栗钙土和红黏土，呈微碱性和碱性，土壤有机质含量低，结构差，抗蚀能力弱。研究区于1985年开始进行初级土地复垦工作，迄今已经持续30余年。南排土场属于土石混排类型，其中粒径≥50 mm的岩石占46.98%，粒径为5～<50 mm的砾石占15.48%，粒径<5 mm的砾石、土砂占37.54%，采用“草-灌-乔”的植被配置模式进行复垦[28-30]。西排土场、西排土场外扩区和内排土场分别于2000、2005和2010年相继开始复垦工作，主要复垦为草地和林地，种植作物为柠条、紫花苜蓿、小叶杨、刺槐、榆树和油松等。

研究区地处黄土高原东部，冬春季节干旱多风，夏秋季节降水集中，年蒸发量超过降水量的4倍[30]。区域土壤物理风化强烈，侵蚀严重，较为贫瘠，土体干旱，土质偏沙，是黄土高原典型生态脆弱区[14]。在此地区进行大规模露天开采，更是加剧了水土流失。

1.2 样方设置与剖面建立

综合考虑研究区气候与土地覆盖情况，于2018年5月在中煤平朔安太堡露天煤矿南排土场西侧共设置5个样点进行实地采样,于2019年5月在内排土场、西排土场和西排土场扩大区分别设置3个样点进行实地采样，并采用GPS记录样方中心点坐标，同时记录样点高程、地形、坡度、坡向以及植被类型和生长状况。在4个排土场分别选取3个样点，尽量保证选取的12个样点处于近似高程，同时环刀取样法只适用于纯土壤剖面，因此只选择纯土壤剖面样点，选取的样点为南排土场(S1、S3、S4)、西排土场(W1、W2、W3)、西排土场扩大区(WK1、WK2、WK3)和内排土场(N1、N2、N3)，具体分布情况见图1。

排土场平台倾斜角度在0～2°之间，直接在平台上设置10 m×10 m的大样方；在排土场边坡上的大样方宽度设置为10 m，样方长度根据边坡坡度确定，使其水平方向上也为10 m。在大样方内随机设置1个1 m×1 m的小样方，清理土壤表面植物残留物和杂质，在小样方内制作土壤剖面，剖面平均深度为60 cm。采用环刀(体积为100 cm3)在每个土壤剖面上分层采集0～10、>10～20、>20～30、>30～40、>40～50和>50～60 cm土层样品，同时对各土层样品进行编号、保存。对采集的土层样品立即称重，并记录土样湿重。回到实验室后将土样放入105 ℃烘箱烘烤8 h至恒重，重复3次，计算土壤容重。土壤容重计算公式为

(1)

式(1)中，ρ为土壤容重，g·cm-3；V为环刀体积，100 cm3；m为环刀质量，g；M为环刀和烘干土样总质量，g。

土壤孔隙度为土壤中孔隙体积与土壤体积的比值，其计算公式为

(2)

式(2)中，n为土壤孔隙度；2.65为土壤密度平均值，g·cm-3。

1.3 探地雷达参数确定与数据处理

探地雷达是利用电磁波的反射来探测地下介质的一种无损探测技术。探地雷达主要由主机和天线组成，发射天线发射的电磁波在传播过程中遇到介电常数存在差异的介质时，电磁波在2种介质的分界面处发生反射和能量衰减，使电磁波双程走时、振幅和波形产生相应改变，进而解译出目标体形态、空间位置和结构[19-20]。土壤剖面挖好后，在探地雷达探测面上，将管径3、4和6 cm的钢管分别打进剖面的不同深度，然后在地表利用探地雷达进行探测。使用探地雷达对大样方和小样方分别进行探测，探测时探地雷达设定电磁波在土壤中平均传播速度为0.1 m·ns-1，介电常数为9。探测过程中采用整体增益和分段增益2种方式处理雷达数据，首先将整体增益调整到20 db,测量3次后取最优结果，然后将分段增益调整到合适区间，同样探测3次后取最优结果。所用探地雷达型号为LTD-2600，频率为900 MHz，由中国电波传播研究所提供。采用IDSP 7.0雷达图像处理软件处理雷达图像数据。

根据钢管标定深度和传播时间，确定电磁波在不同土壤层中传播速度。介电常数计算公式为

(3)

式(3)中，ε为介电常数；c为光速，3.00×108m·s-1；v为电磁波信号在介质中传播速度，m·s-1，即土层深度与电磁波传播时间的比值。

2 结果与分析

2.1 重构土壤容重总体特征及变化规律

将野外调研所得数据整理分析得到不同排土场重构土壤容重，进一步计算其均值与标准差(表1)。

表1 采样点重构土壤容重及土壤孔隙度

W1样点重构土壤容重均值最小，为1.20 g·cm-3；N2样点重构土壤容重均值最大，达1.75 g·cm-3，两者相差0.55 g·cm-3。同一排土场不同样点重构土壤容重均值之间存在差异，但差异较小；不同排土场重构土壤容重均值之间存在明显差异，尤其以西排土场(W1)和内排土场(N2)之间差异最大(图2)。

采用标准差表示样点不同深度重构土壤容重差异性。在0～60 cm土层，S1样点标准差(0.24)最大，这是因为S1样点0～10 cm土层重构土壤容重较小，远小于样点平均值(1.51 g·cm-3)，而>10～60 cm土层重构土壤容重明显增加，尤其是>20～50 cm土层重构土壤容重均超过1.55 g·cm-3，该样点不同土层重构土壤容重差异明显；WK3样点重构土壤容重标准差(0.03)最小，这表明该样点不同深度重构土壤容重无明显差异。

由图2可知，分析不同排土场不同土层重构土壤容重变化发现，各排土场0～10 cm土层重构土壤容重均值由大到小依次为内排土场、西排土场扩大区、西排土场和南排土场。各排土场>10～20、>20～30、>30～40、>40～50和>50～60 cm土层重构土壤容重均值由大到小依次为内排土场、西排土场扩大区、南排土场和西排土场。

2.2 利用探地雷达图像定性分析重构土壤容重差异

土壤容重与土壤孔隙度有关，土壤孔隙度越大，容重越小。土壤孔隙度越大，土壤中空气越多，土壤对电磁波的吸收越不明显，电磁波振幅也越大。同时，电磁波在传递过程中，遇到不同质地时，探地雷达信号会出现明显振幅反射[21]。可通过土壤孔隙度不同导致探地雷达信号图上的差异，反演重构土壤容重的差异。由表1可知，样点W1土壤孔隙度均值最大，样点N2土壤孔隙度均值最小，两者相差0.21。所有样点土壤孔隙度标准差均小于0.10，其中样点S1土壤孔隙度标准差(0.09)最大，样点WK3(0.01)最小，说明样点WK3土壤均一性较高。所有样点中只有样点S3和W1土壤孔隙度均值大于0.50，标准差分别仅为0.03和0.02，说明这2个样点土壤相对疏松。

采用IDSP 7.0雷达图像处理软件对各样点雷达图像进行校正零偏、调节零点、调整增益、整体法背景消除和滑动平均等处理，最后进行RGB调色，颜色越深表示探地雷达信号振幅越大(图3)。大振幅信号越多，土壤孔隙度越大，重构土壤容重越小。在0～10 cm土层，南排土场雷达图像(对应时间为0～1.25 ns)与内排土场(对应时间为0～2.00 ns)之间差异显著，南排土场雷达图像具有更多大振幅信号，而内排土场雷达图像颜色相对一致。在10～30 cm土层，西排土场雷达图像(对应时间为2.50～4.20 ns)与西排土场扩大区(对应时间为2.50～4.50 ns)和内排土场(对应时间为2.10～4.50 ns)之间差异显著，相对于西排土场扩大区和内排土场，西排土场雷达图像有较多大振幅信号，而西排土场扩大区和内排土场雷达图像颜色相对均一。

2.3 介电常数反演重构土壤容重

在内排土场(N1、N2和N3)、西排土场扩大区(WK1、WK2和WK3)以及西排土场(W1、W2和W3)各样点分别采用直径为3、4和6 cm钢管测定土壤介电常数，在南排土场(S1、S3和S4)各样点使用直径为4和6 cm钢管测定土壤介电常数。在雷达图像振幅原始数据基础上，提取电磁波传播走时(图4)。39组实测重构土壤容重数据见表2。

表2 不同深度探地雷达探测参数

由于土壤体积含水率差异对介电常数影响较大，选择土壤体积含水率小于15%的19组重构土壤容重实测值，采用式(3)计算得到不同土层的介电常数，并利用各土层重构土壤容重实测数据与对应介电常数构建拟合关系模型，拟合分析结果见图5。

通过拟合关系可知，不同排土场重构土壤容重与相应介电常数均呈正相关关系，其拟合关系R2=0.508，P=0.001，拟合结果可信。将19组介电常数代入拟合关系模型反演重构土壤容重(表3)。由表3可知，模型反演所得重构土壤容重略高于实测值，由2种方式得到的重构土壤容重差值绝对值最大为0.227 g·cm-3，最小为0.01 g·cm-3。

表3 重构土壤容重实测值与模型反演结果对比

3 讨论

3.1 不同排土场重构土壤容重差异原因分析

研究区不同排土场重构土壤容重大小差异明显。有研究[14,30-31]表明，重构土壤容重会随复垦年限的增加而呈下降趋势。笔者研究中，不同排土场0～60 cm实测重构土壤容重均值由大到小依次为内排土场、西排土场扩大区、南排土场和西排土场，南排土场在复垦时，0～20 cm深为表土，20 cm以下土壤均质性差，重构土壤容重明显增加，同时植被类型为刺槐和柠条，根系主要分布于0～25 cm土层，导致40～60 cm土层重构土壤容重较大，使得南排土场样点重构土壤容重均值大于西排土场。就不同排土场同一深度重构土壤容重而言，0～10 cm土层重构土壤容重均值随复垦年限增加而减小，由大到小为内排土场、西排土场扩大区、西排土场和南排土场；这是由于随复垦年限增加，地表植被生长发育并进行自然演替，导致表土细根逐渐增加，特别是草本植物根系非常密集，植物残渣丰富，增大了土壤孔隙[30]，从而降低了重构土壤容重。各排土场>10～20、>20～30、>30～40和>50～60 cm土层重构土壤容重均值由大到小为内排土场、西排土场扩大区、南排土场和西排土场。内排土场各土层重构土壤容重均值均最大，这是由于内排土场复垦年限最短，植被生长演替情况较差，土壤由粉(砂)壤土和砂质壤土组成，因土地在平整过程中被大型机械车辆压实，导致土壤呈片状、块状，自然条件下容重发育消失[30]，且后期复垦容重修复较慢所致。

3.2 重构土壤容重的探地雷达图像识别

根据雷达图像上的大振幅信号或颜色深浅差异可以定性分析土壤质地均一性，并大致显示出某范围内重构土壤容重差异。电磁波在遇到2种电性差异较大的介质时，在2种介质分界面上会发生反射和能量衰减，这在探地雷达信号图上表现为大振幅信号[27]，这为基于探地雷达信号图定性分析重构土壤容重提供了理论依据。在土壤体积含水率较低条件下，不同深度重构土壤容重不同时，介质的电导率或介电常数也不同。

笔者研究中，为减少土壤体积含水率对介电常数的影响，选取土壤体积含水率较接近的典型样点进行定性分析。W2样点>10～20、>20～30和>30～40 cm土层土壤体积含水率均为16%左右且差异较小，大振幅信号数量较少，雷达信号分布较均匀，实测结果也显示重构土壤容重差异较小；W3样点>10～20、>20～30和>30～40 cm土层土壤体积含水率均为15%左右且差异较小，而>20～40 cm土层雷达图像有较多大振幅信号，实测结果也表明>20～40 cm土层重构土壤容重存在较大差异；WK1样点>10～20、>20～30和>30～40 cm土层土壤体积含水率均为18%左右且差异较小，而>20～40 cm土层雷达图像有较多大振幅信号，这与该土层重构土壤容重实测值差异较大一致。综上所述，利用雷达图像对重构土壤容重差异进行定性分析时，大振幅信号与重构土壤容重差异之间存在相关关系，大振幅信号较多时，重构土壤容重差异较小。但土壤体积含水率或土壤质地差异较大时，利用该原则得到的结论有一定偏差，需要根据实际情况进行调整。

3.3 重构土壤容重-介电常数拟合效果评价

在影响介电常数的因素中，土壤体积含水率对介电常数的影响较大[25-26]，因此笔者所选采样时间均为非雨季，且尽量选择处于非饱和状态的土壤，使重构土壤容重成为影响电磁波能量衰减的重要因素之一，避免土壤体积含水率差异过大对试验结果产生影响。拟合结果显示，当土壤体积含水率小于15%时，重构土壤容重(y)与介电常数(x)之间呈正向线性相关，拟合函数为y=0.054 3x+1.019 2。笔者通过多次实地测试发现，土壤体积含水率小于15%条件下，重构土壤容重对土壤介电常数的影响仍十分明显。由拟合结果可知，土壤体积含水率小于15%条件下，重构土壤容重随介电常数的增大而增大，拟合关系R2为0.508 3，P值为0.001，这说明对于晋北露天煤矿区不同复垦排土场，重构土壤容重与介电常数均存在正相关关系，这与其他对介电常数与土壤容重关系的研究结论[24]一致。

笔者在采集相关数据时，对土壤质地影响的考虑不足，存在某一采样土层砂质土与黏质土混合的情况，这增大了拟合曲线与拟合函数的误差，未来研究可以考虑分别采集不同土质数据，以提高拟合函数精确性。同时，笔者研究未深入讨论土壤体积含水率高于15%条件下重构土壤容重与介电常数之间的关系，如何通过技术手段有效减少过高体积含水率对介电常数的影响，使重构土壤容重与介电常数拟合关系误差减小，这也是未来的研究方向。

3.4 探地雷达在反演重构土壤容重方面的可行性

早期有研究人员采用电磁波波速和振幅反演土壤容重，认为土壤容重与电磁波波速呈反比关系，与电磁波最大振幅呈负相关关系,且建立了土壤容重-介电常数数学模型，土壤容重和介电常数呈正相关关系[24]。综上所述，在土壤非饱和状态下，利用介质电导率参数可在探地雷达信号图上定性分析重构土壤容重相对差异，利用介质介电常数-重构土壤容重模型可定量分析重构土壤容重真实值。但多数研究只停留在实验室验证阶段，实际应用较少，尤其在针对矿区复垦地重构土壤这类复杂情况时，基于探地雷达的重构土壤容重分析往往达不到预期效果。笔者研究结果表明，当土壤体积含水率差异较小时，通过分析雷达大振幅信号分布情况可以判断不同深度重构土壤容重之间是否存在较大差异，雷达图像有较多大振幅信号时，说明该深度重构土壤容重差异明显，大振幅信号越多，重构土壤容重差异越大。

当土壤体积含水率小于15%时，探地雷达可实现对重构土壤容重的定量分析，利用雷达图像上电磁波传播到设定深度处的双程走时，可以得到电磁波在该土壤深度范围的传播速度，将传播速度代入经验公式即可得到相应介电常数，再根据笔者所得拟合关系，即可得到所测范围重构土壤容重估测值。笔者研究思路可用于与研究区具有相似气候特点的其他露天矿区土地复垦质量检测工作中。但在推广使用时，对于土壤质地存在较大差异和土壤体积含水率大于15%的地区，该方法具有局限性，利用该方法得到的重构土壤容重估测值可能与实际值存在较大差异。

4 结论

(1)晋北露天矿区不同复垦排土场及同一剖面不同深度重构土壤容重差异明显。各排土场0～10 cm土层重构土壤容重由大到小为内排土场、西排土场扩大区、西排土场和南排土场，重构土壤容重随复垦年限增长而增加。

(2)针对晋北露天矿区不同复垦排土场，根据雷达大振幅信号分布情况可定性分析重构土壤容重相对大小。电磁波大振幅信号越多，土壤孔隙度越大，重构土壤容重越小；雷达图像颜色越均一，说明该范围重构土壤容重差异越小。

(3)针对晋北露天矿区不同复垦排土场，在土壤体积含水率小于15%条件下，可通过土壤介电常数反演重构土壤容重，重构土壤容重随着介电常数增大而呈线性增加趋势。