黄土露天矿区重构土壤体积含水率表征与反演

罗古拜，曹银贵，2①，白中科，2，况欣宇，王舒菲，宋 蕾（.中国地质大学（北京）土地科学技术学院，北京 00083；2.自然资源部土地整治重点实验室，北京 00035）

中国黄土露天矿区位于降雨少、蒸发量大、地下水水位深的干旱半干旱区。由于机械压实、结构破坏和水土流失等问题，矿区排土场复垦土壤达不到原有土地的生产力［1-2］，并且复垦地土壤质量需要长时间恢复，复垦后13 a土壤也不能恢复到未扰动土壤质量水平［3］。土壤体积含水率是制约复垦土壤质量的重要因素［4］。在立地条件较差的排土场，土壤体积含水率直接决定植被的生长状况，所以对矿区复垦地土壤体积含水率的研究尤为重要。有研究指出，复垦后土壤持蓄调节水分能力有所提高，0～20 cm土层含水率随复垦年限增加而明显增加，20～40 cm土层含水率增加不明显［5-6］。

国内外逐渐展开探地雷达（GPR）技术在土壤探测方面的应用研究［7-10］，尤其在土壤水探测方面的应用。GPR能够准确测定土壤体积含水率，并反馈含水率时空变化的差异［11］。GPR探测土壤水分有一定的可行性，但其探测精度、分层含水量的确定等方面的研究还需加强［12］。重构土壤的介电常数与含水率间存在明显的函数相关性，应用GPR测定重构土壤的含水率在方法上是可行的［2］。前人的研究主要是验证已有的Topp模型能否用于揭示矿区复垦地土壤体积含水率。笔者以中煤平朔安太堡露天煤矿南排土场为研究对象，通过剖面实测和GPR探测2种方法获取南排土场不同地形条件下重构土壤剖面体积含水率，一是分析不同剖面位置、坡向土壤体积含水率的总体差异和深度差异；二是识别重构土壤体积含水率差异在GPR图像上的表征；三是比较基于GPR的Topp模型和拟合模型（介电常数-土壤体积含水率）模拟土壤体积含水率与实测土壤体积含水率之间的差异；四是讨论不同方法所得土壤体积含水率差异的原因，并尝试建立一种基于矿区典型特征的介电常数-土壤体积含水率关系模型。通过对排土场不同位置及同一位置不同深度处土壤体积含水率的特征进行分析，可为矿区排土造地工艺和土地复垦方案提供依据。

1 研究区概况与研究方法

1.1 研究区概况

研究区位于山西省朔州市平鲁区（图1），是中煤平朔安太堡露天矿早期的一个外排土场，排土时间为1985—1989年，海拔高度为1 360～1 465 m，边坡坡度为20°～40°。排土结束后其岩土容量达1.16×108m3［13］，属于土石混排类型，其中粒径≥50 mm的岩石占46.98%，5～＜50 mm的砾石占15.48%，＜5 mm的砾石、土砂占37.54%［14］。地表覆盖以黄土、红土和红黄土为主。1990年开始进行复垦，采用草-灌-乔的植被配置模式。该研究区所在的原地貌为黄土丘陵，土壤侵蚀严重，同时冬春季节风大风多，地表干燥，是黄土高原典型的生态脆弱区。

1.2 样方与剖面建立

试验在5月份进行，避免了降雨对实验结果的影响。根据高度、坡度和坡向及样方压实情况，在排土场的不同位置建立5个样方（3个平台和2个边坡），样方编号分别为S1、S2、S3、S4和S5。在排土场平台上样方的大小设置为10 m×10 m；在排土场边坡上样方的宽度设置为10 m，样方的长度根据边坡的坡度来确定。在样方内挖重构土壤剖面，剖面编号分别为 P1、P2、P3、P4和 P5。重构土壤剖面的平均深度为100 cm，个别剖面因障碍层太浅，开挖困难，最终开挖深度为60 cm，各样方的基本特征如表1所示。

1.3 GPR探测与参数确定

1.3.1 GPR探测

剖面挖好以后，在GPR探测面上把管径为3和5 cm的钢管分别打进剖面的不同深度处（图2），然后在地表利用GPR进行探测。探测时设定的速度为0.10 m·ns-1，介电常数为8。

1.3.2 介电常数

介电常数的计算公式为

式（1）中，ε为介电常数；c为光速，3.00×108m·s-1；v为电磁波信号在介质中的传播速度，m·ns-1，在该研究中是土层深度与电磁波传播时间的比值。

1.3.3 土壤体积含水率

（1）土壤体积含水率经验值

