不同充填基质复耕土壤含水率反演与特征分析

王 阳，周 涛，程 琦，宋桂芳,张世文*

(1.安徽理工大学 地球与环境学院，安徽 淮南 232001；2.安徽理工大学 空间信息与测绘工程学院，安徽 淮南 232001；3.安徽理工大学 安徽省高潜水位矿区水土资源综合利用与生态保护工程实验室，安徽 淮南 232001)

随着我国城镇化步伐的加快，耕地后备资源不断减少，人地矛盾愈演愈烈，加之人类不合理的耕地利用方式，耕地质量日趋下降，耕地供需平衡压力日益加大，土地复耕已成为落实耕地占补平衡的一项有力举措。其中，土壤含水率是判断复耕、复种效果指标体系中的重要一环[1]。快速、精准、无损探测复耕土壤含水率对复耕农业水资源管理、复耕效果监测以及复耕农田水文研究具有重要意义。目前，含水率测量主要使用烘干法，虽然取得的结果准确性较高，但烘干法效率低下，大范围探测较为困难。而中子射线法和TDR等方法操作复杂且需要进行事先标定[2]。卫星遥感等空间信息处理技术测量含水率的范围较大(通常是对大于1公顷的地块)，对于构造较为复杂、充填基质多样的新增复耕耕地，无法满足复耕质量检测精度的要求。探地雷达作为一种新型的近地微波物探技术，逐渐运用于资源环境监测，对于含水率的探测较为快速和准确，既弥补了烘干法在效率方面的不足，也能够对充填多样、层次复杂的土壤含水率产生响应。

近年来，探地雷达在土壤含水率探测方面的研究进展较快，同时分析方法已逐渐从波速、距离估计法预测向电磁波属性特征分析转变。Wang等[3]利用探地雷达对青藏高原季节性冻土区进行含水率探测，结合振幅分析、波形分析以及电磁波衰减属性等分析方法，研究得出振幅属性以及瞬时频率等属性与含水率间具有一定的相关性。Shen等[4]通过分析电磁波不同特征属性与有机碳、含水率的响应关系，研究得出所有特征属性都与土壤含水率呈现显著相关，其中振幅分析中能量属性与含水率相关性最高。Fabio等[5]通过室内实验分析探地雷达谱峰频率分量与土壤中黏土含量的响应关系，实验结果表明电磁波谱峰频率分量与土壤中黏土含量具有较好的响应关系，随着黏土含量的增加，频率峰值逐渐减小。Ferrara等[6]在自然条件(野外环境)下通过研究电磁波信号的早期振幅包络平均值与不同深度土壤含水率的响应关系，结果表明早期振幅包络平均值能够准确反演出0～50 cm范围内的土壤含水率。罗古拜[7]利用GPR反射波法探测结果以及烘干法结果进行建模，建立了介电常数和含水率的模型。乔新涛等[8]采用包络振幅平均值对黄土高原复垦矿区土壤进行探测，对黄土区土壤水分反演精度较高。程琦等[9]使用AEA法对以粉煤灰为充填基质的复垦土壤进行探测，发现电磁波特征波段与耕层含水率存在较强响应，反演拟合效果较好。目前研究大多基于自然土壤下的含水率探测，而对于不同充填基质下复耕重构土壤的探测鲜少。复耕重构土壤因受到人为干预导致其土壤环境较正常土壤更为复杂，内部充填基质不再均一，雷达信号受影响因素颇多。而新构土壤含水量的研究对于掌握矿区土壤水土保持生态修复十分关键。郭婷婷[10]以矿区重构土壤含水率为研究对象，结果表明不同粒径的粉煤灰、煤矸石呈现不同的土壤水分下渗规律，从而影响土壤田间持水量。胡振琪等[11]在自然条件下开展不同夹心土层位置的黄河泥沙夹层式充填复垦，结果表明土壤水分累计入渗量随夹层位置的变化而变化，HYDRUSHY-1D结合含水量变化反演土壤水力参数精度较高。不同的土壤重构、充填方式对于土壤含水率保持至关重要。

