基于曲线拟合的探地雷达层状介质参数反演算法

刘仁杰 兰 天 杨小鹏 李奕璇

（1.北京理工大学信息与电子学院雷达技术研究所，北京 100081；2.北京理工大学重庆创新中心，重庆 401120）

1 引言

探地雷达（Ground Penetrating Radar，简称GPR）因其实时性、无损性、高效性被广泛应用于地质工程［1］、建筑监测［2］、地下管道电缆探测［3-4］、军事探测［5-6］、植物检测［7-8］等各个领域。在道路检测方面，探地雷达技术是利用高频电磁波对地下地质结构进行无损测量，相较于取芯法的高成本、低效率与破坏性，探地雷达技术逐渐成为地下结构的主流勘测方法。探地雷达发射的电磁波，在地下结构的传播过程中，遇到分层界面、地下裂缝等电特性差异的地方会发生反射、折射、绕射等现象［9］，在雷达接收端形成单次测量的A-scan 数据，通过改变收发间距或是移动收发位置，得到多次测量后的B-scan 数据。层状介质的参数反演方法主要是从B-scan 数据入手，基于共中心点（common middle point，简称CMP）的模型结构，对回波数据的电性参数进行提取与处理，反演出地下的介电参数与厚度。

Chen 与Colby 分别在文献［10］、［11］中介绍了电磁波在传播过程中的波形构成，包含了球形空气波、球形地波、横向波与消逝波四类，文献［10］强调了横向波在未来探地雷达技术的潜在应用，文献［11］给出了直耦波的传播方程。袁明德在文献［12］中简要概述了在B-scan 图像中由空气耦合波与地面滑行波推导得到光速与地下介质中波速的方法，即曲线斜率为波速的倒数，但是这类关于曲线拟合的方法仅局限于浅层的参数反演，不适用于多层介质下的参数反演。除了通过B-scan 曲线拟合的方式反演地下介质结构信息，常用的反演方法还包括常规的CMP方法［13-15］、速度谱［16-17］、反射系数的谱反演［18-19］等，Kao 在文献［14］中给出了共中心点的探地雷达模型，通过构造电磁波在分层介质中的传播路径，实现对介质参数的反演，但是论文方法对于相对介电常数的反演误差很大。Wang 在文献［16］中使用MIMO 体制的探地雷达收集CMP数据，利用速度谱方法，实现对地下介质速度的有效反演。

针对B-scan 图像中的反射波曲线，本文提出了一种基于曲线拟合的参数反演算法，通过对反射波的位置拟合以及电磁波的垂直入射模型，提取得到相对介电常数与厚度的关系曲线，获得等效相对介电常数与厚度值，经过Dix公式［20］，将均方根速度转变为平均层速度，反演得到地下介质的参数信息。

2 电磁波波形构成

当发射天线端置于介质层表面时，根据电磁波的传播形式，可以分为球形空气波、球形地波、横向波与消逝波四类。如图1 所示，波源S 处发射而出的电磁波以光速在空气中形成了A 所示的球形波，在介质中以较慢速度形成了B，两类球形波跨越介质分界面时分别形成了横向波C 与消逝波D。收发天线之间存在的空气耦合波主要由球形空气波构成，地表滑行波由球形地波与横向波一同构成。

在单个A-scan数据中往往会存在不同类型的波形，对于单层地下介质，接收端的接收波形包含有空气耦合波、地表滑行波、单次以及多次反射波，当地下存在多层介质时，根据分层界面上下介质相对介电常数的不同，会出现界面上下电磁波传播速度的差异，由费马定理的最短时间原理，即光线传播的路径是需时最少的路径，满足一定收发间距时，电磁波会在传播速度较快的一侧出现滑行波的现象。

3 曲线拟合与斜率估计

3.1 空气耦合波

针对电磁波的传播形式展开研究，如图2 所示的空气耦合波，由采集CMP 数据时收发间距的步进改变，可以构建电磁波传播时间的方程

其中d为收发天线的间距，c为在空气中的电磁波传播速度，t为电磁波由Q1传播至Q2的时间。

由式（1）可以看出，以收发间距d为自变量，以传播时间t为因变量，此方程在B-scan 图像中体现为一条直线，直线的斜率为空气中电磁波传播速度c的倒数。

3.2 地表滑行波

如图3 所示的地表滑行波，电磁波在空气-地面分界层下的介质中传播，可以依照式（1），得到相同形式的电磁波传播时间方程

其中，此时的vm为在介质中的电磁波传播速度，有

2.6 BRMS1、CD82、KISS1、TIMP1在侵袭性乳腺导管癌中mRNA表达降低 与非侵袭性乳腺导管癌相比，在侵袭性乳腺导管癌中抑制乳腺癌转移基因BRMS1、CD82、KISS1、TIMP1 mRNA 表达量均降低（P＜0.05），见图6。

由上可知，空气耦合波与地表滑行波在B-scan数据中具有相同的分布形式，通过两者曲线的拟合与斜率的求取，可以分别得到空气中的光速c与介质中的波速vm，反演出地表浅层的相对介电常数。

3.3 介质层反射波

针对接收回波信号中的反射波，如图4所示，建立电磁波传播时间方程