前人研究确定了GPR土壤体积含水率的经验公式，目前最为常用的是Topp模型［15］。Topp模型反演土壤体积含水率经验公式为

式（2）中，φ为土壤体积含水率，%；ε为介电常数。

（2）土壤体积含水率反演值

根据实测土壤体积含水率与介电常数之间的关系，构建两者之间的拟合关系模型，再根据关系 模型反演土壤体积含水率。

图1 中煤平朔安太堡露天煤矿南排土场Fig.1 South dump site of the Pingshuo Antaibao Open-Pit Coal Mine of China Coal Energy

表1 各个样方及剖面的基本特征Table 1 Basic features of each sample and profile

1.4 土壤样品采集与测试

在剖面内每10 cm采集土壤环刀样，称重后立即带回试验室烘干，在105℃烘箱中烘干至恒重，称取干重。根据湿重和干重确定土壤体积含水率。

2 结果与分析

2.1 土壤体积含水率总体特征

各剖面0～60 cm深度内土壤体积含水率的平均值和标准差如表2所示。各个剖面土壤体积含水率的均值相差较大，土壤体积含水率均值最高的是P5剖面，其值达到18.56%。P5剖面位于S5样方内，该样方位于排土场半阴坡平台，其植被类型以乔木刺槐为主。P2剖面土壤体积含水率均值也高达17.29%，该剖面位于S2样方内，属于阴坡处的边坡样方，其坡度在22°左右，植被类型以乔木刺槐为主。土壤体积含水率均值最低的是P3剖面，其值为4.38%。P3剖面位于S3样方内，该样方位于排土场阳坡，其坡度在33°左右，植被类型以灌木柠条为主。

P1和P4剖面均为排土场平台样方S1和S4内的剖面，分别位于顶部平台与半阳坡平台，土壤体积含水率均值在10%左右。基本上可以看出，阴坡、半阴坡处的土壤体积含水率较高，顶部平台、半阳坡处的土壤体积含水率较低，而阳坡处的土壤体积含水率最低。

图2 样方内剖面深度标定Fig.2 Depth calibration in the soil profile

表2 各剖面不同深度处土壤体积含水率Table 2 Soil volumetric water content at different depths of each profile %

2.2 土壤体积含水率深度差异

从各剖面土壤体积含水率的标准偏差来看（表2），在0～60 cm深度内，P5剖面土壤体积含水率的标准偏差最大，达4.84%，主要是由于0～20 cm深度内土壤体积含水率较高，大于剖面土壤体积含水率均值。P1剖面土壤体积含水率的标准偏差较大，其值达到4.17%，其情况与P5相反，主要是0～20 cm土层土壤体积含水率较小，低于剖面土壤体积含水率均值。P2剖面土壤体积含水率的标准偏差也达3.88%，主要是＞20～50 cm土层土壤体积含水率较小，低于剖面土壤体积含水率均值。P3和P4剖面土壤体积含水率的标准偏差较小，分别为2.22%和1.39%。通过对所有剖面不同深度处的土壤体积含水率的变化分析可知，P3和P4剖面各层土壤体积含水率变化均较小，P1和P5剖面＞20 cm土层的土壤体积含水率变化较小，而P2剖面在＞10～20和＞50～60cm深处呈现2个明显的峰值。

从5个剖面不同深度处土壤体积含水率差异来看（图3），在0～60 cm整体深度内，P1、P4剖面差异不显著，P2、P5剖面差异不显著，但是P1、P4剖面与P2、P5剖面及P3剖面间存在显著性差异（P＜0.05）。从不同的层位来看，在0～20 cm深度内，P3剖面与P2、P5剖面差异显著，但是与P1、P4剖面差异不显著，而P1、P4剖面与P2、P5剖面差异不显著。在＞20～60 cm深度内，P1剖面与 P2、P4、P5剖面差异不显著，而与 P3剖面差异显著。P2、P5剖面差异不显著，而与P3、P4剖面差异显著。P3剖面与P1、P2、P4、P5剖面差异都显著（P＜0.05）。

图3 各剖面不同深度处土壤体积含水率差异Fig.3 Differences in soil volumetric water content atdifferent depths of each profile

2.3 剖面各层土壤介电常数确定

