基于探地雷达的土体构型无损探测方法研究

2022-02-09姚喜军陈晓东

干旱区地理(汉文版) 2022年6期

吴全，姚喜军，陈晓东，赵敏，赵欢，云浩

（1.内蒙古自治区国土空间规划院，内蒙古呼和浩特 010000；2.内蒙古达智能源科技有限公司，内蒙古呼和浩特 010000）

土壤是地球上所有生态系统的自然支撑和生物屏障。近年来随着工业化、城市化进程的加剧，人与土地之间的矛盾也日益突出，土壤质量和可持续发展已成为社会经济发展和土壤学、农学、环境科学等研究的热点［1］。土壤质量检测方法通常有宏观和微观两种，宏观即卫星遥感探测法，是利用卫星遥感技术探测土壤表层属性，对土壤表层含水率、全氮和全磷含量、盐渍化、重金属含量等参量进行检测［2-7］，能以不同时空尺度不断地观察地表情况，提供地表特征信息［8］，适用于空间覆盖范围较广的土壤表层信息监测，但监测数据的空间分辨率往往偏低；微观即点测法，是通过剖面挖掘、打土钻、取样方、测土样等方式获得有限的点测数据［9］，来推断整个面或三维空间的土壤质量，测量准确，但仅适用空间覆盖范围偏小的检测。从空间上来讲，遥感探测法与点测法不够全面，亟需在中观上建立一种客观、快速、准确的土壤质量检测方法。

探地雷达是一种利用高频电磁波的反射探测目的体即地质现象的物探方法，亦称地质雷达，具有高效、快速、连续、无损、低成本和高分辨率成像等特点，已成为公路路面、机场道面、隧道、水坝和堤防等层状结构基础工程设施无损检测技术的重要组成部分［10-12］。近年来，随着探地雷达技术的快速发展逐渐应用到土壤检测领域中，成为宏观探测与微观探测之间尺度转化的纽带。如吴志远等［13］利用探地雷达早期信号振幅包络法平均值方法和时域反射技术分别探测降水前后野外农田表层的土壤含水率，结果表明利用探地雷达波早期信号振幅包络值能够获得与时域反射技术实测精度相近的黏性土壤含水率；侯晓冬等［14］将探地雷达用于土壤污染物探测中，以土壤介电常数为纽带，通过神经网络技术，建立介电常数、雷达回波信号和污染土壤多个参数之间的联系，能够实现污染土壤中多种污染物含量的检测；也有部分学者将探地雷达用于土壤深度和分布的探测，如赵艳玲等［15］使用探地雷达在室外进行复垦土壤层次无损探测的模型实验，证明探地雷达能够有效探测土层厚度和结构。以上研究，第一从探地雷达回波信号着手，建立土壤与雷达波形数据的关系；第二研究土壤的单一属性，如含水率、层厚、分布、污染物含量等。本文提出了一种使用探地雷达快速测量土壤土体构型的无损探测方法。该方法分别以探地雷达波形及其图像为研究对象，结合探地雷达波形数据分析算法、土壤介电常数推演、基于主成分的图像白噪声估计、支持向量机等技术，研究土壤探地雷达波形与土壤物性之间的关系，以探测土壤剖面的三大主要属性（层次、层厚、土质），综合地对土壤土体构型进行研判，并编制了一套快速识别土壤土体构型的信息系统。野外探测试验表明该方法实现了对试验基地区域内土壤土体构型快速地、无损地探测。

1 研究区概况

