基于探地雷达的煤矿开采区地表土壤含水率变化研究

吴志远，尹尚先，马丽红

(华北科技学院 安全工程学院，北京 东燕郊 065201)

0 引言

地表土壤水分在控制植被生长、调节地表生态环境等方面起着重要的作用，特别是在干旱地区，地表土壤水分是防止土壤沙漠化的关键因素之一[1-3]。地处毛乌素沙地边缘的神府东胜矿区，是中国目前已探明的煤炭储量最丰富的地区，该地区煤矿的开采对保证国家优质动力煤供应及促进国民经济发展具有重要意义。但该区域生态环境脆弱，蒸发强烈、降雨量较低、地表砂质土壤疏松，地表土壤表现出水分少、水分存储能力差的特点[4-7]。又由于地下煤层开采，打破了本来已十分脆弱的生态环境，致使开采区发生塌陷，导致本来已十分脆弱的生态环境进一步恶化[8-10]。因此了解煤矿开采前后地表土壤水分含量及其分布变化的主导因素是进行开采区环境治理及生态恢复的关键环节。

探地雷达(Ground penetrating radar，GPR)已被证明能够较好地获得地表土壤水分含量的变化情况[11-12]。其具有探测面积大、探测速度快、无损探测等优点，避免了传统方法如中子仪、时域反射仪法(Time-Domain Reflectometry，TDR) 等只能点测、费时费力的缺点，且没有辐射的危害。与遥感探测方法相比，探地雷达具有探测精度高、受植被干扰小等优点[13-18]，非常适合土壤水分的探测。

1 材料与方法

1.1 研究区概况

神东矿区地处中国西北部，矿区属半干旱大陆性季风气候，干燥少雨，大风频繁，多年平均降水量为194.7～531.6 mm，如图1所示，多年年平均蒸发量为2 297.4～2 838.7 mm。由于气候干旱、地表径流较少，导致该地区地下潜水面深度一般在30 m以下。同时由于地下煤层开采导致地表出现大量裂缝，进一步加剧了该地区生态的恶化情况。

图1 研究区2013、2014年8月份降雨量分图

所选研究区位于神东煤矿大柳塔矿区，如图3所示。研究区主要属风沙堆积地貌，上部由砂层覆盖，砂层厚度0～20 m不等。区内高程最大变化值为35.2 m，地势中间高两边低，植被主要为沙柳，由于前几年煤矿公司的绿化，研究区内分布着少量的杨树。

1.2 仪器及探测方法

本次探测利用GR探地雷达系统对研究区进行数据采集。GR探地雷达系统主要由一个主机和两个天线盒子组成，两个天线盒子一个为发射天线，一个为接收天线，如图2所示。天线的中心频率为400 MHz天线，雷达参数采样点数为2048，空间分辨率0.2 m。

图2 地表直达波方法探地雷达工作原理图

在开采区中心位置选取一条900 m长的L1测线进行雷达探测及取样测试；同时在未开采区内选取了一条900 m长测线L2进行雷达探测及取样以判断两次探测过程中降雨及季节变化等因素是否对研究区地表土壤含水率产生影响。各测线与正北方向呈30°夹角，如图3所示；大范围的土壤含水率变化由小范围内土壤含水率变化引起，为了了解小面积范围内土壤含水率的变化规律及其影响因素，在开采区选取了 40 m×40 m 正方形区域进行了雷达探测，如图3中红色方块所示。

图3 研究区地理位置图，红色星号为研究区位置

2013年8月30日用探地雷达对两条900 m测线进行了探测，以确定收发天线的最优天线距及地下煤层开采前地表土壤的水分含量及分布情况。同时对测线进行钻孔取样，取样间隔为25 m，取样工具为高10 cm，体积为100 cm3的环刀，利用烘干法在室内对样品含水率进行测试,以对比分析雷达探测含水率的精确性。然后在开采区较平坦位置选出一个40 m×40 m的正方形区域进行雷达探测及取样。正方形区域内的探测在直角坐标系内进行，在X、Y方向每隔4 m布置一条测线。由于受沙柳、杨树等植被的影响，会在探测时对测线位置做较小的改动。2014年8月30日探地雷达探测及取样与2013年方式相同。

探地雷达探测在钻孔取样之前，但时间相差很短，以保证含水率探测的准确性。在当天把样品称重后放入实验室KLM29-CW-01烘箱内，在105°C的条件下烘8～12个小时直至恒重，在烘干后取出土样并计算烘干后土样干质量，然后计算出水分含量，最后计算出土壤体积含水率。