根据标定深度和传播时间确定电磁波在不同土壤层中传播的速度，利用式（1）计算出不同层介质的介电常数（表3）。如图4所示，在P1剖面中2根钢管标定的深度分别是19和43 cm，GPR电磁波传播至标定深度处的时间分别是2.50和7.10 ns，同时根据2处标定深度的时间差确定19～43 cm处电磁波传播的时间为4.60 ns。在剖面P2中存在明显的土壤和煤矸石分层界面，并且分层信息在图像中明显，可以确定的分层深度分别是31、33、43和47 cm；选定31和47cm作为标定深度，GPR电磁波传播至标定深度处的时间分别是6.69和9.06 ns，同时根据2处标定深度的时间差确定31～47 cm处电磁波传播的时间为2.37 ns。在P3剖面中2根钢管标定的深度分别是13和30 cm，GPR电磁波传播至标定深度处的时间分别是1.38和3.85 ns，同时根据2处标定深度的时间差确定13～30 cm处电磁波传播的时间为2.47 ns。在P4剖面中2根钢管标定的深度分别是22和33 cm，GPR电磁波传播至标定深度处的时间分别是3.60和3.90 ns，同时根据2处标定深度的时间差确定22～33 cm处电磁波传播的时间为2.30 ns。在P5剖面中2根钢管标定的深度分别是20和35 cm，GPR电磁波传播至标定深度处的时间分别是4.37和7.32 ns，同时根据2处标定深度的时间差确定20～35 cm处电磁波传播的时间为2.95 ns。

表3 各剖面不同深度处GPR的探测参数值Table 3 Detection parameter values of GPR at different depths of each profile

图4 基于GPR图像的各剖面深度标定Fig 4 Depth calibration of each profile based on GPR image

2.4 土壤体积含水率经验公式计算

对比剖面分层土壤体积含水率的实测值与GPR土壤体积含水率Topp模型计算出的分层土壤体积含水率，分析各重构土壤剖面不同分层处土壤体积含水率的差异。由于采样实测的值是每10 cm分层采样测试所得，因此实测值均为10 cm深度的整数倍，而剖面标定的深度非10 cm深度的整数倍（表4）。为了2种测量方式的深度对应，假定标定深度内的介质与其对应的10 cm整数倍深度内的介质相对均一。通过对比发现，经验公式计算出的土壤体积含水率略高于采样实测出的土壤体积含水率；2种方式测量出的土壤体积含水率差值的绝对值最大为3.20%，最小为0.13%；实测值与经验值的平均偏差率为13.42%。

表4 各剖面不同深度处采样实测与GPR反演土壤体积含水率对比Table 4 Comparison of sampling and GPR inversion of soil volumetric water content at different depths of each profile

2.5 土壤体积含水率介电常数反演

把剖面各层采样的土壤体积含水率实测数据和对应层的GPR介电常数进行拟合，并构建拟合关系模型。在进行拟合分析之前，假定采样深度处的土壤体积含水率与标定深度处的体积含水率接近。一共有15组土壤体积含水率实测值与介电常数值，根据对多个拟合结果的筛选与比对，最后选择13组土壤体积含水率实测值与GPR介电常数进行拟合分析，拟合结果如图5所示。其拟合关系的R2值高达0.950 5，拟合效果好。此外，土壤体积含水率与介电常数呈三次函数关系式，与GPR土壤体积含水率的经验公式表现出同样的关系趋势，并且回归阶数均为3。

将15个介电常数值代入图5中的拟合关系模型进行土壤体积含水率反演，其反演结果如表5所示。通过对比发现，模型反演出的土壤体积含水率要略低于采样实测出的土壤体积含水率；2种方式测量出的土壤体积含水率差值的绝对值最大为3.11%，最小为0.23%；实测值与反演值的平均偏差率为9.83%。

图5 土壤体积含水率实测值与介电常数的拟合关系Fig.5 Fitting relationship between measured values of soil volumetric water content and dielectric constant

3 讨论

3.1 不同剖面土壤体积含水率差异的原因

矿区排土场不同位置及同一剖面不同深度处重构土壤剖面体积含水率差异明显（表5）。即土壤体积含水率呈阴坡/半阴坡＞顶部平台＞半阳坡/阳坡的规律性变化。山西省地处中纬度大陆性季风气候区，干燥期大于雨雪期，年光照时间较长。阳坡土壤水分蒸发量大于阴坡和平台，土壤体积含水率小于平台和阴坡。这与曹银贵等［16］对黄土矿区复垦地土壤物理性质分异与机理以及张建军［17］对山西省固县林场野生连翘不同坡向生长状况的研究结论相同。

表5 各剖面不同深度处采样实测与Topp模型反演土壤体积含水率对比Table 5 Comparison of sample volume measurement and model inversion soil volumetric water content at different depths of each profile