因此，为探究不同充填基质土壤环境下探地雷达含水率探测问题，本文提出了一种利用零点调节、一维滤波、道间均衡、背景消除、小波去噪、FIR滤波、Hilbert变换等信号处理流程。通过构建复耕耕地不同充填基质下土壤含水率与振幅属性的关系模型，并开展野外试验探究模型适用深度以及不同充填基质下土壤含水率差异规律。然后选取最优地统计插值方法对复耕区水分空间变异性进行了分析。研究结果可在复耕农业水资源管理、复耕效果监测以及复耕农田水文研究等方面提供理论与技术支撑。

1 实验设计与方法原理

1.1 实验设计

实验开展在安徽省淮南市安徽理工大学试验区，属亚热带季风性气候，土壤质地为黏土，土壤容重偏大，介于1.3～1.6 g/cm3。实验时间为6～7月，实验选择在最后一次降雨的3 d后开展。设置2个1 m×1 m的实验样方，为还原真实复耕耕地地下充填物的复杂与多样性，填埋物采用粉煤灰、建筑用红砖、瓷砖等。A样方以瓷砖和红砖分别充填；B样方以粉煤灰和红砖碎料分别充填，覆土厚度均为40 cm。具体充填方式、方位、雷达测线及取样点如图1所示。雷达天线中心频率为500 MHz-Shield，时窗设置为50ns。由于充填层厚度要大于最小垂直分辨率(500 MHz中心频率换算后垂直分辨率约为5 cm)，故充填基质厚度统一设置为20 cm厚。充填完毕后，等待土层的自然沉降，于一个月后对试验区进行雷达探测。雷达探测后重新开挖坡面，在每种充填模式下延测线方向探测4个层次的土壤容重和含水率，分别为0～10 cm、10～20 cm、20～30 cm和30～40 cm。为保证数据的准确性，当天进行实验室烘干和测量记录。

1.2 原理与方法

(1)雷达探测原理。探地雷达(Ground Penetrating，GPR)是对地下物体发射高频电磁波(1 MHz～10 GHz)的一种物探方法。由于电磁波在传播过程中遇到不同介质时会发生波的折射和反射，当以耕层土壤为研究对象时，会因为不同的充填基质土壤介质产生不同的反射信号。

电磁波速度与介电常数的关系式为

(1)

本实验电磁波垂直入射，表达式为

(2)

(2)数据预处理。在得到每个剖面的探测图像后，先用Reflexw软件对每个图像都进行零点调节、一维滤波、道间均衡、背景消除、小波去噪等方法处理。零点调节是对雷达剖面图像进行上下移动，移动到大气与地面的分界面，切除地上部分，保留地下部分。一维滤波的作用主要是压制干扰信号提高信噪比，提取地下介质的响应特征信号。道间均衡通过道间的相关性进行信号的加强，对于道间不相关的信号，处理后信号减弱，此方法对雷达剖面图像中相关性较好的弱信号有极大的增强作用。背景消除是为了凸显异常信号，消除背景噪声和由于阻抗不匹配产生的驻波干扰，达到压缩水平干扰信号，提取有用信号的目的。小波去噪也是用于压制噪声信号，提高有用信号。水分反演部分预处理主要是利用Matlab软件，发现第一个正半周期是信号与土壤介电特性相关性最高的部分，能够使电磁波信号的信噪比最大化，并使浅层界面反射的干扰最小化。因此只需要进行去直流漂移处理与时间零点校正，其中去直流漂移目的是使有效信号不受漂移现象的影响。

(3)FIR(Finite Impulse Response)滤波。FIR滤波系统函数为H(z)：

(3)

(4)

h(n)=ω(n)×hd(n)，

(5)

根据时域相乘关系与频域卷积关系得

