土壤含水率是影响土地复垦的主要因子之一，准确监测土壤含水率时空分布对生态修复有重要意义。传统的土壤含水率测定方法通过开挖土壤剖面或钻孔取样，实验室测量或通过便携式仪器如时域反射仪(time-domain reflectometer，TDR)等获取，虽然这些传统探测方法能够得到精确的地下土壤的结构特征和土壤含水率，但仅限于小尺度范围内^[1-2]。要得到较大区域范围的土壤含水率时空分布不仅费时费力，而且还会出现以点代面的情况。近年来，探地雷达(ground penetrating radar，GPR)已被应用于中尺度范围内土壤含水率的探测研究中，其中探地雷达地面波法、折射波法等探测土壤含水率均需要计算波速，进而获取介质的介电常数，再由Topp公式等经验公式计算含水率值^[3]。常见的经验并不适用于所有类型的土壤，而且会造成误差的二次传递，从而影响含水率测量的精度。相比于雷达波波速分析求土壤含水率方法，Pettinelli等^[4]提出了一种利用平均振幅包络倒数(average envelope amplitude，AEA)方法，计算认为雷达早期信号的平均振幅包络与浅部地层介电常数呈负相关关系，从而计算出土壤的介电常数，再由Topp经验公式计算含水率值^[5]。Di Matteo等^[6]利用AEA法分析室内模型和野外数据的早期信号振幅属性的变化，证明该方法能够很好地反映土壤含水率的变化；Pettinelli等^[4]在野外对AEA方法进行了验证，结果表明雷达波早期振幅包络与介电常数有着高度的相关性，但并没有对该方法的适用条件进行阐述。

在国内，吴志远等^[7]利用AEA方法在野外进行地表含水率探测，找出与介电常数相关系数大的振幅包络倒数平均值，利用Topp公式计算出土壤含水率并与TDR和钻孔取样探测的地表含水率进行对比，表明使用探地雷达AEA方法能够能得到与TDR及取样烘干实测含水率相近的土壤含水率。乔新涛等^[8]采用AEA法对复垦区农田地块土壤含水率进行反演，并与雷达测线上TDR法测量所得土壤含水率结果进行对比，认为研究的4种关系模型中对数模型的拟合结果最好。

前人研究是利用雷达信号包络平均值求出介电常数，然后再采用Topp公式等经验公式进一步估算含水率，这容易将误差进行二次传递。本文是在前人的研究基础上，通过物理模拟试验研究第一峰值信号、第一正半周期和第一半周期3种典型的GPR早期信号^[6]平均振幅包络值和土壤含水率的直接相关性，拟合找出3种早期信号中相关性最强的早期信号，快速获取大范围内土壤含水率值。

1 方法原理 1.1 数值模拟试验

利用Giannopoulos^[9]开发的开源软件GprMax建立土壤模型控制体积含水率(volume moisture content，VMC)分别为5、10、15、20和25 cm³/cm³，激励源选择主频为400 MHz的雷克子波，模拟出不同含水率的土壤对应的雷达剖面，选取单道信号振幅，分析土壤含水率的变化对雷达波早期信号振幅的影响。

1.2 土壤类型及物理模拟试验 1.2.1 土壤类型测试

利用紫外可见分光光度计测量试验土样颗粒粒径，黏粒(< 2 µm)、壤粒(2 ~ 20 µm)、砂粒(20 ~ 2 000 µm)分别为4.5%、48.5%、47%。对照国际制土壤质地分类标准^[10]确定试验土样为壤土，再根据中国土壤分类系^[11]，结合农业生产方面分类，试验土样应属六级分类制中土类红壤。

1.2.2 土壤含水率物理模拟

共设置7种不同含水量的物理模型，模型箱大小为1.2 m × 0.8 m × 1 m。分别用300、400、900 MHz 3种不同频率的中高频天线沿测线方向(图 1)对模型箱土壤含水率采用共偏移距方法进行探测，目前，国际上应用GPR探测土壤含水率时采用225 ~ 900 MHz的中高频天线^[12]，225 MHz以下的低频天线尺寸较大，分辨率较低不适合土壤含水率物理模拟试验，中高频天线在保证探测精度的同时也兼顾一定的探测深度。测完后在模型剖面方向打标处从上至下间隔10 cm依次用环刀连续取10个样以及对应的平行样。取样结束后从模型箱里取出土样用洒水壶喷水，并搅拌均匀，再装入模型箱按照上述试验数据采集取样步骤完成其他组试验。第一组试验对采集的原状土样进行探测。由于水量人为控制，后面几组试验视土壤湿度控制喷水量，采集完7组不同含水量的雷达剖面，最后将所取样品带回实验室用烘干法测样品的土壤含水率，以及土壤容重等参数。

