探地雷达波振幅与土壤含水量关系的数值模拟

李蕙君，钟若飞

1.首都师范大学资源环境与地理信息系统北京市重点实验室，北京10048 2.首都师范大学城市环境过程与数字模拟国家重点实验室培育基地，北京10048 3.首都师范大学三维信息获取与应用教育部重点实验室，北京100048

探地雷达波振幅与土壤含水量关系的数值模拟

李蕙君1，2，3，钟若飞1，2，3

1.首都师范大学资源环境与地理信息系统北京市重点实验室，北京10048 2.首都师范大学城市环境过程与数字模拟国家重点实验室培育基地，北京10048 3.首都师范大学三维信息获取与应用教育部重点实验室，北京100048

通过建立探地雷达波振幅与土壤介电常数之间的定量关系可反演土壤含水量.在采用时间域有限差分法建立土壤分层模型的基础上，模拟探地雷达电磁波在均匀分层土壤介质中的传播特性，分析土壤介电常数与雷达回波振幅之间的变化规律，证实了回波振幅随着介电常数的增大而减小，总结出介电常数与振幅之间的幂函数关系式.最后通过室内土壤模型对此关系式进行了验证，结果表明由此关系式反演出的土壤含水量值与土壤三参数仪测定的含水量值相一致.

含水量；探地雷达；数值模拟

传统的土壤水分测定方法均是定点测定，其中包括烘干法、中子水分仪测定法、时域反射仪（time domain reflectometry，TDR）测定法、伽玛射线法、电容器传感法等［1］.这些方法仅适用于单点测定土壤水分的情况，而在大区域内测定土壤水分时存在效率低和破坏土壤环境等问题.遥感技术的发展弥补了大区域范围土壤水分监测的空白，它在提供土壤水分日数据的同时也可以实现大区域水分的动态变化监测，这对宏观尺度上的土壤水分监测具有重要的意义.然而，该方法易受植被/雪覆盖、土壤温度、地形、地表粗糙度等的影响，且空间分辨率低，只能监测0.05 m深的土壤水分［2-3］.因此，针对中小尺度区域，传统方法和遥感方法都有一定的局限性，迫切需要一种新的技术手段来推动土壤水分测定的发展.

基于电磁波原理的探地雷达，其中心频率介于10 MHz～3 GHz之间，具有快速高效、探测深度深、分辨率高且对周围环境无破坏的优点，可对中小区域土壤水分进行快速无损探测.目前，利用探地雷达探测土壤水分的方法主要有地面波法、钻孔雷达法.文献［4-5］采用地面波法计算了土壤表层的含水量，并利用空间插值方法分析了一个60 m×60 m区域的水分变化；文献［6］对比分析钻孔雷达方法估计的土壤水分值与中子仪法测定的水分值，发现两种方法所得土壤体积含水量的均方根误差为0.03 m3·m-3.在地面波方法中，雷达天线的设置方式有单偏移距和多偏移距，不同设置方式采样深度不一，且在单偏移距中需要确定最佳天线距离，比如文献［7］利用450 MHz天线并分别采用此两种方式探测，得出最佳天线距离在1.5～2.0 m之间，但在实际应用中仍需标定.另外，这种方法首先需要确定地下反射层的深度，依据深度-时间计算电磁波传播速度，从而得到土壤介质的介电常数，然后反演出含水量，所以有学者通过预埋反射体来确定反射深度，或者应用于有明显反射层的区域［8］，而这无疑增加了探测难度也限制了其应用范围.钻孔雷达法是将发射天线与接收天线分别放置于两个相隔一定距离的钻孔中，记录电磁波从发射天线到接收天线的旅行时间，探测两个钻孔之间土壤介质的含水量，此法可探测土壤剖面的垂直水分变化，在雷达信号衰减相对较大的区域可使用.比如文献［9］在一个垃圾填埋地发现地面波法不能有效测出其含水量，钻孔法却表现得较好，但是操作起来耗时耗力.另外，文献［10］在比利时黄土带状分布区利用地面分离天线配置方式，采用全波形反演方法进行了土壤水分值的探测并生成连续变化的土壤水分剖面图.