试验选择2 个典型区域，分别位于内蒙古自治区标准耕作制度分区中Ⅱ5河套平原区的土默特左旗和Ⅱ2 黄土丘陵区的达拉特旗（图1）。土默特左旗位于内蒙古自治区中南部、黄河北岸，属于典型大陆性半干旱季风气候。降水年际变化大，且年内分配不均，7—9 月降水量占全年降水量的68.8%，年平均降水量约350 mm。试验采样区位于土默特左旗（109°54′44.3″E，40°19′46.0″N）附近，分别设置3 个采样点，采集的典型土壤剖面构型为壤土/黏土（表示上层为壤土下层为黏土，后文描述的土体构型均以此方式表示）、砂壤土/壤土。达拉特旗地处内蒙古自治区鄂尔多斯高原北端，黄河中游内蒙古段“几字湾”南岸［16］，属于典型的半干旱大陆性气候，受季风环流影响，达拉特旗春季降水少，夏季降水集中，且容易出现大到暴雨天气，该地区降水季节分配不均，旱涝灾害频繁，年平均降水量约300 mm。试验采样区位于达拉特旗的库布其沙漠（109°53′11.4″E，40°16′41.7″N）附近，设3 个采样点，采集的典型土壤剖面构型为砂土/砂壤土、通体砂。

2 数据与方法

2.1 数据采集

选取壤土/黏土、砂壤土/壤土、砂土/砂壤土、通体砂4 种具有明显差异的土体构型作为试验对象，在无作物覆盖的耕地开展雷达测线布置和土壤剖面点测双盲测试。每种土体构型分别布置了10 条长10 m的平行测线。

点测采样：在每条测线的中点进行土壤剖面挖掘，土壤剖面描述主要依据土壤颜色、土壤松紧度、粒径分布等特征判定，测量各质地层次的深度和厚度，利用质地层次排列顺序及层次厚度初步判定土体构型，取各层土壤样品，进行实验室化验测定土壤指标参数。

探地雷达采样：设备采用青岛中电众益智能科技发展有限公司的GER-10 型脉冲式探地雷达（图2），因试验区的耕作土层厚度约为1 m，结合土壤介电特性受土壤质地影响的最小范围，选取900 MHz的发射天线，雷达具体参数及指标见表1。

图2 GER-10型脉冲式探地雷达Fig.2 GER-10 pulsed ground penetrating radar

表1 GER-10型脉冲式探地雷达技术参数及指标Tab.1 Technical parameters and indexes of GER-10 pulsed Ground Penetrating Radar

探地雷达沿平行测线往返探测数据，一条测线探测4 次。现场测量开始前首先使用Pereometer V.7型介电常数仪、SU-LB型土壤水份仪对样地进行探测，根据探测值对雷达的采集参数进行设定。当遇到坑洼地区时，截断采集数据避免空气波对图像的影响。一条测线探测1次约生成4～6组雷达波形图。参数设定的内容包括时间窗口大小、扫描样点数、每秒扫描数、A/D转换位数、增益点数等内容［17］。参数设置的是否合理影响到记录数据的质量，至关重要。

（1）探测深度与时窗长度

探测深度的选择尤为重要，选得太小会丢掉重要数据，选得太大会降低垂向分辨率。一般选取探测深度（H）为目标深度的1.5倍。文中探测的土壤深度为1 m，探测深度设置为1.5 m。根据H和介电常数（ε）确定采样时窗长度（Range，单位为ns），公式如下：

以达拉特旗的砂土/砂壤土为例，该区域为农耕土，ε为11，Range 应选40 ns，时窗选择略有富余。

（2）A/D采样分辨率

雷达的A/D 转换有8 Bit、16 Bit、24 Bit，由于探测深度在1.5～5 m之间，采集速度平缓，选择16 Bit。

（3）采集模式及扫描样点数

采集模式分点测模式、距离模式、时间模式，试验时选择距离模式，该模式下雷达按照相等的距离间隔采集雷达数据。扫描样点数Samples 有128、256、512、1024、2048 可选，对于雷达中心频率为Fa、时窗长度为Range，Samples 应满足下列关系：

Samples ≥10-8×Range×Fa （2）

为保证高的垂向分辨，在容许的情况下Samples尽量选的大些。文中雷达主频率Fa 为900 MHz，Range 为40 ns，每一道扫描样点数大于360，选择接近值为512。

（4）扫描速率