1.3 地表直达波方法探测土壤含水率原理

电磁波在地下传播时对土壤介电常数的变化非常敏感，含水量较小的变化将会引起介质的介电常数改变， 进而引起电磁波在土壤中的传播速度。地表直达波在地表的最上部传播, 它的振幅随距离衰减很快， Du 和 Rummel认为在没有明确地下反射层的条件下， 地表直达波是测定浅表层土壤含水率的最佳方法[19-21]。在GPR 多偏移距剖面上很容易识别出地表直达波。地表直达波走时与天线距离之间呈线性关系，如图4所示，直线的斜率即是地表直达波速度。

图4 探地雷达多偏移距方法示意图

但是在实际操作过程中GPR 多偏移距测量会花费大量时间和人力，不适合应用于大范围探测。因此，Du 和Sperl 提出用以下3 个步骤确定土壤水含率 ：①用多偏移距方法测量, 确定采用不同天线距离时所对应的地面波大致到达时间。②选择一个天线距离, 使地面波与空气波和反射波明显地区分开。③使用这个天线距离做共偏移距测量, 将地面波旅行时间的变化与土壤介电常数的变化联系起来。

地面波走时tGW( s) 、天线距离x( m) 和土壤介电常数之间最直接的关系由下式给出:

(1)

式中，tAW为空气波走时，s。

Topp等人在1980年提出了如今已较为常用的“εθ”模型[22]，其公式如下：

θ=-5.3×10-2+2.92×10-2ε-5.5×10-4ε2+4.3×10-6ε3

(2)

其中θ为体积含水率，ε为相对介电常数。

2 最优天线距确定

如图5所示，L2测线在横向上穿越了9个沙丘(分别位于25、142、250、300、400、525、650、775、875处)，测线内最大高差为19.9 m。2013年和2014年8月30日分别利用多种探地雷达方法沿L2测线进行了浅部地层土壤含水率探测，同时每隔25 m进行一次取样，取样点共37个。遵循天线分离距内雷达波信号干扰最小、子波最清晰的原则，2013年和2014年从多个不同天线距的探地雷达测量效果中找出了两个最优天线分离距，分别为102 cm和114 cm。

3 煤矿开采对地表含水率变化影响

3.1 未开采区2013年与2014年地表土壤含水率对比

图5给出了L2测线2013、2014年探地雷达及取样实测30 cm以上砂层含水率的粗略分布图。从图中可以看出，探地雷达反演含水率比取样实测含水率要高，2013、2014年雷达探测所得含水率平均值分别为0.0837 cm3/cm3、0.0802 cm3/cm3比2013/2014年取样实测含水率0.0636 cm3/cm3、0.0548 cm3/cm3分别大0.0201 cm3/cm3和0.0254 cm3/cm3。造成这种情况的原因有两个：一是于L2测线地表地形起伏较大，雷达在上下坡时贴地不紧也会致使探测结果存在误差；二是两种方法探测的尺度不同。同时地表土壤含水率越小，探地雷达与实测含水率的误差也就越大。

图5 L2测线探地雷达和钻孔取样实测含水率对比图。右侧为使用归一化高程处理后的测线内各点的相对高度。钻孔取样深度为30 cm。

2013年L2测线雷达探测土壤含水率与实测土壤含水率的相关性分析图如图6所示，从图中可以看出探地雷达探测含水率比实测含水率略大，两种方法所探地表土壤含水率的相关系数为0.76，具有较强的相关性。表明利用探地雷达能够获得与实测相近的土壤含水率。

从探地雷达两次探测结果中发现，L2测线内地表土壤含水率变化与沙丘分布关系密切。沙丘脊部地表土壤体积含水率较低(0.05～0.07 cm3/cm3)，其中顶部水分最小；沙丘底部体积地表土壤含水率较高(0.07～0.012 cm3/cm3)，沙丘迎风坡地表土壤含水率略高于沙丘背风坡土，这与吕贻忠等[23]的研究结果相同。同时可以看出在测线最西北端海拔较低处，含水率比其他地区略高。

图6 L2测线2013年探地雷达探测含水率与实测含水率相关性图

两次探地雷达反演平均体积含水率基本相同，2013年(0.0837 cm3/cm3)比2014年(0.0802 cm3/cm3)稍高，这说明两次探测过程中未开采区地表土壤水分基本没有变化，表明两次探测时间内地表土壤水分基本不受降雨等因素的影响。