文献［11-12］发现，探地雷达反射信号频率与含水量相关，即随着含水量增加，中心频率会变低.利用频率步进方法可精确提取出其频率漂移量，增加反演精度.研究发现电磁波在土壤介质中传播时，反射波振幅大小与土壤含水量有密切关系，当土壤含水量增加时，电磁波的振幅变小［13-14］.基于前人的定性研究，本文拟通过数值模拟方法，利用时间域有限差分法建立模型，模拟探地雷达电磁波在匀质分层土壤中的传播特性，分析电磁波振幅与土壤介电常数之间的定量关系，从而依据介电常数与含水量关系模型反演出土壤含水量.

1　相关理论基础

1.1时间域有限差分法

电磁波在介质中的传播服从麦克斯韦方程组，在介质的参数不随时间变化且各向同性的无源区域，麦克斯韦旋度方程可写成

式中，E为电场强度，H为磁场强度，ε为介电常数，σ为介质电导率，µ为磁导率，ρ为计算磁损耗的磁阻率.在直角坐标系中，此旋度方程可写成微分分量式，文献［15］对此微分分量式引入差分格式，则二维TM模式下的差分方程为

式中，Δx、Δy、Δt分别为x、y方向的网格步长以及时间步长，需满足稳定性条件

式中，对于非均匀区域，v选最大值.

1.2介电常数与含水量关系模型

土壤是固体、气体、液体组成的三相混合体，水的介电常数为80，空气的介电常数为1，一般土壤固态物质的介电常数在5～15之间.可以看出土壤的介电常数受其含水量多少的影响很大，说明介电常数可以反演出土壤含水量.

探地雷达反演土壤水分中应用最多的方法是Topp等人［16］提出的经验关系式.Topp等发现相对于土壤质地、土壤容重、土壤温度、含盐量等，土壤含水量对土壤介电常数的影响更大，于是提出了土壤介电常数ε与土壤体积含水量θ之间的关系式

2　GPRSIM数值模型构建及模拟分析

2.1数值模型构建及模拟

探地雷达在介质中传播时影响电磁波电性参数的因素主要有介质的介电常数、电导率和磁导率，此三者的变化会引起电磁波传播的速度、振幅等参数的变化.在探地雷达应用中，不考虑大部分介质的磁导率影响，因此利用探地雷达测定土壤含水量的方法适用于低损耗低盐土壤介质，此类土壤的磁导率一般近似为1［17］.在土壤学的相关研究中，对1 m土层以内的表层土壤水分研究最多，该层土壤受质地、降雨、蒸发、植物蒸腾等影响，土壤水分变化较大.水分含量的多少对电导率有很大的影响，且该影响随着土壤类型的不同而有所不同.比如同一含水量的砂土和黄土，砂土的电导率就低于黄土的电导率.因此，探地雷达电磁波在土壤介质中传播时有损耗和衰减，在电导率高的土壤介质中衰减更严重.

在自然状态下，土壤水分在垂直剖面上是呈梯度性变化的，而土壤含水量的变化会引起介电常数的梯度性变化，于是参考Goodman［18］介绍的模型，将数值模型设置为介电常数呈梯度性变化的土壤分层模型.模型长为2 m，深为1 m，每一层土壤都假设各向均匀，电导率为0.001 S/m，磁导率为1，介电常数为其唯一变化参数［17-18］.本文模拟土壤介电常数在6～24范围的变化，模型分别以1为梯度增加介电常数值，每一个梯度中再分别设置以0.1为梯度变化的介电常数，取平均值为这一梯度内的介电常数值.所建模型如图1所示，平均介电常数为6.5时的模型参数设置如图2所示.

图1　土壤分层模型Figure 1 Soil layered model

图2　模型介电常数设置Figure 2 Model settings of dielectric constant

利用探地雷达反演土壤含水量时，其使用频率带宽一般为10 MHz～1 GHz，而Topp关系式在高频域（500～1 000 MHz）轻质地土壤含水量与土壤介电常数的关系拟合中效果最好［19］，故本文选取天线频率为500 MHz.模型参数设置如下：模型介电常数如图2所示设置；脉冲响应设置为时间窗口40 ns、频率500 MHz、扫描/道512；天线方向响应采用软件默认设置.模拟结果如图3所示（ε=6.5），其中单道波形如图4所示（ε=6.5、11.5、16.5、21.5）.