1.2.3 土壤介电常数物理模拟

利用GPR反射波法，测量电磁波在不同土壤含水率中的传播速度，再根据速度与介电常数之间的关系确定不同含水率土壤的介电常数，将求得的介电常数代入Topp公式并与烘干法测量的结果进行对比。本试验选择一个0.45 m × 0.32 m × 0.25 m的小塑料箱建立物理模型(图 2)，并在塑料箱底部放置一块金属板，以便实现对电磁波的全反射，获取电磁波在已知厚度的土壤层中的双程旅行时，用于计算土壤的介电常数。

1.3 反演土壤含水率原理 1.3.1 反射波法反演土壤含水率

利用GPR反射波法测量不同含水率的土壤介电常数，电磁波的传播速度会随着介质介电常数的增加而减小，当遇到良导体时，电磁波便会发生全反射。雷达波在土壤中的传播速度与土壤的介电常数存在以下关系，即：

$ v = \frac{c}{{\sqrt s }} $

(1)

式中，c为雷达波在真空中的传播速度(0.3 m/ns)，ε为土壤的相对介电常数。

因此，只需要求得电磁波在土壤中的传播速度，便可通过式(1)计算得到土壤的相对介电常数。利用GPR反射波法，在已知测量土壤厚度h的情况下，根据电磁波在介质中的双程旅行时间t，可以根据式(2)求得电磁波在土壤中的传播速度v。

$ v = \frac{{2h}}{{\Delta t}} $

(2)

由公式(1)、(2)便可计算出土壤的介电常数：

$ \varepsilon = {\left( {\frac{c}{v}} \right)^2} = {\left( {\frac{{c \times {\text{d}}t}}{{2h}}} \right)^2} $

(3)

将公式(3)计算得到的介电常数带入Topp公式(公式(4))等θ–ε经验公式便可反演得到土壤体积含水率。

$ \theta {\text{ = }} - {\text{5}}{\text{.3}} \times {\text{1}}{{\text{0}}^{ - {\text{2}}}}{\text{ + 2}}{\text{.92}} \times {\text{1}}{{\text{0}}^{ - {\text{2}}}}\varepsilon - {\text{5}}{\text{.5}} \times {\text{1}}{{\text{0}}^{ - {\text{4}}}}{\varepsilon ^{\text{2}}}{\text{ + 4}}{\text{.3}} \times {\text{1}}{{\text{0}}^{ - {\text{6}}}}{\varepsilon ^{\text{3}}} $

(4)

式中：θ为土壤体积含水率(cm³/cm³)，ε为土壤介电常数。

1.3.2 AEA法反演土壤含水率

电磁波在介质中传播时，振幅会受到周围介质中电磁波性质的影响，随着传播距离z增加，其振幅A相对于初始振幅A₀呈指数衰减^[13]，具体计算公式如下：

$ A=A_{0}\text{e}^{–αz} $

(5)

式中：α为衰减常数，$ \alpha = \frac{\sigma }{2}\sqrt {\frac{\mu }{\varepsilon }} $；σ为介质电导率(S/m)，μ为介质磁导率(H/m)，ε为介电常数。

对于低损耗介质α是与频率无关的，因此由上式可以看出介电常数ε对振幅影响很大。水的介电常数为81，土壤的介电常数一般为3 ~ 25，利用土壤和水的介电常数差异得到不同土壤含水率对应的电磁波早期振幅信号，通过希尔伯特变换^[14]得到振幅包络，设GPR得到的原始连续信号为x(t)，对其进行希尔伯特变换，其变换公式如下：

$ H(t) = x(t) \times h(t) = \frac{1}{\pi }\int_{ - \infty }^{ + \infty } {\frac{{X(\tau )}}{{t - \tau }}} d\tau $

(6)

式中：$ h(t) = \frac{1}{{\pi t}} $。

信号数据经过希尔伯特变换后，相位谱要作90°相移，相当于一次滤波，以x(t) 为实部，以H(t) 为虚部，构建解析信号^[15]：

$ f(t) = x(t) + iH(t) $

(7)

则振幅包络为：

$ E(t) = |f(t)| = \sqrt {{x^2}(t) + {H^2}(t)} $

(8)

求得的振幅包络为正值，它使GPR单道信号变得简化，让振幅更容易被解释。将变换后的数据用matlab程序计算得到雷达各单道振幅包络值，选取取样附近的单道数据求得平均振幅包络倒数。

2 结果与分析 2.1 数值模拟结果

模拟出的雷达波早期信号振幅的变化情况以及通过希尔伯特变换后得到的振幅包络如图 3所示，从中可以看出，当土壤含水率较低时，微小的含水量变化就会导致振幅有很大的波动，而当土壤含水率较高时，即使含水量变化很大但振幅的波动并不明显。因为土壤含水率越高，电磁波衰减越快，这也就导致了GPR实际探测精度降低^[16-20]。