3.2 开采区2013年和2014年地表土壤含水率对比

3.2.1 L1测线2013年和2014年地表土壤含水率对比

L1测线位于开采区中间位置，测线内主要有四沙丘(400、550、700和875 m)，主要位于测线西北端，测线内最大高差为22 m。图7为L1测线2013、2014年探地雷达及取样实测砂层含水率的粗略分布图。

使用GPS对L1测线高程测量发现，开采后地表沉降值在1.27～2.92 m之间，平均下降了1.95 m。从图7中可以看出，2014年探测含水率(平均为0.099 cm3/cm3)比2013年探测含水率(平均为0.083 cm3/cm3)大，推测这主要是由于近年来矿区内植树造林增强了地表土壤的蓄水能力。L1测线内沙丘顶部地表土壤水分较低，如400 m、550 m处体积含水率分别为0.063 cm3/cm3、0.065 cm3/cm3；沙丘底部地表土壤体积含水率较大，如475 m、625 m处，地表土壤含水率分别为0.131 cm3/cm3、0.135 cm3/cm3，沙丘迎风坡(375 m、525 m处)土壤水分相比背风坡(410、600 m)含量略高。这与未开采区的L2测线地表土壤含水率的变化规律基本相同。

图7 L1测线探地雷达含水率分布图。右侧为使用归一化高程处理后的测线内各点的相对高度。钻孔取样深度为30 cm

3.2.2 正方形区域2013年和2014年地表土壤含水率对比

大范围内土壤含水率变化主要由小范围内土壤含水率变化引起。本次研究在开采区内靠近沉降中心位置较平坦地区选择了一个40 m×40 m的正方形区域以研究小范围内开采区土壤含水率的变化规律。正方形区域内植被主要为沙柳和杨树。图8(a)、(b)中的红点为杨树根茎位置，绿色沙柳根茎位置(沙柳较密集)，黑色圆圈为树冠覆盖直径，其余部分为裸砂地。图8(a)、(b)展示了2013年8月26和2014年8月26日各种环境因子条件下正方形区域浅部砂层体积含水率分布图。

图8 利用克里格插值法所得 2013、2014年两次探测含水率平面分布图。注：虚线为含水率为6 cm3/cm3等值线，实线为含水率为8 cm3/cm3等值线，红色实线为归一化等高线。

从图8(a)、(b)中可以看出2014年探地雷达所得体积含水量(平均为0.096 cm3/cm3)明显比2013探测所得体积含水率(平均为0.082 cm3/cm3)高，这与长测线体积含水率探测结果一致，但与长测线不同的是，正方形区域内两次探测所得的含水率的空间分布十分相似，含水率只是在数量上有所变化。从图中可以看出植被基本都分布在含水率较高的区域(黑色实线范围内，含水率>8 cm3/cm3)，其中含水率最高区域位于沙柳覆盖区，如25/24，36/33处，含水率达到0.12 cm3/cm3；含水率较低区域一般都为裸砂地(黑色虚线范围内,含水率<0.06 cm3/cm3)，如4/0、16/0、8/12处，含水率小于0.04 cm3/cm3。另外，两年探测正方形区域内地表沉降平均为1.84 m，沉降较为均匀。

4 结论

利用探地雷达对中国西部煤矿开采区开采前后地表土壤含水率进行了探测，利用地质统计方法分析了土壤含水率的变化情况，并得出以下结论：

(1) 干旱半干旱地区使用探地雷达能够得到与取样烘干精度相近的土壤含水率，且探地雷达具有快速、连续的探测优点，这是传统方法所不能做到的，在以后大范围土壤含水率探测中探地雷达应该得到更加广泛的应用。

(2) 在长测线范围内地表土壤含水率及其分布主要受控于海拔高度和沙丘地形，未开采区地表土壤体积含水率2013、2014年两次探测变化不大，分别为0.084 cm3/cm3、0.082 cm3/cm3，而开采区中心位置地表土壤含水率在地下煤层开采后显著增大，2013、2014年两次探测分别为0.099 cm3/cm3、0.083 cm3/cm3。

(3) 40×40 m正方形区域研究结果显示，2014年探测地表土壤含水率明显大于2013年，两次探测分别为0.096 cm3/cm3、0.082 cm3/cm3。在含水率较高的区域基本都分布在植被覆盖区，其中含水率最高区域位于沙柳覆盖区，裸沙地地表土壤含水率较低。