图3　模拟结果（ε=6.5）Figure 3 Simulation results（ε=6.5）

2.2模拟结果分析

如图3所示，探地雷达电磁波在介电常数变化的界面发生反射，在探地雷达图像上出现明显的分层，能反应土壤分层现象，其中首尾两个层面在图像上显示较强，中间相对较弱.这是因为首尾界面处上下两层介质的介电常数差异较大，回波振幅幅度变大.由图4可以看出，当介电常数为6.5时，雷达回波振幅值均较大.随着介电常数增加，雷达回波振幅值明显减小，这是因为土壤介电常数与含水量相关性最大的缘故.随着含水量增加，介电常数变大，含水量对雷达波的吸收增强，雷达波衰减变强，反射回波信号就较弱.如图4（a）所示，中间相对平稳的4个振幅值是在介电常数6.3到6.7之间以0.1的等差值变化的，变化幅度小，差异不大.因此，本文以此4个振幅的平均值作为当前土壤介电常数下对应的振幅值，再分别计算每一介电常数下的振幅，结果如表1所示.

图4　单道波形图Figure 4 Single channel wave

利用SPSS软件进行回归分析，回归分析选用幂模型.R2为0.997时拟合效果最好，得到拟合曲线如图5所示，拟合公式为

3　室内模型验证

3.1模型建立与数据采集

为了验证模拟公式的有效性，另外建立一个室内实体模型（长1 m、宽25 cm、高20 cm），材料为塑料泡沫板，并用防水塑料布水平隔成3块区域，所用土壤为野外地表风干砂壤质土.此类土壤广泛分布于华北平原，土壤中砂粒、粉粒、黏粒比例适当，通气透水性好，土壤呈块状结构，有机质含量低，对探地雷达信号衰减较小.土壤装箱前，先将大块状土体碾碎成细小颗粒土，装入实体模型中并轻微压实.实验时土壤温度为12.9◦C，实体模型如图6所示.分别采用探地雷达和土壤三参数测量仪测定实体模型内土壤含水量.土壤三参数测量仪（简称三参仪）是基于时域反射仪法（TDR法）设计开发的仪器，可同时测定土壤的水分、温度、电导率3个参数［20］.TDR法是基于电磁波与介电常数理论的一种方法，通过电磁波在土壤中的传播速度来计算土壤介电常数，从而依据介电常数-含水量模型反演土壤含水量.TDR法可快速测定体积含水量，且无需校正，可达到误差为1%的精度［1］.

表1　介电常数ε与振幅数据ATable 1 Dielectric constant ε and amplitude data A

图5　土壤介电常数与振幅的关系Figure 5 Relationship between dielectric constant and amplitude

图6　室内模型Figure 6 Indoor model

在实验过程中，由于探地雷达信号存在边界效应，分别在0.20 m、0.25 m、0.30 m、0.35 m设置4个控制点，用土壤三参数仪测定这4个控制点处的含水量并记录；然后采用500 MHz屏蔽天线探地雷达进行测量，见图7.为了设置含水量的变化梯度，在第1次测定后对土壤试验进行人工加水，隔日进行第2次相同土壤水分测定实验.

图7　测量示意图Figure 7 Measurement schematic

3.2数据处理

对探地雷达图像进行预处理，包括去直流漂移、时间归零校正、能量增益、背景去除、带通滤波、平均道，然后提取出4个控制点处的振幅值.本文所使用的模拟软件使用的是delta函数，故软件记录的振幅强度值均在-1到1之间，这与探地雷达数据采集软件记录的振幅强度值有差异.因此，需对提取出的雷达波振幅进行幅度变化（缩小10 000倍），从而模拟软件中的振幅值相适应；然后按照式（1）计算出每个控制点处的土壤介电常数值，再由Topp公式反演出控制点处灌水前后的含水量，所得结果如表2所示.

表2　灌水前后三参仪与GPR分别测定的含水量值Table 2 Soil water content measured by wet sensor and ground penetrating radar before and after migration