即系统每秒采集的数据道数，GER-10 有32、64、128、256、512 共计5 种，降低扫描速率可以提高系统的信噪比。当选择距离模式采集数据时，扫描速率自动设置在当前采样点数下的最高扫描速率。如本文主机连接900 MHz 天线，扫描样点数设置为512，则扫描速率自动选择为512 Traces·s-1（道每秒），测距轮的间隔设置为2 cm，那么雷达的最大工作速度为36.86 km·h-1，如果天线的移动速度超过该速度时采集数据水平距离的误差将增大。

2.2 研究方法

2.2.1 图像预处理地下所有物体会产生回波信号，介质的不均匀也会产生回波信号［18］。土壤粒径大小、压实度、含水程度、地面坡度等几何特性的不均匀性，使得探地雷达接收的回波信号包含多种随机噪声和杂波干扰。本文利用均值法、带通滤波法、中值滤波法、指数增益调节等信号处理算法对雷达波数据进行去除背景影响、高频杂波、来自地面发射的直达波、信号增强等预处理，提高信号的信噪比［19］。

2.2.2 基于包络检波的土壤分层信息识别探地雷达发射的电磁波信号遇到不同质地的土壤时会产生不同波形幅度与相位的回波，波形的正负峰分别以白色和黑色标识、以灰阶或彩色标识。这样，同相轴或等灰度、等色线即可形象的表征地下反射面。反射脉冲波形的明显程度，是对探地雷达图像进行地质解释的重要依据［20］。

为抑制干扰噪声的影响，采用基于Hilbert 的包络检波法从探地雷达回波中提取包络信号［21-24］，分析其瞬时相位来确定土壤分界面的位置。探地雷达的回波信号x(t)进行Hilbert变换，公式如下：

式中：X̂(t)为Hilbert 变换后所得结果；x(t)为探地雷达回波信号；t为时间；*表示卷积运算。构建复数域的解析信号x͂(t)，公式如下：

然后对ϕ(t)进行高阶最小二乘拟合得到拟合曲线。经反复试验，当阶数取20时拟合曲线的凹凸拐点与土壤剖面分界线的实际位置几乎完全拟合。因此将该拟合曲线上的凸凹拐点处作为土壤分界面的时域位置。

2.2.3 土壤介电常数与土层厚度计算采用探地雷达回波振幅反演介电常数方法［25-27］，确定结构层的介电常数与回波的波幅关系，公式如下：

εi=εi-1[(1+R)/(1-R)]2i=1,2,…,N（7）

式中：εi为下层介质的介电常数；εi-1为上层介质的介电常数；R为反射系数；N为表层土壤的总层数。R的计算公式如下：式中：A(t)为反射波幅；Am为全反射波幅数值，可以通过金属板试验获得。利用探地雷达测量时，顶层介质为空气，第1 层介质为第1 层土壤。因此可由顶层的空气介电常数，推演出第1层土壤介电常数，依次可推出其他土层的介电常数。而对于介电常数为εi的土壤介质来说，电磁波在土壤中传播的速度（vi）可由公式（9）推导：

式中：c为电磁波在真空中的传播速度，取3×108m·s-1。各土层厚度计算公式为：

式中：ΔHi为第i层土壤厚度；Δti为电磁波在第i层土壤中的双程时延。土层厚度可根据探地雷达在土层中回波时延和电磁波传播速度计算获得。

2.2.4 土质识别方法研究（1）基于主成分分析的土壤图像噪声估计。土壤的粒径、密度、孔隙度、渗透性、磁性等属性不同，生成的回波信号图像包含噪声不同，因此通过估计探地雷达回波图像噪声可间接反映所测土质的属性。陈晓东等［28］对土壤的探地雷达图像进行基于回归分析和主成分分析（Principal component analysis，PCA）的真实噪声均方差的估计，提出一种利用图像噪声估计分析土壤含砂量的方法，并指出回波图像噪声数值大小与土壤含砂量成正比关系。本文利用该方法，来计算土壤的含砂量。

探地雷达所采集到的含噪图像（f）可看作是由不含噪声的真实图像（r）和加性高斯白噪声信号（n）组成，其表达式为:

