0 引言

水是整个生态系统赖以生存的重要资源.土壤含水量关系植物的生长，也关系到地基强度进而影响建筑物的稳定.因而在环境、生态、安全等许多领域，土壤含水量的研究得到越来越多的关注.及时探明土壤含水状况，以便依此做出科学的决策或采取合理的措施具有十分重大的意义(冉弥等，2010).传统探测含水量方法包括烘干法、核技术法、时域反射法、遥感测定法等.其中前三种方法适用于小尺度测量，只能测定单个点的体积含水量，对土壤层有破坏，第四种方法适用于大范围的测量，但分辨率较低，且要求植被覆盖率低.与传统方法相对比，探地雷达方法作为一种获取中尺度的土壤含水量的方法，能够在保证原位测定精度同时，便捷地获取土壤含水量信息，提高测量效率(雷少刚和卞正富，2008).

近年来，探地雷达方法逐渐被应用于浅地表探测中来(Kao et al., 2007；De Pascale et al., 2008；冯德山等，2011；丁亮等，2012；刘澜波和钱荣毅，2015；Cao et al., 2017).高频雷达波在地下介质中的传播速度主要取决于电导率σ和介电常数ε.在探地雷达的工作频率(10~2500 MHz)范围内，20 ℃的水的相对介电常数为80(Huisman et al., 2003；Galagedara et al., 2005)，空气的相对介电常数为1，而大部分土壤固相的相对介电常数为4~9(Lunt et al., 2005).因此，水是决定土壤相对介电常数的最重要的因素.利用探地雷达，可以在得到地层的埋深及厚度的同时，获取地层的介电常数及电导率等信息，进而求出土壤含水量(Topp et al., 1980；Herkelrath et al., 1991, 黄忠来和张建中，2013).随着探地雷达技术的不断发展，其测量精度也在不断提高.因此，探地雷达在评估土壤含水量方面具有广阔应用前景.

目前，应用探地雷达来测量土壤含水量的方法主要有以下四种：① 反射波法；② 地面波法；③ 地表反射系数法；④ 钻孔雷达法.其中反射波法又分为单偏移距法(Single Offset Reflection Methods)和多偏移距法(Multi-Offset Reflection Methods)(Huisman et al., 2003).单偏移距法的优点是在短时间内可以获得大面积的地下潜水面信息.但该方法只能用于获得地面到已知的反射层之间的平均含水量，无法控制深度分辨率.最常用的多偏移距法包括宽角法(Wide Angle Reflection and Refraction，简称WARR)和共中心点法(Common Midpoint，简称CMP).CMP方法工作量较大，但优点在于可以获得同一地点不同深度的电磁波传播速度的变化，进而获得不同深度的介电常数及体积含水量.本文主要研究CMP方法的应用.

本文基于速度分析原理和混合介质介电模型，建立了一套计算土壤含水量的探地雷达数据处理和分析流程.利用提出的数据处理和分析流程对多层水平介质模型的正演模拟数据进行了处理和分析，获得的结果说明了应用此套流程对探地雷达CMP数据估计土壤含水量是适用可行的.然后应用此流程对蒙古国乌兰巴托市水源区域的CMP数据进行了处理和分析，获得了该区域地下水位和土壤含水量信息，并与观测井的信息进行了对比.

1 探地雷达测量土壤含水量的理论基础 1.1 速度分析原理

在探地雷达数据处理流程中，速度是最为重要的参数，直接影响土壤含水度的评估结果.为了准确提取速度, 我们通常先对原始共中心点剖面数据进行滤波、增益等预处理，使剖面中各界面处的反射波更加清晰，叠加速度谱图像中的速度信息更易于拾取.

速度分析的前提是假设介质为水平层状介质，其反射波正常时差具备关系为

(1)