P2剖面和P5剖面是典型的土石混排剖面，剖面0～20 cm以表土为主，＞20 cm是土石混排，表土土壤孔隙大，储水能力强［18］。剖面又属于阴坡，日光照时间少，所以0～20 cm土壤体积含水率较高。P3剖面是边坡又是阳坡，植被类型是柠条，且保水的草本几乎没有，所以其土壤体积含水率较小。P1和P4剖面都在平台，两者土壤体积含水率差异不大，且沿深度方向土壤体积含水率差异也不大。从方差分析的结果来看，因为P1和P4剖面位于平台，P2和P5剖面位于阴坡，P3剖面位于阳坡边坡，区位的不同使 P1、P4剖面与 P2、P5剖面及 P3剖面都存在显著性差异。从不同层位来看，不同深度处土壤体积含水率差异性显著，特别在0～20 cm处最显著。表土层所处环境较为复杂，其土壤体积含水率对蒸发和地形等自然条件反应更为敏感［19］。

3.2 土壤体积含水率特征的GPR图像呈现

如图6所示，对土壤体积含水率差异较大剖面，GPR图像经过相同处理后可观察到土壤体积含水率大的剖面有更多的低频信号。含水率较高的剖面或土层介电常数相对较大，高频信号衰减较快，信号以低频为主。P3剖面和P4剖面沿深度方向土壤体积含水率大小接近，信号图上显示信号均一。P5剖面土壤体积含水率均值远大于P3剖面，探地雷达信号图上P5剖面比P3剖面有更多的低频信号。同一剖面不同深度的剖面含水率较高的地层具有更多的低频信号，如0～20 cm处信号图上部明显具有较多低频信号。P2剖面在＞10～20和＞50～60 cm深度处土壤体积含水率呈现2个明显的峰值，在信号图上也有较好的呈现。对与土壤体积含水率差异较小的剖面或土层，因为GPR分辨率及数据处理技术的局限，在图像上难以直观呈现，国外学者在这方面研究较多［20-22］。

3.3 不同方法所得土壤体积含水率差异的原因

Topp模型和拟合关系模型反演所得土壤体积含水率与实测值的相关系数都在0.90以上，但2种方法所得土壤体积含水率有一定的偏差。Topp模型是土壤体积含水率与介质介电常数的关系模型，目前已被广泛应用于反演土壤体积含水率［2，23-24］。探地雷达区分介质的前提是目标体与周围介质具有一定的电性差异，即介电常数有一定的变化。但若土壤中不只有水这种介电常数大的介质存在，或土壤体积含水率较小情况下由Topp公式反演的土壤体积含水率就会不准确［15］。当使用Topp模型所得土壤体积含水率与实测土壤体积含水率存在较大误差时，需要寻找更合适的模型，或对Topp模型进行修正。Topp模型不可能适用于所有土壤体积含水率的测定［1，8］。TOMER等［25］认为Topp公式适用于沙质土壤，不适用于粗糙的密度较大的碎屑岩石。矿区排土场多属土石混排类型，且有煤矸石层存在，Topp模型可能不适用于此类研究区。实际土壤探测过程中如果能对土壤中的介质类型有一个充分的了解，更有利于Topp模型的使用与修正，或选取更合适的模型。

图6 探地雷达信号图定性分析土壤体积含水率Fig.6 Qualitative analysis of soil water content by ground penetrating radar signal diagram

由表3可知，P1剖面＞20～40 cm、P2剖面＞30～50 cm深度模型反演和采样实测的偏差率较大。由式（1）～（2）及速度、时间、和深度关系可知，土壤体积含水率大小取决于介电常数，介电常数大小取决于电磁波波速度，电磁波速度决定于标定深度和时间的准确性。在标定深度过程中，如果尺子实测钢管埋藏深度与实际深度存在一定差异，所求电磁波在介质中的传播速度就会存在较大偏差，介电常数就存在偏差，导致实测和反演所得土壤体积含水率存在较大偏差。信号处理过程中若钢管埋藏深度处时间拾取不准确也会影响电磁波波速大小，使实测和反演所得土壤体积含水率存在偏差。

4 结论

通过剖面实测和探地雷达探测2种方式，得到了安太堡露天矿区复垦排土场土壤体积含水率，对不同位置及同一位置不同深度处土壤体积含水率特征进行分析，同时对基于GPR获取土壤体积含水率的不同模型进行了探讨。通过分析得到以下结论：

（1）矿区排土场不同位置及同一剖面不同深度处土壤体积含水率大小存在差异性。

（2）Topp模型的应用有一定的适用范围，实际操作中根据土壤类型选取合适模型可提高GPR技术在计算土壤体积含水率时的准确性。

（3）介电常数与土壤体积含水率拟合精度较高，通过探地雷达探测土壤体积含水率的方法是可行的。