(6)

式中,H(ejω)为FIR滤波器频率响应；WR[ej(ω-θ)]为窗谱。

(4)提取GPR信号振幅包络值。探地雷达通过500 MHz-Shield天线对耕层发射电磁波。土壤中含水率的变化是影响雷达反射波振幅的主要原因。利用Hilbert变换后解析信号得到信号振幅包络值。公式如下：

(7)

解析信号为

R(t)=x(t)+iH(t),

(8)

则其振幅包络为

(9)

(5)土壤干容重与体积含水率测定。为验证AEA法反演含水率的效果，在每个样方同一深度，利用环刀均匀提取3份样品，装入密封袋中。提取当天在实验室将样品放置在烘干机内，在105 ℃环境下烘干24 h，烘干后得到土壤干重，取平均值。然后利用以下公式计算出干容重和体积含水率：

(10)

(11)

θv=θm×ρ,

(12)

式中,θv、θm分别为土壤样品的体积含水率和质量含水率；mwet和mdry分别是土壤湿重和土壤干重；v是环刀体积，为100 cm3；ρ为土壤样品的容重，单位为 g/cm3。

2 结果与分析

2.1 基于AEA法的含水率反演模型优选

为了探究雷达电磁波与含水率的响应关系，提取不同测线单道波形，截取特定时窗下雷达信号波，经预处理后，求得振幅包络值，并将测线处附近各测点的振幅包络值取平均值。随机选取测线上各个层次20个样点为建模集，10个样点为验证集，建立多种回归模型如表1所示。从建模集的R2来看，不同拟合模型的精度差异较小，对数拟合模型精度最高，R2为0.939；线性拟合模型精度最低，为0.889。多项式拟合模型RMSE最低，为1.987。从验证集的R2来看，拟合效果均在0.89以上，说明基于0～12 ns的波段提取的信号振幅包络平均值与土壤含水率有较强的响应关系，通过AEA法对含水率进行反演是较为精准的。从验证集RMSE来看，对数拟合模型RMSE最低，为1.707。土壤含水率反演模型优选如图3所示。图3a是振幅包络平均值与实测含水率的4种回归拟合，图3b是含水率实测值与模型预测值对比图，其中横轴表示10个验证集点，纵轴表示含水率。其中，对数拟合模型的预测效果最佳。

表1 AEA-1含水率建模精度

2.2 基于AEA法的含水率反演

为探究AEA法对不同层次土壤含水率的识别能力，使用以上确定的对数模型分别对试验区土壤剖面0～10 cm、10～20 cm、20～30 cm和30～40 cm各层的含水率进行反演，结果如图4所示。由图4可知，从R2来看，0～10 cm土壤含水率反演精度最高，为0.89；30～40 cm反演精度最低，为0.65。以RMSE作为反演精度标准，0～10 cm土壤含水率反演精度最高，为1.82；20～30 cm土壤含水率反演精度最低，为3.50。在20～30 cm和30～40 cm反演中，均存在部分差异过大点位。基于0～4 ns和4.2～8.1 ns的波段提取的信号振幅包络平均值与0～10 cm和10～20 cm土壤含水率有较强的响应关系，通过AEA法对含水率进行反演是可行且较为精准的。8.1～11.6 ns和11.6～15.4 ns波段提取的信号振幅包络平均值与20～30 cm和30～40 cm土壤含水率响应关系较差。

为更直观体现不同层次含水率的差异及反演模型的验证精度，以A样方为例，以对数回归模型作为矩形区域雷达数据反演模型，分别进行各个层次土壤实测值与AEA预测值的插值，其中图示颜色越深表明含水率越高。应用Sufer 14软件进行克里金网格插值得到该实验样方各个层次的含水率平面分布图，如图5所示。图5中分别为烘干法实测含水率分布图和雷达信号反演含水率分布图。根据插值结果，两者含水率分布较为接近，反演效果较好，但是存在部分差异较大现象，其中20～30 cm和30～40 cm两者之间差异性相较于上两层较大。