式中x_i为炮检距，t₀为零偏移距时的双程走时，v_rms为叠加速度.实际数据处理中，不可能从地震记录精确地读取反射时间.水平层状介质条件下，反射波的同相轴为双曲线，可以考虑扫描同相轴，达到速度分析的目的.

地震记录中，反射波双程走时满足关系为

(2)

因此，速度分析时，可以给定一个合适的零偏移距双程走时t₀的范围[t_0min t_0max]和叠加速度范围[v_min v_max]，确定t₀值，对每个叠加速度v_rms按照式(2) 进行扫描；然后对下一个t₀值重复以上过程；这时就可以形成一个关于t₀和v_rms的网格系统，每个网格值都被唯一确定，按网格值大小可以对其进行伪彩色填充，这就获得了速度谱.与真实地下层状介质零偏移距双程反射时间和叠加速度相同的t₀和v_rms确定的网格，其网格值将达到最大，形成能量团聚焦，可以通过拾取其对应的时间和速度，确定层速度和层厚度.实际速度分析中，应根据子波长度，合理确定一个时窗对双曲线轨迹进行扫描.常用的判别准则有平均振幅能量准则、平均振幅准则、非归一化准则、归一化准则和相似系数准则.考虑到探地雷达的信号特点和本次测量的目的，选取对大幅值干扰不灵敏的叠加类准则，即平均振幅能量准则和平均振幅准则.其定义如下.

① 平均振幅能量E定义为

(3)

式中，N为探地雷达数据道数，为叠加时窗内的采样点数，f_{i, j}为地震数字信号.当扫描速度等于均方根速度时，E达到最大值，表现为速度谱上的能量团聚焦.

② 平均振幅A定义为

(4)

理论上与平均振幅能量等效，计算量下降一些.

拾取叠加速度之后，可以使用Dix公式求取层速度，表达式为

(5)

式中v_{i, n}为第n层层速度，v_{rms, n}为至第n层底部的叠加速度，t_{0, n}是至第n层底部零偏移距双程走时，v_{rms, n－1}为至第n-1层的叠加速度，t_{0, n－1}是至第n-1层底部零偏移距双程走时.

第n层介质的层厚度计算公式为

(6)

获取了雷达波在地下传播的层速度后，我们需要将其转化为相对介电常数.在低损耗介质中，雷达波的速度与相对介电常数具备关系为

(7)

式中C为真空中电磁波速度.

1.2 混合介质介电模型

为了根据相对介电常数获得土壤含水量，许多国内外学者都进行了研究(Weiler et al., 1998；巨兆强，2005；Černý，2009).混合介质模型就是表征含水量与相对介电常数关系的一组模型(王春辉，2007).其中比较著名的有：Topp公式、折射系数(CRIM)公式、Malicki公式、Roth公式和Alharathi公式等.

Topp公式(1980)是最常用的由介质的相对介电常数获得含水量的经验公式(Stoffregen et al., 2002；Lu Q，2005).公式中相对介电常数和体积含水量具备关系为

(8)

式中θ为土壤含水量；ε_r为相对介电常数.该公式是对大量的含有不同胶结系数和含水饱和度的土壤样品进行实验后得到的经验公式.但由于实际测量中，土壤往往非常复杂，因此应用Topp公式时可能会出现较大的差别，这需要我们通过实验数据修改公式中的系数，来获得更好的拟合结果.该公式对于质地较粗的土壤有很好的计算结果(Greaves et al., 1996).

另外一个用的比较广泛的是折射系数模型(Complex Refractive Index Model，简称CRIM).这是一个半经验公式，有一定的理论支撑，经常用于电磁测井的数据解释中.它基于介电混合模型，理论上与估算多孔介质声波速度的Wiley Time-verage Equation相类似.折射系数公式的表达形式为

(9)

式中：φ(m³m^-3)表示介质的孔隙度；ε_w、ε_s和ε_a三者分别表示水、土壤颗粒和空气的相对介电常数；α与外加电场的方向有关.当电场方向与介质层平行时α=1，当电场方向与介质层垂直时α=－1，α=0.5代表是各向同性介质.经过整理可得到公式为