由式（4）可以推导出，反射波在B-scan 图像中，其收发间距d与传播时间t构成了双曲线的关系，不妨令一变量u，有u=。此时，u与t构成了线性关系，斜率为介质中传播速度的倒数。

根据B-scan 数据中的反射波双曲线信息，提取得到每一次测量下的反射波传播时间t=[t1，t2，…，tn]，对于每一个固定的层厚d1，通过求取相对应的变量u与t构成直线的斜率，可以得到与层厚一一对应的相对介电常数。于是，根据先验限制的地下厚度范围，假设一系列的介质层厚d1=[d11，d12，…，d1m]，可以拟合得到一系列的相对介电常数，绘制得到层厚d1与相对介电常数εr1一一对应的关系曲线。

当地下存在多层介质分层时，可以将多层介质看成单一均匀的介质层，采取单层CMP 模型下的反射波形式，可以推演得到u与t所构成直线，此时的直线斜率为均匀介质下的平均传播速度的倒数，通过曲线拟合得到的一系列相对介电常数也均为均匀介质下的等效相对介电常数，此时绘制层厚d1与等效相对介电常数εr1的关系曲线。

3.4 平均速度拾取

由于折射点位置的求取涉及一元四次方程的求解问题，并且折射点近似求取的方法存在一定偏差，于是如图5 所示，建立收发间距为0 m 的单层介质CMP模型，即属于垂直入射的传播模型，提取反射波传播时间T。由于本文方法不考虑薄层情况，即时域波形中各层反射波不存在混叠，此时根据零收发间距下采集的A-scan 波形，通过时域波形的幅度与相位变化，可以获取各层反射波的位置，提取得到反射波传播时间T，建立电磁波的传播时间方程

由式（5）与式（6）可知介质厚度d1与相对介电常数εr1一一对应，根据先验已知的相对介电常数范围，可以绘制出厚度d1与相对介电常数εr1的关系曲线。

针对地下存在多层介质分层的情况，可以将多层介质看成单一均匀的介质层，将多层介质模型简化为单层介质模型，采取以上收发间距为0 m的条件，获取介质层厚d1与等效相对介电常数εr1的关系曲线。

由图6 所示，基于曲线拟合与垂直入射模型下的参数关系曲线，其交点位置确定了介质的厚度信息与等效的相对介电常数，从而可以推导得到等效的平均速度，完成速度拾取。平均速度近似为均方根速度，通过引入Dix 公式，可以将均方根速度转换为每层介质的层速度，最后推导得到准确的相对介电常数与层厚参数。

3.5 Dix公式

将地下的多层介质分层等效为单一均匀介质层时，获得的速度近似为均方根速度，在基于水平层和层间速度恒定的假设下，Dix 公式可以将均方根叠加速度转换为平均层速度。因此，为了估计各层的相对介电常数，必须将均方根速度修正为层速度。第n层的平均层速度可用Dix公式计算

3.6 仿真结果

以有限差分时域（FDTD）仿真软件gprMax［21］作为仿真软件对本文所提出的层状介质参数反演算法进行实验验证。采取表1 的实验参数，建立图7所示的三层介质CMP 模型，由传播路径可以将回波信号分为空气直达波、地表滑行波、第一分层界面、第二分层界面、与第三分层界面反射波。

表1 仿真参数Tab.1 Simulation parameters

经过本文的算法流程可以获得如图8所示的仿真结果。图8（a）为基于三层介质的CMP 模型下回波B-scan 图像，通过波形特征获取到达时，绘制得到（b）中各波形曲线的位置，由上至下分别为空气直达波、地表反射波、第一分层界面、第二分层界面与第三分层界面的反射波。对波形位置进行曲线拟合与速度拾取，得到（c）、（d）、（e）所示的不同分层介质下参数关系曲线，交点位置确定了介质的厚度信息与等效的相对介电常数，推导得到等效的平均速度，近似为均方根速度。

由表2 所示的仿真结果可以看出，未引入Dix公式的算法在多层情况下无法有效反演介质的介电常数，这是由于曲线拟合得到的速度近似为均方根速度，不是介质的平均层速度。于是算法引入Dix 公式，将均方根速度转换为介质层速度。经过仿真实验验证，基于曲线拟合的参数反演算法可以将相对介电常数与层厚的反演误差均控制在5%以内，具有不错的反演准确率与反演效率。在实际测量中，由于大部分探地雷达采取发射和接收天线整体安装的形式，因此实测时需要搭配矢量网络分析仪与导轨，放置多个天线，设置测量迹线，进行CMP数据的采集。

表2 仿真结果Tab.2 Simulation results of the model a

4 结论

本文提出的基于曲线拟合的层状介质参数反演算法，打破了仅从空气耦合波与地表滑行波拟合出发的浅层参数反演的方法，而是从B-scan 数据中的反射波位置出发，提取得到相对介电常数与厚度信息的关系曲线，通过垂直入射模型，将多层介质模型等效为均匀单层模型，由Dix 公式将均方根速度转变为平均层速度，反演得到地下介质结构的参数信息。经过电磁仿真软件gprMax的仿真验证，证实本文所提方法具有很高的参数反演效率与反演准确率，相对介电常数与层厚的反演误差均控制在5%以内。