0 引 言
陕北煤田作为我国煤炭资源储量最大的煤田之一，因其可采煤层多、煤层厚、煤质优良、构造简单、开采条件优越而受广泛关注[1-2]。但随着开采强度、开采规模的日益增大，煤炭开采对生态环境的影响也日益突出[3]。为了保护和恢复当地生态环境，有必要对该地区的土壤进行研究。土壤含水率是土壤理化性质中影响植物生长、光合速率、蒸腾作用和气孔导度等的重要条件之一[4]。研究区上的植被主要依靠浅部砂壤中赋存的水资源生存，因此准确表征浅地表土壤尤其是砂壤土的含水量对指导生态环境的恢复至关重要[5]。
目前测量土壤含水率主要有烘干法、中子法、时域反射仪(Time Domain Reflectometry，TDR)法、遥感和探地雷达法等。其中烘干法、中子法、TDR法为点测量，测量精度高，但仅能针对某些点位进行测量且对土壤具有破坏性，无法进行大面积监测[6-7]；遥感技术由于其监测平台的独特性可进行区域性的土壤水分测量工作，并能够对同一地区进行长时间持续观测，但监测深度浅，仅能获得深度10 cm内土壤水分值[8]；相对于以上方法，探地雷达具有快速测定、对土壤无损伤、相对探测深度深等优点，且覆盖面积可达几公里范围[9]。因此其作为适用于中尺度范围土壤含水率的测量方法被逐渐应用[10]。SHAMIR等[11]通过探地雷达研究了以色列典型土壤中1.5 m深度的空间土壤含水率；ALLRED等[12]以探地雷达与TDR结合的方式绘制了整个高尔夫球场果岭砂层含水量；MINET等[13]使用探地雷达获得了2种介质下连续土壤水分剖面；崔凡等[14]基于探地雷达技术针对中国西部砂壤土提出了通过功率谱计算含水率的测量方法。总体上探地雷达测量含水率使用的计算方法有2种：①基于雷达波的频谱分析进行估算；②基于TOPP等[15]在1980年提出的含水率与介电常数之间的经验关系进行计算。由于前者存在深部信号衰减严重导致测量值出现较大误差的问题，因此测量较深部含水率基本采用介电常数法。但实践证明，Topp公式计算不同类型土壤含水率的准确性不一，ROTH等[16]对公式进行了修正，但并未明显提高针对砂质土壤的适用性。文献[17]提出的半理论半经验公式作为同样应用较多的体积含水率计算公式，经JACOBSEN等[18]的验证，可以获取较高的体积含水量测定精度，但针对不同土壤必须首先准确确定其校正参数，因此可由这2种公式进行相互验证。
针对我国西部风积沙区域砂壤介质，通过探地雷达数据与钻孔土样实测含水率值进行综合分析，用钻探方法标定地层界面位置，获取雷达波精确波速并计算出介电常数，再结合取样烘干法测实际含水率，对Topp公式形式的“介电常数-体积含水率值”模型进行了改进，使其更适用于中国西部，尤其是陕北地区风积沙区域砂壤介质的中砂与细砂含水率计算，从而更好地反映地下煤矿开采条件下地表浅层介质含水率分布情况，对保护和恢复当地生态环境具有一定的指导意义。
1 研究区概况
研究区所在的柠条塔矿区南翼位于陕北煤田东北部，属毛乌素沙漠西南缘。柠条塔矿区煤层厚，埋藏浅，地质构造简单，符合中国西部煤矿典型开采条件。矿井采用斜井开拓，综合机械化采煤，全部垮落式管理顶板。根据钻孔揭露，研究区内浅层基本被风积沙与萨拉乌苏组沙层所覆盖，风积层厚度为0～15.6 m，平均厚度为7.7 m，矿物成分主要为石英、长石，黏土矿物较少，粒径以0.05～2.00 mm中细砂为主。全区干旱少雨，蒸发量大，多年平均蒸发量为1 712 mm。
2 含水率计算研究方法
2.1 探地雷达探测原理
探地雷达向地下发射脉冲形式的高频宽带电磁波，电磁波在地下介质传播过程中，当遇到存在电性差异的目标体时便会发生反射，并由接收天线所接收。对接收到的电磁波信号进行处理和分析，根据信号波形、强度等参数来推断目标体的结构、几何形态等，从而达到对地下目标物的探测[19]，探地雷达探测原理如图1所示。
图1 探地雷达探测原理
Fig.1 Detection principle of ground penetrating radar
电磁波在地下传播过程中，介质的相对介电常数与电磁波波速满足以下关系：
(1)
式中：t为电磁波在介质中传播的旅行时间，ns；z为电磁波所遇到的反射面深度，m；v为电磁波在当前介质中的传播速度，m/ns；c为电磁波在真空中的传播速度，m/ns；εr为介质的相对介电常数。
由式(1)可知，通过求取电磁波在当前介质的波速可以获得该介质的介电常数，进而通过Topp公式计算含水率。Topp公式为
(2)
其中：θV为计算体积含水率，%。由此可知准确获取电磁波在当前介质中的传播速度是求取介质含水率的关键。
2.2 实测数据采集
1)测线布置。为避免煤层开采过程对试验结果产生的影响，根据文献[20-21]研究结果(地裂缝超前于工作面开采位置平均10.359 m，而裂缝对于周围土壤含水率影响范围为75 cm以内)，选取工作面中部未受地下采煤扰动影响区域作为研究区，该区域原始浅层砂壤介质中、细砂作为研究对象，并垂直于工作面推进方向布设2条测线，测线1、2分别距切巷2 136 m和1 656 m，2条测线长度均为1 970 m，进行数据采集时工作面已推进至距开切眼610 m处。测线布设位置如图2所示。
图2 研究区测线布置
Fig.2 Survey line layout in study area
2)雷达探测。试验探地雷达采用中国矿业大学(北京)自主研发的GR雷达系统，天线中心频率为100 MHz，采样点数为1 024，时窗为200 ns。使用探地雷达分别沿2条测线进行探测，如图3a所示，以获取探地雷达实测数据。
3)实地取样。同时为对探地雷达获得的含水率进行校准和验证，每条测线每隔200 m设置1个取样点，取样点分布情况如图3所示。采用浅部地层钻探方式取样。以洛阳铲为工具，对每条测线上的测点进行打钻取样，每隔0.5 m进行1次提钻，对钻内砂壤使用环刀取样，为保证砂壤物性不被破坏，土样取出后立即放入铝盒内密封保存。取样工作如图3b所示。
图3 雷达探测与洛阳铲取样现场施工
Fig.3 Site construction photos of radar detection and shallow sampling
4)测定砂壤样品体积含水率。利用电子秤对铝盒内的土样及铝盒进行称重，得到质量m0，之后将其置于105 ℃的烘箱内烘12 h左右直至恒重，称其质量m1。所用环刀规格为ø50.46 mm×50 mm，容积V0为100 cm3，则砂壤体积含水率计算式为
(3)
式中：
为实测体积含水率，%；m0为湿土和铝盒质量，g；m1为干土和铝盒质量,g；ρw为水的密度g/cm3。
2.3 GPR波速标定方法
通过探地雷达数据识别介质反射面，获得雷达剖面中雷达波到达反射界面双程走时，同时结合钻探取样结果对研究区不同介质位置进行标定，如图4所示。在已知介质层位深度与雷达波走时条件下，精确求得其在不同介质中的波速，进而得出不同介质的相对介电常数。
图4 雷达波剖面
Fig.4 Radar wave reflecting surface
3 结果与讨论
3.1 计算含水率与实测含水率误差分析
由于中砂或细砂在部分钻孔中并未出现，因此由GPR波速标定方法获取细砂介质介电常数10组，中砂介质介电常数19组，见表1、表2。通过经典Topp公式对2种砂壤介质进行计算获得体积含水率(θV)，并与这2种介质的实测含水率
进行对比，结果如图5所示。其中，误差率为误差值与实测值的比值。
由表1、表2及图5可以看出，2种粒径砂壤实测含水率平均值与Topp公式计算含水率平均值均存在较大误差。
表1 细砂Topp公式计算含水率与实测含水率误差
Table 1 Errors of fine sand water content between calculated value by Topp formula and actual value measured
取样点εr计算含水率/%实测含水率/%误差L1-16.100.105 6300.034 0660.071 564L1-228.300.153 9290.086 2320.067 696L1-328.200.151 8280.085 1950.066 633L1-427.200.130 3320.064 7310.065 601L1-527.350.133 6150.065 0560.068 558L1-626.350.111 3430.035 7940.075 548L1-727.550.113 6120.062 1700.051 442L1-927.300.137 9590.065 0530.072 905L2-326.050.132 5230.062 2510.070 271L2-928.400.156 0200.101 4560.054 564
表2 中砂Topp公式计算含水率与实测含水率误差
Table 2 Errors of medium sand water content between calculated value by Topp formula and actual value measured
取样点εr计算含水率/%实测含水率/%误差L1-15.650.095 1980.045 8550.049 343L1-125.600.094 0270.048 0640.045 963L1-228.200.151 8280.097 8740.053 954L1-326.800.121 4800.056 3660.065 113L1-426.800.121 4800.047 0370.074 442L1-527.450.135 7910.066 5580.069 232L1-626.450.113 6120.053 3690.060 243L1-726.000.103 3280.043 4790.059 849L1-826.300.110 2050.034 9190.075 286L1-927.500.136 8760.070 8590.066 017L2-28.000.107 9220.046 6010.061 321L2-127.110.147 6010.0692670.078 334L2-227.550.149 7190.088 8290.060 890L2-427.650.105 6300.052 7750.052 855L2-528.000.147 6010.077 3840.070 217L2-626.100.140 1170.074 6520.065465L2-728.100.137 9590.075 1890.062 770L2-828.000.128 3530.064 5310.063 822L2-926.200.147 6010.095 5880.052 013
图5 Topp公式计算值与实测值对比
Fig.5 Comparison of actual value measured and value canculated by Topp formula
其中，细砂误差范围为0.051 4～0.075 5，平均误差为0.066 5，与实测值相比误差率范围为53.78%～211.06%，平均误差率达114.61%。中砂误差范围为0.046 0～0.078 3，平均误差为0.062 5，与实测值相比误差率范围为54.41%～215.60%，平均误差率达106.53%。其精度显然无法满足计算该类地区土壤含水率的要求，因此有必要对Topp公式进行修正，建立更为精确的砂壤含水率计算模型。
3.2 Topp公式修正模型建立与验证
为获得更精确的且适用于西部不同粒径砂壤的体积含水率与介电常数关系模型，选1号测线测试的细砂和中砂数据进行多项式回归分析，结果如图6所示。
图6 1号测线细砂和中砂实测数据回归
Fig.6 Regression of measured data of fine sand and medium sand in line 1
由拟合结果可得相对介电常数与体积含水率关系如下
0-0.963 5εr+0.134 3εr2-0.006 1εr3
(4)
式中：
为由Topp公式修正后的计算体积含水率。为了验证式(4)的准确度，使用该公式对所有2号测线数据含水率重新进行计算并与实测数据进行对比，结果见表3、表4。
表3 细砂拟合公式计算平均值与实测含水率平均值误差
Table 3 Errors of average values of fine sand water content between the results obtained from fitting formula and actual measurements
取样点含水率/%拟合公式计算实测误差L2-320.027 7250.062 2510.034 526L2-920.059 3140.101 4560.042 142
表4 中砂拟合公式计算平均值与实测含水率平均值误差
Table 4 Errors of average values of medium sand water content between the results obtained from fitting formula and actual measurements
取样点含水率/%拟合公式计算实测误差L2-20.056 0000.046 6010.009 400L2-120.038 1670.069 2670.031 100L2-220.047 7610.088 8290.041 068L2-420.049 8440.052 7750.002 931L2-520.056 0000.077 3840.021 384L2-620.027 3690.074 6520.047 283L2-720.057 2830.075 1890.017 906L2-820.056 0000.064 5310.008 531L2-920.026 9910.095 5880.068 5968
由表3和表4可以看出，由Topp公式修正后的计算的体积含水率，细砂误差范围为0.034 5～0.042 1，平均误差为0.038 3，与实测值相比误差率范围为41.54% ～55.46%，平均误差率为48.50%。中砂误差范围为0.002 9～0.068 6，平均误差为0.027 6，与实测值相比误差率范围为5.55%～71.76%，平均误差率为35.18%。与经典Topp公式相比，由修正后Topp公式计算的体积含水率，细砂误差范围减小了68%，平均误差下降了42%，误差率范围减小了91%，平均误差率下降了近58%；中砂误差范围虽然有所扩大，但平均误差下降了近56%，误差率范围减小了近59%，平均误差率下降了近67%。
4 结 论
1) 针对我国西部浅层砂壤介质的中砂与细砂，对探地雷达法测量砂壤含水率所用经典Topp公式进行修正，建立了Topp公式的修正模型。结果表明2种粒径砂壤体积含水率准确度均大幅提高。其中，修正后Topp公式计算的体积含水率，细砂误差范围减小了68%，平均误差下降了42%，误差率范围减小了91%，平均误差率下降了近58%；中砂误差范围虽然有所扩大，但平均误差下降近56%，误差率范围减小近59%，平均误差率下降近67%。
2)Topp公式修正模型的准确性，验证了GPR波速标定方法对中国西部砂壤介质含水率的精确计算方法的合理性，对计算中国西部浅层砂壤含水率具有应用价值。
3) 基于GPR波速标定方法的Topp公式的修正模型可应用于中国西部陕北柠条塔的砂土介质体积含水率计算。由于目标介质的物性参数存在较大差异，本模型的广泛适用性需进一步验证。探地雷达数据与钻孔土样实测含水率的综合分析方法并不受限于研究区介质种类及物性，应用前景良好，可为更多地区生态环境的改良和恢复提供理论依据。
[1] 申 涛,朱占荣.陕北矿区煤炭开采沉陷实测参数分析[J].煤炭科学技术,2019,47(12):207-213.
SHEN Tao,ZHU Zhanrong.Analysis of measured mining subsidence parameters in Northern Shaanxi Mining Area[J].Coal Science and Technology,2019,47(12):207-213.
[2] 侯恩科,童仁剑,王苏健,等.陕北侏罗纪煤田风化基岩富水性Fisher模型预测方法[J].煤炭学报,2016,41(9):2312-2318.
HOU Enke,TONG Renjian,WANG Sujian,et al.Prediction method for the water enrichment of weathered bedrock based on Fisher model in Northern Shaaxi Jurassic coalfield[J].Journal of China Coal Society,2016,41(9):2312-2318.
[3] 范立民,马雄德,蒋泽泉,等.保水采煤研究30年回顾与展望[J].煤炭科学技术,2019,47(7):1-30.
FAN Limin,MA Xiongde,JIANG Zequan,et al.Review and thirty years prospect of research on water-preserved coal mining[J].Coal Science and Technology,2019,47(7):1-30.
[4] TIAN Xiaoming,WANG Kaiyong,LIU Yinghui,et al.Effects of polymer materials on soil physicochemical properties and bacterial community structure under drip irrigation[J].Applied Soil Ecology,2020,150:103456.
[5] AINIWAER Mireguli,DING Jianli,KASIM Nijat,et al.Regional scale soil moisture content estimation based on multi-source remote sensing parameters[J].International Journal of Remote Sensing,2020,41(9):3346-3367.
[6] 卢奕竹.基于探地雷达的土壤水遥感监测地面验证及尺度问题研究[D].北京：中国水利水电科学研究院,2017.
[7] 易小波.西北干旱区土壤含水量时空变化特征及土壤物理性质模拟试验研究[D].杨凌:西北农林科技大学,2017.
[8] 姜 红,玉素甫江·如素力,拜合提尼沙·阿不都克日木,等.基于支持向量机回归算法的土壤水分光学与微波遥感协同反演[J].地理与地理信息科学,2017,33(6):30-36.
JIANG Hong,YUSUFUJIANG Rusuli,BAIHETINISHA Abudoukerimu,et al.Soil moisture retrieval by synergizing optical and microwave remote sensing data based on support vector machine regression algorithm[J].Geography and Geo-Information Science,2017,33(6):30-36.
[9] KLOTZSCHE Anja,JONARD François,LOOMS Majken Caroline,et al. Measuring soil water content with ground penetrating radar:A decade of progress[J].Vadose Zone Journal,2018,17(1):541-549.
[10] 聂俊丽.基于地质雷达技术的采煤对浅部地层含水量影响规律研究[D].北京：中国矿业大学(北京),2014.
[11] SHAMIR O,GOLDSHLEGER N,BASSON U,et al.Mapping spatial moisture content of unsaturated agricultural soils with ground-penetrating radar [J].International Archives of the Photogrammetry,Remote Sensing & Spatial Information Sciences,2016,41(XL1-B8):1279-1285.
[12] ALLRED Barry,FREELAND Robert,GROTE Katherine,et al.Ground-penetrating radar water content mapping of golf course green sand layers[J].Journal of Environmental and Engineering Geophysics,2016,21(4):215-229.
[13] MINET Julien,WAHYUDI Agung,BOGAERT Patrick,et al.Mapping shallow soil moisture profiles at the field scale using full-waveform inversion of ground penetrating radar data[J].Geoderma,2011,161(3-4):225-237.
[14] 崔 凡,陈柏平,吴志远,等.基于探地雷达功率谱和雷达波振幅包络估算砂壤含水率[J].农业工程学报,2018,34(7):121-127.
CUI Fan,CHEN Baiping,WU Zhiyuan,et al. Soil moisture estimation based on GPR power spectrum and envelope amplitude in sand loam[J].Transactions of the Chinese Society of Agricultural Engineering,2018,34(7):121-127.
[15] TOPP G C,DAVIS J L,ANNAN A P.Electromagnetic determination of soil water content:measurements in coaxial transmission lines[J].Water resources research,1980,16(3):574-582.
[16] ROTH C H,MALICKI M A,PLAGGE R.Empirical evaluation of the relationship between soil dielectric constant and volumetric water content as the basis for calibrating soil moisture measurements by TDR[J].Journal of Soil Science,1992,43(1):1-13.
[17] HERKELRATH W N,HAMBURG S P,Murphy Fred.Automatic,real‐time monitoring of soil moisture in a remote field area with time domain reflectometry[J].Water Resources Research,1991,27(5):857-864.
[18] JACOBSEN Ole,SCHJONNING P.Comparison of TDR calibration functions for soil water determination[C]//:Time domain reflectomety applications in soil science.Denmark:Danish Institute of plant and Soil Science,1995:25-33.
[19] 杨 峰,彭苏萍.地质雷达探测原理与方法研究[M].北京：科学出版社,2010.
[20] 胡振琪,王新静,贺安民.风积沙区采煤沉陷地裂缝分布特征与发生发育规律[J].煤炭学报,2014,39(1):11-18.
HU Zhenqi,WANG Xinjing,HE Anmin.Distribution characteristic and development rules of ground fissures due to coal mining in windy and sandy region[J].Journal of China Coal Society,2014,39(1):11-18.
[21] 台晓丽,胡振琪,陈 超.风沙区采煤沉陷裂缝对表层土壤含水量的影响[J].中国煤炭,2016,42(8):113-117.
TAI Xiaoli,HU Zhenqi,CHEN Chao.Effect of mining subsidence fissure on moisture of surface soil in aeolian sand area[J].China Coal,2016,42(8):113-117.