(10)

一般情况下，我们设定ε_a=1，α=0.5.公式(10) 即可写为

(11)

该公式对于高频段(MHz~GHz)的介电常数效果较好.

2 基于CMP数据土壤含水量估测

在实际情况中，通过正演模拟或实际测量所获得的原始共中心点剖面往往不能直接用来做速度分析.这是由于雷达波在地下传播的过程中，其能量会随深度的增加而衰减，因此原始共中心点剖面中往往只能看出第一个反射界面，而深部的反射波信息则会由于能量过低而显示不出来.因此，在速度分析之前，我们需要对原始共中心点剖面进行线性增益处理，使剖面中各界面处的反射波更加清晰.

需要提出的是，在处理野外实测数据时，由于噪声的影响，除了线性增益外，往往还需要进行时移改正和去除背景噪声处理.

图 1为CMP数据的处理流程图.通过正演模拟或实际测量获得原始CMP剖面后，我们首先对其进行线性增益处理，让有效波更加突出.然后进行速度分析，根据实际数据选取合理的时窗对双曲线轨迹进行扫描，利用平均振幅能量准则及平均振幅准则提取出叠加速度.随后根据Dix公式将叠加速度转换为层速度和层厚度，再利用雷达波的速度与相对介电常数之间的关系得到相应各层的相对介电常数，最后运用混合介质介电模型理论计算土壤含水量(Donoho and Johnstone, 1994；Nakashima et al., 2001；Turesson，2006).本文中主要应用Topp公式及孔隙度为39%下的CRIM模型计算各层的体积含水量.

2.1 多层水平介质模拟

图 2a所示为一个水平距离5.2 m，深度4.0 m的四层水平介质模型.由上到下分别为干砂岩、湿砂岩、饱和砂岩和基底，相对介电常数依次为5.0、8.0、11.0和15.0.模型中包括三个反射界面.运用雷达波仿真模拟软件GprMax2.0进行对模型进行正演，设定其的工作主频300 MHz，波源为Ricker子波，时间步长0.0236 ns，空间步长0.01 m，满足CFL条件(Courant et al., 1967).正演天线的初始偏移距为0.1 m，每次各向两侧对称地移动0.05 m.正演提取的共中心点剖面数据包含46道，中心点位于模型2.5 m处.对该模型进行上述处理得出结果如图 2所示.

应用混合介质模型求取土壤含水量：给定砂、水和空气的相对介电常数分别为ε_s=4.0、ε_w=81、ε_a=1.分别运用Topp公式及孔隙度为39%下的CRIM公式来计算土壤含水量，计算结果如表 3所示，垂向含水量剖面如图 3所示.

结果综合分析：

(1) 速度分析

该模型包含四个水平层，三个水平界面.在地下传播的过程中，雷达波的能量会随深度的增加而衰减.我们能从原始共中心点剖面图 2b中看出：反射能量振幅随深度增加而明显衰减.线性增益处理之后，三个水平界面均能很明显的反映出来，如图 2c.提取速度谱后我们得到了三个能量团，如图 2d.平均振幅能量曲线图 2e，表征不同反射时间平均振幅能量变化情况.随后根据Dix公式得到层速度曲线图 2f，其速度分析结果如表 1所示.经过分析误差我们可以发现，表 2中介电常数的相对误差最大为10.1%，层速度的相对误差最大5.4%，分析结果与正演参数的拟合效果较好.通过该模拟实验可以看出，所设计的速度分析算法可以很好的提取出叠加速度，并且该速度谱算法能够很好地降低多次波的影响，反演结果误差较小，效果可靠.这表明，应用探地雷达反射波法可以提取地下介质层速度信息，可以对土壤含水量进行评估.

(2) 含水量分析