对比分析三参仪测定土壤含水量值与GPR反演的土壤含水量值，如图8所示.灌水前试样土壤含水量约为0.13 m3·m-3，灌水后土壤含水量增加至0.20 m3·m-3左右.两种方法测定的数据在控制点处含水量的变化趋势一致，利用探地雷达数据由式（1）反演出的含水量值与三参仪测定的含水量值吻合性较好，其均方根误差为0.014 m3·m-3.

图8　灌水前后三参仪与GPR测定含水量对比Figure 8 Comparison of soil water content measured by wet senor and ground penetrating radar before and after irrigation

4　结论

土壤介电常数主要受其含水量的影响，于是可通过求取土壤介电常数来反演出土壤含水量.本文通过构建土壤分层模型，模拟探地雷达电磁波在土壤介体中的传播特性，证实雷达回波振幅值随着土壤介电常数的增加而减小，这与已有研究结果相一致［14-15］.然后回归分析得到土壤介电常数值与雷达回波振幅值之间的幂函数关系式.最后建立室内实体模型，并与土壤三参数仪探测的数据进行对比来验证其有效性.结果显示，由式（1）反演出的含水量值与土壤三参数仪测定的含水量值吻合得较好，均方根误差为0.014 m3·m-3，说明此研究得出的土壤介电常数与振幅的关系式有一定的合理性.

需要说明的是：

1）在本文实验是在结构单一质地均匀的土壤模型中模拟，忽略了实际情况下土壤的其他属性特征，比如土壤电导率受土壤含水量的影响很大.随着含水量增加，电导率会变大，而探地雷达电磁波在高电导率介质中衰减.因此，在实际应用中应考虑电倒土壤其他物理性质的影响.

2）在本次实验中，模型箱体本身对雷达波信号有影响，导致探测过程中存在边界效应，且土壤试样准备不理想，导致得到的探地雷达数据杂波信号较多.振幅提取时，是在预先试验几次后选择了干扰较小波形较好的道，所以结果较为理想.那么，在实际应用过程中，探地雷达数据解译中本身存在的不确定性和主观性，将增加振幅精确提取的难度，从而影响最终的反演结果.

3）在自然状态下，土壤结构复杂，组成成分多样，对探地雷达电磁波有衰减吸收反射，可能导致反射振幅失真.探地雷达数据处理时振幅的恢复和提取精度会直接影响到最后含水量反演的精度.文中使用的介电常数和含水量关系模型是TOPP经验关系模型，在实际应用中，公式中的系数也需要依据试验地的土壤物理性质进行校正.

因此，本文所得只是一个初步结果，后期将以更多实验数据来完善此方法.在下一步的实验研究中，将考虑更多的影响因子，诸如不同土壤质地、不同土壤孔隙度、不同天线频率等.在实际应用中，需依据自然区域的土壤状况对式（1）进行改进.

［1］Pariva D，Ashi Q，Ruchi B.A review of the methods available for estimating soil moisture and its implications for water resource management［J］.Journal of Hydrology，2012，458：110-117.［2］Wigneron J P，CalveT J C，Pellarin T.Retrieving near-surface soil moisture from microwave radiometric observations：current status and future plans［J］.Remote Sensing of Environment，2003，85：489-506.

［3］Jackson T J，Schmugge J，Engman E T.Remote sensing applications to hydrology：soil moisture［J］.Hydrological Sciences Journal，1996，41：517-530.

［4］Huiiman J A，Sperl C，BouTen W，VersTraTen J M.Soil water content measurements at diferent scales accuracy of time domain reflectometry and ground-penetrating radar［J］.Journal of Hydrology，2001，245：48-58.

［5］Huisman J A，Snepvangers J J J C，BouTen W，Heuvelink G B M.Mapping spatial variation in surface soil water content comparison of ground-penetrating radar and time domain reectometry［J］.Journal of Hydrology，2002，269：194-207.

［6］Alumbaugh D，Chang P Y.Estimating moisture contents in the vadose zone using crossborehole ground penetratin radar：a study of accuracy and repeatability［J］.Water Resources Research，2002，38（12）：45_1-45_2.

［7］Galagedara L W，Parkin G W，Redman J D.Field studies of the GPR ground wave method for estimating soil water content during irrigation and drainage［J］.Journal of Hydrology，2005，301：182-197.