2.3 不同充填基质下水平、垂直向含水率规律分析

(1)不同充填基质垂直向土壤含水率规律。为探究基于AEA法的不同充填基质下垂直向水分演变规律，于两个样方各个充填基质上布设测线，以确定的反演模型进行反演并揭示规律。运用Origin软件绘制不同基质测线中3道雷达电磁波所反演的含水率如图6所示。从图6a中可以看出，在沉降稳定后，不同充填基质下土壤各层次含水率均遵循随着深度的增加，含水率上升这一趋势。从图6b正态曲线来看，0～10 cm的随机变量正态分布较为密集，σ较小；30～40 cm的随机变量正态分布较为疏松，σ较大。各种充填模式下的土壤含水率遵循深层的离散程度比浅层离散程度大。从A样方与B样方相同的充填基质红砖可以看出，由于沉降的作用，B样方的含水率要高于A样方的含水率。通过对比不同充填基质下不同深度土壤含水率，以粉煤灰为充填基质的土壤含水率在各层均大于其他充填基质，其他不同充填基质下的含水率差异较小。

(2)不同充填基质水平向土壤含水率规律。为探究基于AEA法的不同充填基质下水平向水分演变规律，于两个样方各个充填基质上布设成散射形由内向外测线。提取测线上雷达信号电磁波段，以确定的反演模型进行反演并揭示规律。并将反演的结果以插值的形式表现，不同充填基质下表层土壤含水率分布图如图7所示，颜色越深表示含水率越高。从图7中可以看出，不同充填模式下，样方周围的土壤含水率大体上遵循距离样方越远，含水率越低的规律。其中，从A样方的瓷砖与红砖充填可以看出，以红砖为充填基质的土壤含水率要高于瓷砖；从B样方的红砖充填与粉煤灰充填可以看出，以粉煤灰为充填基质的土壤含水率高于红砖。从测线向土壤含水率变化规律来看，以瓷砖作为充填基质的土壤含水率除充填位置中心外，与周边土壤含水率无明显差异。以红砖为充填基质的土壤含水率比周边土壤含水率高，延测线方向，呈现先减少后增加，最后与周边土壤含水率无差异。以粉煤灰为充填基质的土壤含水率比周边土壤含水率高，越往外含水率越少。通过含水率等值线分布可以看出，A样方以红砖为充填基质的土壤含水率等值线范围大约为40 cm；B样方以粉煤灰为充填基质的土壤含水率等值线范围大约为50 cm，以红砖为充填基质的土壤含水率等值线范围与A样方红砖相近，约为30 cm。

3 讨论

探地雷达通过主机及天线收发电磁波，通过截取时窗、时间、频率对土壤进行无损探测。在本次研究中发现500 MHz雷达电磁波信号特征时窗波段与土壤体积含水率具有较强的相关性，通过振幅包络值与含水率的相关性拟合可以建立反演模型。通过随机选取的建模集与验证集，分别建立多种回归模型，对数拟合模型精度最高，R2为0.939。运用以上确定的模型对不同深度土壤含水率进行反演，基于0～4 ns和4.2～8.1 ns的波段提取的信号振幅包络平均值与0～10 cm和10～20 cm土壤含水率有较强的响应关系。从正态曲线疏松程度σ来看，30～40 cm的随机变量正态分布较0～10 cm大，0～10 cm反演精度较好，通过AEA法对含水率进行反演是可行且较为精准的，R2分别为0.89和0.80；8.1～11.6 ns和11.6～15.4 ns波段提取的信号振幅包络平均值与20～30 cm和30～40 cm土壤含水率响应关系较差，RMSE均在3.20以上，与张金珠等[12]研究结果一致。