依据速度分析得到的结果，我们分别运用了Topp公式、CRIM模型计算正演层状介质模型含水量.从表 3可以发现，孔隙度为39%的条件下，CRIM模型计算的含水量与Topp公式计算的含水量基本一致.图 3为垂向含水量剖面，含水量随着深度的增加而变大，这与设计的模型含水量变化一致.

2.2 野外实测数据处理与分析

待测水源区位于蒙古国首都乌兰巴托市区的南部，该区为全市居民提供生活用水，地下水丰富.测量位置在10号水井附近，测量时水井已停止泵浦十几小时并且地下水位已趋于稳定.此时，用时域反射仪(Time Domain Reflectometry，简称TDR)测量地表土壤速度获得的结果为0.159 m/ns，从井中观测到的地下水位为6.58 m(Liu and Sato, 2012).

共中心点数据由瑞典ramac探地雷达系统采集(见图 4).其中天线选用的是中心频率为100 MHz的非屏蔽天线，时窗为267.13 ns，时间域采样间隔为0.2609 ns，采样点数为1024，道数为91道，天线的初始偏移距为0 m，每次各向两侧对称地移动0.1 m.所测得的原始数据见图 5a.

结果综合分析：

(1) 速度分析

此次测量中，探地雷达主频很高，地下介质可以被认定为是低损耗介质，雷达波速传播过程中不会发生改变.首先对原始数据进行时移改正和去除背景噪声处理，得到剖面如图 5b所示，线性增益后得到图 5c所示剖面，从图 5c中可以看出第一个反射同相轴的零偏移距双程走时大约为40 ns，所以速度谱中大于40ns的能量团才有意义.经过速度分析后绘制出速度谱图像如图 5d所示，从图中可以看出速度谱中有三个有意义的能量团，对应三个反射界面，这对应于平均能量曲线图 5e中的三个能量峰值.拾取速度谱能量团极值对应的零偏移距双程走时和对应的叠加速度，利用Dix公式计算层速度和层厚度，如图 5f.其中, 第一层速度为0.163 m/ns，这与时域反射仪(TDR)测量获得的结果0.159 m/ns相比，误差仅为2.5%，证明由探地雷达方法获取的速度是可靠的.根据层速度计算介质的相对介电常数，结果如表 4所示.由介电常数可以计算土层的含水量，如表 5所示.

(2) 含水量分析

根据计算出的相对介电常数, 分别利用Topp公式和孔隙度为39%下的折射系数公式求出土层含水量，如表 5所示，从表中可以看出深度越深，含水量越高.从垂向含水量剖面图 6可以看出，第三层含水量突然升高，这可以解释为第二个反射界面是潜水面(Lu et al., 2007；Liu and Sato, 2012)，对应于速度谱中89.8 ns的能量团.可以得出潜水面的深度为6.66 m，这个结果与水井停止泵浦十几小时后水位的观测结果6.58 m(Liu and Sato, 2012)相比，相对误差为1.2%，这也表明获得的土壤含水量信息是可信的.

3 结论 3.1

本文讨论了探地雷达测量土壤含水量的基本原理，探地雷达数据的预处理流程.设计了平均振幅能量和平均振幅速度分析程序，并将这一速度分析算法应用于探地雷达正演数据和实测数据.正演结果表明速度分析算法的误差较小，可以用于对实测数据的分析.

3.2

应用本文设计的数据处理流程，对蒙古国首都乌兰巴托南部10号水井处于非泵浦时的探地雷达实测CMP数据进行分析，将土层划分为三层.求取的土壤含水量，在第三层突然升高，可以认为是潜水面，得出其深度为6.66 m，这个结果与水井停止泵浦十几小时后水位的观测结果6.58 m，相对误差为1.2%，这也证明了探地雷达方法得出的土壤含水量信息是可信的，说明探地雷达方法对确定土壤含水量具有很大的应用潜力和研究价值.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！