式中：n为均值为零、方差为σ̂n2 的加性高斯白噪声信号。依据文献［29］的分析，将图像数据f、r变换到PCA域中得到：

式中：λmin(·)为求取协方差矩阵的最小特征值；Cf为f的协方差矩阵；Cr为r的协方差矩阵；σ̂n2 为加性高斯白噪声信号的方差。由于真实图像的细节信息不丰富且具有信息冗余特性，r的数据主要分布在有限主分量上，研究时可以忽略，即存在如下公式：

由此，可以用高斯噪声的方差来估算探地雷达回波图像的方差值，即公式（12）简化为如下公式：

对每一层土质图像都进行噪声估计，将计算的方差值作为该层土质识别的特征参数之一，其数值大小可间接反映土壤含砂量大小。

（2）土壤含水率。在降水或灌溉与土层深度相同的条件下，不同土质因渗透性能存在差异其土壤含水率不同，因此将土壤含水率作为土质识别的另一特征参数。本文使用经典的Topp 多项式［30］来计算εi的第i层土壤含水率（θ），εi的计算方法参见公式（7），θ的计算公式如下：

（3）支持向量机的土质识别。以每类土质图像所计算的σ̂n2、θ、平均灰度、三阶中心矩、一维熵（H）（表示图像中灰度分布聚集特征所包含的信息量）等组成样本特征向量并建立对应的土质特征库，将土质特征库数据带入支持向量机进行训练并识别各土层的土质。

3 结果与分析

3.1 室内试验分析

在分析研究阶段，2019 年8 月设计并建设了土壤分析试验仓。针对土体构型检测设计了5层土壤分层试验，第1～3层和第5层分别铺设粒径为大于2 cm、等于2 cm、1～2 cm、1 cm 大小的砂壤土，第4 层铺设粒径0.5～1 cm 黏土层。采集的原始图像、预处理图像与分析结果图像如图3 所示，图中纵轴为土壤深度，水平纵轴为雷达扫描道数，下文以此图像分析土壤的土体构型。

将图3b中每道雷达数据沿行径方向取均值，并进行归一化和去直流分量处理，得到如图4 所示黑色土层轮廓线，对其进行Hilbert变换后得到红色上包络线、绿色下包络线。

图4中包络线的瞬时相位信号最小二乘拟合后得到瞬时相位拟合曲线（图5），利用微积分获取曲线凸凹拐点（图5 中绿色圆圈），凸凹拐点所在位置即为各土层分界点，最终求得土壤分层图像（图3c）。

图3 探地雷达图像Fig.3 Image of ground penetrating radar

图4 包络分析曲线图Fig.4 Abacuses of envelope analysis

图5 瞬时相位拟合曲线及凸凹拐点图Fig.5 Instantaneous phase fitting curve and convex and concave inflection point plot

表2为室内土壤土体构型试验数据结果。其中包括采用传统仪器测量的结果，简记为测量值，通过本方法计算出的结果简记为计算值。从表2可以看出，本文所述方法与实测数据相比，介电常数的相对误差最小为6.16%，最大为9.57%，层厚的相对误差最小为7.58%，最大为9.13%；土壤含水率的相对误差最小为8.46%，最大为9.69%。另外，土层的雷达回波图像噪声方差越大，其对应的土壤含砂量越大，且各土层的图像噪声方差差异较大，土壤剖面构型总体识别也较为准确。

表2 土壤属性检测结果Tab.2 Results of soil property testing

3.2 土体构型知识库的创建与系统搭建

建立了涵盖土壤的地域信息（如所属旗县、坐标范围、气候信息）、属性指标（土质、含水率、容重、实际剖面构型属性）、雷达探测信息（探地雷达设置参数、探测时间、探线走向、探测波形图）、图像多特征融合信息（噪声方差、平均灰度、三阶中心距、一维熵等）的土体构型特征知识库。建立的基于土壤容重为1.4～1.6 t·m-3，土壤含水率为4%～7%的典型沙地、砂壤土/壤土、壤土/黏土、砂土/砂壤土7 种土体构型知识库样，是智能分析土体构型的基础支撑。

