土壤含水率与GPR信号属性关系数值模拟研究

伊天宇

(中国海洋大学海洋地球科学学院，山东青岛266100)

实际生活中，应用时域反射法检测物质含水量已较为普遍。目前，探地雷达正演模拟通常假设目标体及其周围介质是均匀的，然而物质天然存在着不均匀性，雷达波在不均匀介质中的传播特性会发生改变。特别是，介质中的水会影响整体介质的介电常数。因此，应用探地雷达技术探测物质含水情况具有良好的地球物理前提[1]。

在物理模拟方面，目前已有利用GPR探测土壤含水率的实例。马福建等[2]利用GPR探测人工模型得到了反射系数与含水率之间的定量表达式。吴信民[3]提出相位移时间的概念，并建立了土体含水率与相位移时间关系的理论模型。而在数值模拟方面，方慧[1]定性的研究了介质含水率与雷达信号的关系，其中用到了雷达振幅属性。总的来说，运用GPR探测土壤含水率定性研究居多，而定量研究相对较少。

基于此，本文采用随机过程理论建立土壤含水平稳随机介质模型，并根据雷达振幅和含水率差值来讨论雷达信号属性与土壤含水率之间的关系。振幅是雷达信号中绝对意义的量，含水率差值是相对意义上的量。只要能够得到浅层土壤含水率，即可通过关系式推导下一层土壤含水率。

1 基本理论

1.1 水的电性特征

在时间域中，电磁波在介质中的传播速度v由因子k决定[1]：

实验表明介电常数还随频率变化，用Debye公式来描述变化规律：

式中：εS、ε∞——直流和极高频下的介质介电常数；

frel——弛豫频率。

Debye模型表明：在低频和高频状态下，介质的介电常数近似为实数，频散不明显，能量耗散较小。

1.2 Topp公式

本文主要基于Topp[4]通过实验给出的土壤介电常数与土壤含水率之间的关系来建立模型：

式中：εb——土壤整体介电常数；

θ——土壤含水率。

2 方法

2.1 建模思路

含水土壤中，水的存在将会影响土壤的电学特性，即会影响土壤的介电常数和电导率。由于水随机分布在土壤中，该固液-双相介质属于随机介质，所以根据式及随机介质理论建立含水土壤平稳随机介质模型。具体实现如下：

（1）含水土壤体规模：0.5m×1m×0.5m，如图1所示；

（2）模型分层如图1所示；

（3）每层赋予背景含水率，共设计8种含水率并代入式。计算，如表1所示；

表1 土壤含水率与相对介电常数对应数据

（4）实现层内含水率随机变化[5]，逐层求取y方向上每层网格内的含水率：

式中：f(x,z)——含水率；

f0——每层的背景含水率；

δf(x,z)——随机扰动的含水率；

σ2(x,z)——具有一定均值和方差的二阶平稳的空间随机过程。

（5）产生随机变化的介电常数体，将模型每个位置的含水率代入式。即可求得该位置的介电常数，从而构建含水土壤体介电模型。

2.2 含水土壤体的算法实现

一般情况下随机过程不存在傅里叶变换，但可从自相关函数的功率谱来描述，其算法原理如下：

（1）选择自相关函数φ(x,z)，这里选用高斯型椭圆自相关函数[6]，其表达式为：

（2）对自相关函数做二维傅氏变换，得到φ(kx+kz)；

（3）用随机数发生器产生在[0,2π]服从均匀分布的二维随机场Θ(kx,kz)；

（4）计算随机功率谱函数：

（5）计算R(kx,kz)的二维反傅氏变换得到=(kx,kz)；

（6） 计 算 均 值μ=E[,z)]和 方 差d2=E[((x,z)-μ)2]并规范化，即使μ=0,d2≈0.03；

（7）再由式得到模型。y方向单层网格点随机扰动建立结果如图2（a）所示，含水率结果如图2（b）所示（第一层含水率15.58%，第二层含水率20.23%），相对介电常数结果如图2（c）所示。

2.3 正演模拟参数设置

探测方式为剖面法，天线间距为20cm，以5cm为探测间距沿测线移动，测线选取为模型中线。采用三维时域有限差分法对含水土体随机介质进行正演模拟[7]。

为满足稳定性条件，通常剖分单元每一个边的长度不得大于发射天线主频波长的10，并且要求时间步长的设计满足下式：

式中：Δx、Δy和Δz——剖分单元的长度；

c——光速。

考虑天线主频，选取空间网格dx=0.001m，时间步长dt由式计算来得到，以满足稳定性关系。

震源采用500MHz雷克子波。

2.4 数据处理