［8］LunT I A，Hubbard S S，Rubin Y.Soil moisture content estimation using ground-penetrating radar reflection data［J］.Journal of Hydrology，2005，307：254-269.

［9］April Y，Richard G Z，Edward A M.Estimating water content in an active landfll with the aid of GPR［J］.Waste Management，2013，33：2015-2028.

［10］MineT J，Wahyudi A，BogaerT P.Mapping shallow soil moisture profles at the feld scale using full-waveform inversion of ground penetrating radar data［J］.Geoderma，2011，161：225-237

［11］BenedeTTo A.Water content evaluation in unsaturated soil using GPR signal analysis in the frequency domain［J］.Journal of Applied Geophysics，2010，71：26-35.

［12］Nicolaescu I，van Genderen P.Performances of a stepped-frequency continuous-wave ground penetrating radar［J］.Journal of Applied Geophysics，2012，82：59-67.

［13］Schmalz B，LennarTz B，WachsmuTh D.Analyses of soil water content variations and GPR attribute distributions［J］.Journal of Hydrology，2002，267：217-226.

［14］Chanzy A，Tarussov A，Judge A，Bonn F.Soil water content determination using a digital ground-penetrating radar［J］.Soil Science Society of American Journal，1996，60：1318-1326.

［15］Yee K S.Numerical solution of initial boundary walue problems involving Maxwell's eqation in isotropic media［J］.IEEE Transactions on Antennas and Propagation，1966，14（3）：302-307.

［16］Topp G C，Davis J L，Annan A P.Electromagnetic determination of soil water content：measurements in coaxial transmissions lines［J］.Water Resource Research，1980，16：574-582.

［17］朱安宁，吉丽青，张佳宝.基于探地雷达的土壤水分测定方法研究进展［J］.中国生态农业学报，2009，17（5）：1039-1044. Zhu A N，Ji L Q，Zhang J B.Research progress on soil moisture measurement via groundpenetrating radar［J］.Chinese Journal of Eco-agriculture，2009，17（5）：1039-1044.（in Chinese）

［18］Dean G.Ground-penetrating radar simulation in engineering and archaeology［J］.Geophysics，1994，2（59）：224-232.

［19］Robinson D A，Jones S B，WraiTh J M.A review of advances in dielectric and electric conductivity measurements using time domain reflectometry［J］.Vadose Zone Journal，2003，2：444-475

［20］赵春江，王成，侯瑞峰.土壤三参数测量方法［J］.现代科学仪器，2007，5：100-104. Zhao C J，Wang C，Hou R F.Preliminary study of the three soil parameters measurement methods［J］.Modern Scientifc Instruments，2007，5：100-104.（in Chinese）

（编辑：管玉娟）

Numerical Study on the Relationship between Amplitudes of Ground Penetrating Radar
Wave and Water Content in Soil

LI Hui-jun1，2，3，ZHONG Ruo-fei1，2，3
1.State Key Laboratory Incubation Base of Urban Environmental Processes and Digital Simulation，Capital Normal University，Beijing 100048，China 2.Beijing Key Laboratory of Resource Environment and Geographic Information System，Capital Normal University，Beijing 100048，China 3.Key Laboratory of 3D Information Acquisition and Application，Ministry of Education，Capital Normal University，Beijing 100048，China

Water content in soil can be derived from the relation between changes in ground penetrating radar signals and the soil's dielectric constant.Using the fnite diference time domain（FDTD）method，a model of layered soil is established to simulate ground penetrating radar electromagnetic waves in the soil.The transmission characteristics of layered soil are analyzed using the model.The results indicate that a functional relation exists between the soil dielectric constant and the amplitude，i.e.，the wave amplitude decreases with the increased electric constant.With an indoor soil model，the obtained relationship is verifed：the computed water content agrees with the wet sensor data.

water content in soil，ground penetrating radar，numerical simulation

P631.3，S152.7

0255-8297（2015）01-0041-09

10.3969/j.issn.0255-8297.2015.01.005

2014-12-02；

2014-12-28

国家科技支撑计划课题基金（No.2012BAH34B01）；国家自然科学基金（No.41371434）资助

钟若飞，教授，博导，研究方向：车载激光测量，E-mail：zrfsss@163.com