根据GPR早期信号对应的时间窗口段内的所有单道振幅包络值求平均，然后选取第一峰值信号、第一正半周期和第一半周期3种早期信号进行研究。如图 4所示，利用gprMax模拟这3个GPR早期信号对应的时窗内单道波形及相应的希尔伯特振幅变换。

2.2 物理模拟结果 2.2.1 土壤体积含水率

利用烘干法测得土样的质量含水率，根据土壤质量含水率、土壤容重和土壤体积含水率之间的关系，利用式(9)计算得到土壤体积含水率，并作出土壤体积含水率分布图(图 5)。

$ {\theta _v} = {\theta _w} \times {\rho _b} $

(9)

式中，θ_v为土壤体积含水率(cm³/cm³)，θ_w为土壤质量含水率(g/g)，ρ_b为土壤容重(g/cm³)。

2.2.2 土壤介电常数

图 6为土样雷达剖面及第500道单道图，t₀为直达波双程旅行时，t₁为电磁波到达金属板位置双程旅行时，$ \Delta t $为电磁波在厚度为h的红壤土中传播的双程旅行时。

通过探测不同含水率土壤得到不同的雷达剖面，提取取样附近雷达剖面单道数据读取电磁波在不同含水率土壤中传播的双程旅行时，根据公式(3)得到不同含水率土壤的介电常数，将计算得到的介电常数带入Topp公式得到土壤预测体积含水率。介电常数试验相关参数见表 1。

从图 7A中可以看出，Topp公式整体高估了土壤的体积含水率，对实测土壤含水率和Topp公式计算得到的土壤含水率进行误差分析（图 7B），平均相对标准差为1.069%。说明Topp公式多适用于描述颗粒粒径较大的土壤类型(如砂土)的介电常数与体积含水率的关系，此外Topp公式需要计算介电常数来反演土壤含水率，这会造成误差的二次传递，这也是Topp公式高估了土壤体积含水率的一个重要因素。

2.2.3 AEA法试验结果

首先对模型箱内采集的雷达数据通过零点校正、背景去噪、一维滤波等常规处理后，提取单道波形，对每种天线采集的7组数据对比分析。如图 8所示，通过归一化处理，可以看出不同含水率土壤早期振幅信号起跳点不同，含水率越低，起跳点越早，振幅越大，如图中黑色虚线代表模拟试验第一组对原状土探测所得的雷达波单道振幅^[21]。随着时间的推移，电磁波信号衰减到越来越弱，对应的雷达波单道振幅信号越来越小，振幅不再有较大波动，雷达波单道振幅曲线接近重合，意味着无法监测深部土壤含水率。

对实际探测单道波形进行希尔伯特变换，从图 8可以看出通过希尔伯特变换后的振幅变得更加简化。读取第一峰值信号、第一正半周期和第一半周期对应的时窗范围，并进行振幅包络计算，得到平均振幅包络倒数。选取第二组数据作为验证，拟合其他6组雷达波振幅包络倒数与烘干法测得的土壤体积含水量之间的关系^[22]，如图 9所示。

根据图 9可以直观地看到第一正半周期内的振幅包络的倒数与土壤体积含水率拟合效果最好，表 2给出了第一正半周期内振幅包络的倒数与土壤体积含水率的拟合关系式。

将第二组试验得到的平均振幅包络倒数代入表 2中的拟合关系式得到土壤体积含水率，并与烘干法测得的土壤含水率进行比较，作出3种早期信号内土壤体积含水率误差分析图，从图 10可以看出，第一正半周期内误差最小，体积含水率平均相对标准差为0.483%，因此可以利用第一正半周期振幅包络倒数计算红壤土的体积含水率。

3 结论

1) 对于3种不同频率的天线，研究的3种时段的早期信号中第一正半周期内雷达波平均振幅包络的倒数与土壤体积含水率有较好的线性关系，300、400、900 MHz天线探测数据计算的土壤含水率与烘干法测得的土壤含水率标准差分别为0.506%、0.289%、0.656%，平均相对标准差为0.483%。

2) 土壤介电常数试验结果表明Topp公式高估了土壤的体积含水率，实测土壤体积含水率与Topp公式计算得到的体积含水率平均相对标准差为1.069%。结果表明利用AEA法研究第一正半周期内雷达波平均振幅包络倒数得到的土壤含水率较为精确。

3) 本文研究不通过估算土壤介电常数就可以直接利用GPR第一正半周期振幅包络倒数快速计算土壤体积含水率，对于农业生产、土地复垦具有重要的指导意义。