搭建基于探地雷达的土体构型快速检测系统，具有雷达波形预处理、数据分析、智能检测和实验报表自动生成等功能。波形预处理包括取值去直达波、中值滤波、高斯滤波、增益调整等；数据分析主要包括数据导入及筛选、分析检测、指标计算、结果呈现、辅助功能；智能检测从探地雷达数据中分析土壤的土体构型，并给出各项指标；实验报表包括当前实验分析数据报表和历史分析报表。

3.3 野外试验分析

2020年9月在土默特左旗和达拉特旗采集土壤土体构型（壤土/黏土、砂土/砂壤土）进行野外双盲试验。图6a与图6b为壤土/黏土的分层结果图像和剖面图，图6c 与图6d 为砂土/砂壤土的分层结果图像和剖面图。试验证明，本方法计算的分层层数、土质与实际测试的结果相同，壤土/黏土的层间厚度误差在2～6 cm 之间，砂土/砂壤土的层间层厚误差在1～6 cm之间，被检测土壤的土体构型与实际完全相符，准确率约为94%。具体的土体构型试验结果见表3。

表3 土体构型试验结果表Tab.3 Soil configuration test results of Tumed Left Banner

图6 探测结果对比图像Fig.6 Comparison image of probe results

4 讨论

采用本方法野外探测的4 种土体构型，土壤层次均为4 层，除障碍层外每层土壤厚度大约20～30 cm，相邻层的土壤质地具有相似性，土壤层次探测准确；由于探地雷达采用的介电常数为仪器探测均值，输出的探测深度略有偏差，相比人工检测土层厚度误差在6 cm以内；每层土壤图像的噪声估计差异明显，土质检测结果与人工基本一致，总体土体构型识别较为准确。

本研究从室内控制实验出发，分析土壤物性与探地雷达波形及其图像间定性与定量的关系，提出利用探地雷达波形的纵向梯度信息分析土层层次，以雷达回波振幅反演土壤介电常数进而推到每层土壤层厚，以图像噪声方差来估计土壤含砂量的理论，然后以野外双盲探测的方法验证了该方法在土壤土体构型检测中的适应性。第1次提出把探地雷达波形图像作为研究内容，从图像分析的角度来识别土壤质地，丰富并拓展了此类研究项目的数据源和解决问题的思路；创新性地提出将传统人工检测仪表与探地雷达相结合进行土壤物性探测，以土壤的层次、厚度、土质为切入点综合地分析土体构型的工程思想。与以往发表的论著不同，本文结合探地雷达波形图像分析方法，从综合的角度给出了一种检测土壤土体构型的方法，并通过室内控制实验和野外实验进行了验证，具有较大的工程意义；相比人工挖剖面分析方法，既不损坏原有土壤结构，又极大地提高了检测速度和范围。

研究发现，由于使用的探地雷达收发信号是脉冲电磁信号，雷达回波数据受地表植被覆盖物、地下植物根系等因素影响较大，其距离地面的真空距离大于一定范围时，将影响采集的信号精度，为此实际探测时应避开植被覆盖时间或地点；探地雷达回波信号受水的干扰很大，特别是处于水饱和状态时土壤的其他特性在图像中均被掩盖，野外实际测量应避开雨天并在地表比较干燥时进行测量，参考地表土壤黏度、水分等情况合理设置雷达探测参数。同时，发现土壤障碍层含水率相对较高，依据探地雷达回波信号与水的关系，可以用这一特性探测土壤障碍层或者水管管道漏水点检测等。

本文已经证明在典型区域使用此方法可以快速、无损、准确地探测土壤土体构型，但我国区域辽阔，土壤土体构型差异较大，雷达发射信号在这些土壤中的穿透能力、受干扰情况也不同，具体扩展到其他区域的试验需下一步工作中进行深入研究。