2.4.1 最优频率选取

Topp公式在高频域(500～1000MHz)轻质土壤含水量与土壤介电常数的关系拟合中效果最好[3]，故本文选取天线频率为500MHz，如图3（截取部分图像）所示。

2.4.2 层速度求取

依次将各含水率（共8种）设置在第一层，第二层设置为与第一层含水率不同。根据第一层与第二层界面反射波时至与震源时间差值t（即双程走时）及图像上每层深度h与雷达天线距d来求取层速度：

因而有对应关系见表2。

2.4.3 振幅求取

本文研究的振幅处理方式为Hilbert变换[8]，该变换是研究瞬时包络、瞬时相位及瞬时频率的重要手段。

对于任意连续实函数x(t)的Hilbert变换定义为：

其中H表示Hilbert变换，∗表示卷积。x(t)与其Hilbert变换是正交的。Hilbert变换相当于通过一个冲击响应为的线性网络，即Hilbert滤波器。

对于x(t)，若其Hilbert变换为R(t)，则有瞬时包络：

图4为其道Ex分量包络图。

表2 土壤含水率与层速度对应数据

3 结果与分析

3.1 土壤含水率与层速度

由前面的表2可知，层速度对土壤含水率非常敏感，即介电常数依旧是对土壤含水率的指示性指标，利用MATLAB将表2数据进行回归分析可得层速度与土壤含水率的关系模型，可以看出呈现明显的负相关线性关系，拟合曲线见图5。

回归分析选用线性模型，R2=0.9927，拟合非常好，表达式为：

式中：θ——含水率；

v——层速度，ms。

由图5可以看出土壤含水率与层速度存在线性模型关系，含水率越大，层速度越小。理论上，如果地下存在明显反射层，应用宽角法等测取表层土壤层速度，便可很方便准确的得到表层土壤的含水率。

3.2 土壤含水率差值与初至反射波振幅

分别将含水率设置在第一二层，保证下一层比上一层的含水率要大，类似于图3，可得到一二层分界面的反射波，再进行Hilbert变换求取振幅信息，这里选用Ex分量进行处理，处理得到的含水率差值与反射波振幅的数据，见表3。利用MATLAB软件进行回归分析，拟合曲线见图6。

表3 含水率差值与Ex振幅对应数据

回归分析选用线性模型，R2=0.9939，拟合很好，并由拟合参数得出表达式为：

式中：θ1、θ2——相邻2个土壤层的含水率，%；

A——振幅。

由图6可以看出含水率差值与反射振幅存在明显的正相关线性关系。随着含水率差值的增加，界面反射振幅也会跟着增加。反射振幅具有相对意义。由表3可知，反射振幅对含水率差值非常敏感。

在实际野外土壤中，水分是连续分布的，只要地下有明显的分界面，可以先利用宽角法等[9]对分界面上覆土壤进行速度测量，再代入式便可得到分界面上覆土壤的整体含水率。要继续递推往下的土壤含水率，需要从GPR雷达记录中提取反射序列，明确反射界面。由于重力的因素，含水率往往是越往深越大，故可依据反射振幅与含水率差值的关系逐层向下求取含水率差值，与表层含水率依次相加便可得到不同深度的含水率值。

4 结论

介质含水率变化对雷达信号的影响在多数情况下被认为是干扰，由于水具有特殊性质，使电磁波在介质的传播过程更加复杂。通过GPR信号属性研究土壤含水率，重点探讨了界面反射振幅与土壤含水率差值可能存在的关系。界面反射振幅与土壤含水率差值都与土壤含水相邻两层相关，都属于相对意义上的量，在实际中有更好的应用前景。

应用界面反射振幅与土壤含水率的关系式，只要得到表层的土壤含水率，依照反射系数分辨清楚反射振幅，并应用合适的振幅补偿，就能将层深引起的振幅衰减降低，从而突出含水率差值对反射振幅的影响，进而逐层换算含水率差值，由浅到深，逐层计算，进而得到深层土壤的含水层。

但在实际中，要想应用GPR信号属性探测土壤含水率，土壤中必须有一个明显的标志反射层，如果土壤呈近乎均匀分布，探测土壤含水率的工作也难以为继。正常情况下，土壤含水率越往深越大。但在正演模拟中，若下一层土壤含水率比上一层土壤含水率低，雷达信号便会出现极性反转的情况，这表明GPR的其他信号属性也有可能较好的描述土壤含水率。

总的来说，应用雷达信号属性测量土壤含水率是可行的，有着广阔的应用前景。其无损、非接触、高效的特性表现出其它检测介质含水率的传统技术所不具备的优势。