[1] J.A.Huisman, S.S.hubbard, J.D.Redman et al. Measuring soil water content with Ground Penetrating radar: A review [J]. Vadose Zone, 2003,15: 476-477.

[2] L.Weihermuller, J.A.Huisman, S.Lambot, et al. Mapping the spatial variation of soil water content at the field scale with different ground penetrating radar techniques. Journal of Hydrology 2007, 340:205-216.

[3] Zhengfu Bian, Shaogang Lei, Hilary,Inyang, et al. Integrated method of RS and GPR for monitoring the changes in the soil moisture and groundwater environment due to underground coal mining[J]. Environ Geol, 2009, 57:13-14.

[4] 宋亚新.神府-东胜采煤塌陷区包气带水分运移及生态环境效应研究[D].北京：中国地质科学院，2007.

[5] 赵红梅.采矿塌陷条件下包气带土壤水分布与动态变化特征研究[D]. 北京：中国地质科学院，2006：131-132.

[6] 杜善周.神东矿区大规模开采的地表移动及环境修复技术研究[D].北京：中国矿业大学(北京)，2010.

[7] 姚国征.采煤塌陷对生态环境的影响及恢复研究[D].北京：北京林业大学，2012.

[8] Ing Kratzseh H.Mining Subsidence Engineering [J].Environ Geo1 Watersci, 1986, 8(3):133-136.

[9] Booth CJ.Groundwater as an environmental constraint of longwall coal mining [J]. Environ Geol，2006，49:796-803.

[10] 臧荫桐.采煤塌陷后风沙土理化性质变异性研究[D].呼和浩特：内蒙古农业大学，2009.

[11] J AHuismana, C Sperlb, W Boutena, et al. Soil water content measurements at different scales accuracy of time domain reflectometry and ground- penetrating radar [J]. Journal of Hydrology, 2001,245:48- 58.

[12] Hubbard Susan. Mapping the volumetric soil water content of a California vineyard using high-frequency GPR ground wave data[ J] . Geophysics, 2002, 21:552-559.

[13] Johnson.L., Nemani,R, Pierce,L, et al. Toward the improved use of remote sensing and process modeling in California’s premium wine industry[J]. Proceedings of IEEE International Geoscience and Remote Sensing Symposium. 2000, 14: 363-364.

[14] AnitaTuresson. Water content and porosity estimated from ground-penetrating radar and resistivity[J]. Journal of Applied Geophysics, 2006, 58:108-110.

[15] K.Grote, S.Hubbard, Y.Rubin. Field-scale estimation of volumetric water content using ground-penetrating radar ground wave techniques [J]. Water Resources Research, 2003,15:145-147.

[16] I.A.Lunt,S.S.Hubbard,Y.Rubin. Soil moisture content estimation using ground-penetrating radar reflection data [J]. Journal of Journal of Hydrology, 2005, 18: 307:267.

[17] Y.G.N.S.Wijewardana, L.W.Galagedara. Estimation of spatio-temporal variability of soil water content in agricultural fields with ground penetrating radar[J]. Journal of Hydrology, 2010, 391: 24-25.

[18] Yanfang Qin, Xi Chen, Kefa Zhou, et al. Ground-penetrating radar for monitoring the distribution of near-surface soil water content in the Gurbantunggut Desert [J]. Environ Earth Sci, 2013,70:2891-2892.

[19] Du S.,Rummel P. Reconnaissance studies of moisture in the subsurface with GPR[C]//Fifth International Conference on Ground Penetrating Radar. 1994: 1241-1248.

[20] Du S. Determination of water content in the subsurfacewiththe ground wave of ground penetrating radar[D]. PhD Thesis Ludwig-Maximilians-Universitat, Munchen, 1996,117.

[21] Sperl C. Erfassung der raum-zeitlichen Variation des Boden-wassergehaltes in einem Agrarokosystem mit dem Ground-Penetrating Radar[D]. PhD Thesis Technische Universitat München, München, 1999: 182.

[22] TOPP G C, DAVIS J L, ANNANA P. Soil water content: measurements in coaxial transmission lines [J]. Water Resource, 1980, 16:574- 582.

[23] Lv Yizhong ,Hu Kelin , Li Baoguo. The spatio-temporal variability of soil water in sand dunes in mowusu desert. [J]. Acta Pedologica Sinica, 2006, 43(1):153-154.