在更深层的土壤含水率反演中，AEA法的适用性依旧需要进行进一步的研究。随着土壤深度的加深，该法对含水率的反演精度逐渐减弱。在先前的研究中发现，影响AEA方法探测深度的深浅主要有两方面因素，一方面是与探地雷达的主频信号有关，而这种情况可以通过频谱补偿来提高最终的反演精度。比如程琦等[13]利用CZT法进行频谱细化处理，根据增强后的频谱信息进行含水率的反演，结果表明，反演精度有了大幅度的提升；而另一方面，时窗的选取划分也会对含水率的反演精度造成影响。由于在不同深度的反演过程中，在时窗划分中存在一定的误差，无法根据实测厚度进行准确的时窗拾取，而随着误差的传递，会引起时窗选取的偏差不断加大，使得最终的精度大幅度下降，所以在实际探测中，若需要获取较为精确的反演效果，可以在探测时取一些点位进行标定处理；或者可以采用不同方法对含水率进行联合反演，比如崔凡等[14]根据不同深度的土壤分别采用AEA法和AMRA法进行含水率的联合反演，从而获取较高的精度。

本文利用AEA法对不同充填类型的土壤含水率进行了反演并揭示规律。实验结果表明，不同充填基质下土壤各层次含水率均遵循随着深度的增加，含水率上升这一趋势，与崔凡等[15]和张世文等[16]的研究结果一致。样方周围的土壤含水率大体上遵循距离样方越远，含水率越低的规律，与Cui等[17]的研究结果一致。以粉煤灰为充填基质的土壤含水率比周边土壤含水率高，越往外含水率逐渐减少，与Tripathi等[18-19]的研究成果一致。在本次实验中，A样方以红砖为充填基质的土壤含水率等值线范围大约为40 cm；B样方以粉煤灰为充填基质的土壤含水率等值线范围大约为50 cm，以红砖为充填基质的土壤含水率等值线范围与A样方红砖相近，约为30 cm。由此可以看出粉煤灰充填的土壤含水率等值线范围大于红砖和瓷砖，Algeo等[20]的研究成果可作为佐证。张敬凯等[21]对采煤塌陷地带干旱期和丰水期的塌陷区和非塌陷区土壤含水率进行研究，与本实验B样方含水率呈现的规律较为一致，夏季雨水较为充沛，水分入渗比周边无充填基质土壤多，故从雷达测线上呈现距离样方越近含水率越高现象。该法对于不同充填基质的覆土土壤均有较好的适用性，虽含水率的分布有所不同，仍得到了较高的精度。

4 结论

本文通过GPR对不同充填基质各层土壤进行探测，得出以下结论：

(1)通过对雷达数据预处理，进行Hilbert变换，提取特征波段，求得振幅包络平均值，建立与含水率的响应关系。分别建立两个实验样方不同充填类型下的不同层次土壤与AEA-1的关系模型。通过烘干法的实测值与AEA-1预测值对比，RMSE均在2.01以下。结果表明AEA法适用于不同充填类型的复杂土壤环境，这为通过GPR快速无损监测新增耕地土壤含水率提供了一种可行途径。

(2)对比不同充填基质各层次AEA法和实测法所得到的结果发现，无论何种充填基质，AEA法反演的0～10 cm、10～20 cm含水率与实测值R2均大于0.80，但在20～30 cm和30～40 cm含水率与实测值R2均小于0.80，对于较深层土壤含水率反演建模仍需进一步研究。

(3)不同充填基质下各层土壤含水率均遵循随着深度增加，含水率上升这一趋势；以粉煤灰为充填基质的土壤含水率大于红砖和瓷砖。以样方中心沉降处为起点向外延伸的测线反演结果显示，以红砖、粉煤灰为充填基质的土壤含水率比周边无充填基质的土壤高。通过测线含水率反演分析等值线图，红砖为充填基质的土壤含水率等值线范围大约为40 cm，B样方以粉煤灰为充填基质的土壤含水率等值线范围大约为50 cm。