利用探地雷达频谱反演层状介质几何与电性参数

黄忠来，张建中

1厦门大学信息科学与技术学院，厦门 361005

2中国海洋大学海洋地球科学学院，青岛 266100

1 引 言

探地雷达（GPR）具有分辨率高、无损、快速等优点，是探测近地表介质的主要技术之一.利用GPR，不仅可以得到地层的埋深和厚度，也可以获取地下反射面的反射系数以及地层的介电常数和电导率等参数，为此，人们在正演［1－4］和反演方面［5－9］做了很多研究.目前层状介质参数的反演方法大体可分为时域反演和频域反演两类.在时域反演算法中，Chien等［10］利用地震勘探中的共中心点方法（CMP），求取地下两层介质的厚度和介电常数，张蓓［11］根据路面层状结构，在时域对探地雷达波进行反演，得到层的复介电常数和厚度，但是这两种方法都要求在时域能够分辨出不同反射层的回波；另一些算法则要求有先验信息作为基础，如Huang等［12］采用参数校验方法，通过在介电常数计算式中引入校验参数来求取水平路基的层厚和介电常数，但需要事先利用钻孔得到地下介质的有关参数.在频域反演方面，通过建立频域正演模型，不断改变模型参数，使模型计算的回波频谱与实际数据频谱一致，获取地层的厚度和电性参数.Qin等［13］通过求取电磁波在层状介质中的传输函数，反演地层的厚度和广义反射系数，但不能求得地层的电性参数；Minet以及 Lambot等［14－16］利用 Green函数建立目标函数，反演厚度和电性参数，但参数过多时，只有给出目标层的先验信息后才能取得较好的结果.谱反演算法［17］是最近提出的利用地震记录的频谱，反演地层深度、厚度和反射系数的新方法.与时域反演方法相比，谱反演算法对薄层（如油气储层）有较强的适应性.本文在谱反演算法［17］的基础上，结合电磁波在层状介质中的传播规律和GPR信号频谱特征，提出了一种利用地下层状介质GPR回波信号频谱，同时估算介质几何参数和电性参数的频域反演算法；通过分析模型参数对GPR频谱的影响规律，提出了分别单独估计一个或几个参数的分步反演方法，并以此结果作为进一步同时反演所有参数的初始值.最后，用理论模型的模拟数据和实测数据对方法进行了测试.

2 频谱反演方法原理

GPR 测 量 目 标 往 往 处 于 天 线 场 远 区［18－19］，且GPR多发射高频电磁脉冲.利用傅里叶变换，电磁脉冲可以分解为一系列不同频率的谐波，而这些谐波的传播都可以近似为平面波的传播形式.可见GPR的理论基础是平面波在介质中的传播规律［20］.因此，本文将雷达波作为平面波处理.对于GPR发射的电磁波来说，地下介质实际上是非理想的有耗介质，电磁波在传播时会逐渐衰减，影响传播的因素除了反射系数、透射系数外还有衰减常数.当介质为非理想介质时，电导率σ≠0，等效介电常数为复数，此时波数传播常量的定义为［21］

图1 探地雷达波在地下传播路径示意图Fig.1 Propagation path of GPR wave in layered media

图1所示的是包含4个反射系数的5层介质模型，虚线表示地面，其z轴坐标定为0，地下第i个反射层的深度用di表示，从上到下各层的介电常数分别为ε1、ε2、ε3和ε4，各层电导率分别为σ1、σ2、σ3和σ4，空气的介电常数为ε0，电导率为σ0.w（t）为雷达发射波，w0（t）、w1（t）、w2（t）、w3（t）和w4（t）分别是在地表和地下各反射面的一次反射回波.这里，我们假设介质为非磁性介质，并仅考虑雷达波入射角接近0°，即单天线形式（Monostatic mode）雷达或发射天线和接收天线之间的距离很小的情况.为清晰起见，图中将入射波角度放大，并将回波路径分开表示.第i层和i＋1层之间反射面的反射系数ri，i＋1和透射系数τi，i＋1分别为［22－23］：

且τi，i＋1＝1＋ri，i＋1.

式中，zi为波在第i层单向传播的距离，α0＝0.可以看到，广义反射系数是电磁波在介质中传播距离以及介质电性参数的函数.

定义了广义反射系数后，我们构造用于反演目标层深度、厚度以及电性参数的代价函数.以下考虑地下含有两个水平反射界面的三层介质情况，如图2a所示.若上、下广义反射系数分别为r1和r2，将地面作为时间零点，上界面到地面的时间距离为t0，下界面到地面的时间距离为t1，中间层的时间厚度

图2 三层介质反射系数模型Fig.2 Reflection coefficients of a three－layer model

为T，则反射系数序列可以表示为

如果将分析窗口的时间零点放在中间层的中点位置，那么反射系数序列的表达式变为

对上式进行傅里叶变换，得到该反射系数序列的频谱：

将它表示成三角函数形式：

令re＝（r1＋r2）/2，ro＝ （r1－r2）/2，则［17］：

其实部和虚部分别为

若雷达发射子波频谱为w（f），接收到的反射回波频谱为s（f），则地下介质的反射系数序列的频谱为

用雷达记录反射系数序列实部和虚部分别与褶积模型反射系数序列实部和虚部之差的绝对值之和，定义的代价函数为

其中，abs表示绝对值；fH和fL分别是所用频谱的频率上限和下限，将根据实际数据中所含噪声的分布情况确定，即尽量选取噪声较少的频带.该代价函数的未知量包括目标层上反射面的时间位置t0，目标层的时间厚度T，广义反射系数对的偶分量re和奇分量ro.由于广义反射系数是电性参数的函数，把代价函数中的广义反射系数奇、偶分量用电性参数替换，此时代价函数可表示成O（ε1，ε2，…，εn，σ1，σ2，…σn，t0，T）.通过求解代价函数的最小值问题，就可以得到地下各层的时间位置、厚度和电性参数，从而实现地下介质几何参数和电性参数的同时反演.

3 反演问题的解法

3.1 全局优化算法

对于求解像（17）式的多参数代价函数的极小值问题，我们采用在模拟退火算法上改进的随机爬山法［24］.算法步骤为

（1）产生初始未知参数向量X＝｛x1，x2，…，xi，…，xN｝.每个参数的值可以在规定的范围内随机生成，也可以根据先验信息给定.记初始参数向量对应的代价函数为O.

（2）在｛x1，x2，…，xi，…，xN｝中随机选取一个参数进行修改.例如选取的是第i个变量xi，则修改后xi变为，且＝xi＋sign·d·rand.sign是随机选取的正负符号，－1或1.rand为0～1之间的随机数，d为修改步长，用以控制修改参数的速度和精度.对于本文的几何和电性参数，需给定不同的d值.

3.2 分步反演方法

通过分析有关参数与GPR回波信号频谱之间的变化特征，提出了分别确定未知参数的分步反演方法.

考察广义反射系数序列的频谱：

其幅度谱为

对上式求导并令导数为0，即，

解之得

f表示幅度谱极大值点和极小值点的频率，这些极大值点和极小值点分别是周期性分布的，相邻极小值点频率之差，即幅度谱凹陷周期为

（22）式表明，幅度谱凹陷周期由时间厚度T决定，而与其它参数无关.例如当T分别取T1＝6.67ns和T2＝3.33ns，其它参数保持不变时，幅度谱凹陷周期分别为Δf1＝150MHz和Δf2＝300MHz，如图3所示.利用这一点，我们就可以在不考虑其它参数的情况下，仅利用反射系序列幅度谱凹陷周期来确定T.当层的厚度很小时，反射系数幅度谱凹陷周期将很大，若超出所设定的反演频带，就不能使用（22）式估算T.这时，可以找出幅度谱的第一个极大值点频率，峰值频率fp，再利用（21）式估算T，在式（21）中令n＝2，即得到T＝1/（2fp）.

从（21）和（22）式可以看出，随着层厚度的减小，反射系数幅度谱的第一个峰值频率和幅度谱凹陷周期将不断增大.即使当层厚度小于调谐厚度时，反射系数幅度谱也有这样明显的变化规律，从而可以利用频谱分辨或反演小于调谐厚度的薄层.

图3 幅度谱凹陷周期与层时间厚度T的关系Fig.3 Relationship between amplitude spectral notches period and layer′s temporal thickness T

将相位谱分为两个部分，第一部分不包含参数t，第二部分不包含反射系数的奇偶分量：

其中，

由于ro/re与频率f无关，所以相位谱∠r（f）随频率f变化的快慢由t0和T决定，且θ1随频率f的变化速度比θ2慢2πft0.当T按上述方法被确定后，就可以利用相位谱的变化快慢估算t0.为说明ro/re对于相位谱的影响，我们随机给定一组参数向量，并使其它参数保持不变，只改变ro/re的值，相应的相位谱如图4所示.Δf是相位发生跳变前持续的频率间隔，基本不随ro/re变化，这说明ro/re比值主要影响相位谱的幅值，而对相位谱变化快慢的影响很小.这样，在确定T参数后，就可以单独估算t0.

图4 相位谱随ro/re值的变化Fig.4 Variation of phase spectrum with the value of ro/re

考察相位谱中不含参数t0的项θ1.由式（3）可以得，当地下有三层介质时，

由于r1和r2只与介电常数有关，α2由第i＋2层介质的电导率以及介电常数决定，z2由T以及第2层的介电常数决定.即当T和t0被确定以后，影响相位谱的只有3层的介电常数和第2层的电导率.这样，仅利用相位谱就可以反演三个介电常数和中间层的电导率.

通过上面分析可知，不同参数对反射系数序列频谱的影响是不同的.T单独决定了幅度谱凹陷周期大小，t0和T共同影响相位谱的变化快慢，t0和各层的介电常数以及中间层电导率则完全确定了相位谱.因此，可以这样分别确定这些参数：首先根据幅度谱的凹陷周期确定厚度参数T，再根据相位谱的变化速度反演时间位置参数t0，然后根据相位谱反演3层的介电常数和中间层的电导率，最后反演上层的电导率.采用随机爬山法分别反演t0、三个介电常数以及中间层电导率.

这里用理论数据来说明分步求取参数的有效性.图5和图6中的实线是利用正确参数向量生成的反射系数序列的频谱；点线是利用估算的参数向量生成的反射系数序列的频谱.其中，图5（a，b）中的点线分别是由随机生成的一组参数向量计算出的反射系数序列频谱的幅度谱和相位谱，与理论频谱相差较大；图5c中点线是利用幅度谱凹陷周期估计T后计算的振幅谱.可以看到在其它参数保持为原随机生成的参数值的情况下，估算的T使幅度谱凹陷周期与理论一致，当然相位谱仍然相差很大，如图5d所示.图6（a，b）中的点线是在图5估算的T的基础上，利用相位谱估算t0，且其它参数不变时计算的幅度谱和相位谱.可以看出，估算的t0值使相位谱更接近理论值，且t0的改变并不影响幅度谱的凹陷周期.图6（c，d）中的点线是由前面确定的T和t0值以及继续反演三层的介电常数以及中间层的电导率后计算的幅度谱和相位谱，这时计算的频谱与理论频谱基本一致，其小的差别是因为上层的电导率没有估算而随机选取的缘故.可以看出，这样求出的介质参数值与理论参数值比较接近，计算的频谱也与理论频谱比较接近.这种分步求取参数的好处是可以极大地减少每次反演的未知量个数，降低收敛到局部最小值的几率，提高运算效率和反演结果的可靠性.特别是能为在没有先验信息情况下的多参数全局优化算法提供接近真解的初值.

4 实验例子

分别利用理论模型合成数据和实测资料对本文方法进行了测试和应用.理论模型采用楔形目标模型，以测试反演方法对不同厚度目标层的有效性.实测资料是收集到的公路路面GPR检测数据.

4.1 理论模型实验

使用GPRMax［25］建立了楔形目标层模型并合成了理论GPR数据.发射雷克子波中心频率为300MHz；上中下三层的介电常数ε1、ε2和ε3分别为9、16和9，电导率σ1、σ2和σ3分别为0.02、0.01和0.02；楔形层上界面的反射系数为－0.25，下界面反射系数为0.25；楔形层上界面距离地表0.38m，楔形层的最大厚度为0.125m，是雷达发射波在楔形中波长的1/2，1/4波长的调谐厚度位于第25道处.模型及楔形层上下反射面位置、楔形厚度以及各层电性参数的反演结果如图7和图8所示.可以看到，当目标厚度小于调谐厚度时，依然可以得到正确的反演结果.

现在以模型第35道为例，说明采用分步反演结果作为初始值的效果.图9给出了分别采用不同初始值时，代价函数随迭代次数的变化.实线为初始值采用分步反演结果，当迭代次数为1200次时就已收敛，对应的代价函数值是0.2834；点划线的初值是随机生成的，当迭代次数超过2500时才收敛，对应的代价函数是0.4269.显然，用分步反演结果作为初值时，反演收敛得更快且收敛于更小的代价函数值.表1列出了参数理论值、用作初值的分步反演结果和同时反演结果，图10是用该反演结果计算的幅度谱（a）及相位谱（b）与理论频谱之间的对比.

从上述反演结果可以看到，当楔形层的厚度大于1/8波长（即图7中第37道位置）时，本文方法可以给出较准确的上反射面时间位置和楔形层的时间厚度，电性参数也基本都在理论值附近；当楔形层厚度小于1/8波长后，反演结果误差较大.总的来说，深度和厚度参数的反演结果优于电性参数的结果.对目标层的时间厚度T的反演较准确，这是因为幅度谱凹陷周期大小主要与它相关，而与其它参数的联系不密切，分步反演给出了非常接近真值的初始值.同样，使用与反射面位置参数t0相关而与其它参数不相关的相位谱属性，使t0的反演结果也比较准确.介电常数既影响相位谱，又影响幅度谱.当需要反演的介电常数个数较多时，反演结果就容易陷入局部最优，这也是介电常数的反演结果没有T和t0结果准确的主要原因.相比之下，电导率反演结果不够准确，因为它决定雷达波的衰减程度，主要体现在雷达信号的幅值变化上，但幅值不仅与电导率有关，还与反射面上、下介电常数的差值以及层厚度有关系.

表1 采用分步反演数据为初始值的参数反演结果Table 1 Parameter results of inversion starting with initial values given by the staged inversion method

为了进一步说明本文分步反演方法的性能，将本文方法与参考文献［11］中的频域反演方法做了比较.文献［11］中的方法是针对地下三层介质，在地下第一层介质厚度、介电常数和电导率已知的情况下，反演第二层厚度和两个反射界面的广义反射系数.因此，设计了中间有一薄层的三层模型.第一层厚度为40cm，中间薄层厚度为3cm，约为电磁波在层中波长的1/7.从上到下三层的相对介电常数分别是9，25，9，电导率为0.02，0.01，0.02.发射天线中心频率为300MHz.分二种情况对这两种方法进行比较：方案一是按照文献［11］方法已知第一层的厚度和电性参数，随机给定其余参数的初值进行反演；方案二是随机给定包括第一层介质参数的所有参数的初值进行反演.两种方法反演过程的代价函数随迭代次数的变化曲线如图11所示，反演出的中间薄层的厚度见表2.可以看出，当第一层的厚度和电性参数已知时，两种反演方法都能很快收敛，且反演结果与理论值误差较小.当所有参数初始值随机选取时，两种方法达到收敛所需要的迭代次数增加，本文方法大约在迭代1200次时收敛，这时反演的薄层厚度与理论值的误差为0.26cm，文献［11］方法则大约在迭代2300次时才收敛，且代价函数仍然比前者大很多，反演的薄层厚度与理论值的误差为0.35cm.可见，在没有先验信息获得较好初值的情况下，本文的分步反演方法可以提高反演精度和效率.

表2 用两种方法反演薄层厚度（单位：cm）Table 2 Thin－layer thickness inverted using the two algorithms

4.2 实际资料应用

图11 两种反演方法代价函数随迭代次数的变化Fig.11 Cost functions versus iteration number for the two algorithms

第一个实测GPR数据来自广东某公路的混凝土路面检测资料.公路设计为水泥稳定碎石基层加水泥混凝土路面，采用C30混凝土，水泥层设计厚度为26cm，基层设计厚度为22cm.通常水泥层介电常数为5.5，电导率约为0.01～0.05；碎石介电常数约为5～15，电导率为3～20.实际钻芯位置在K0＋000－K1＋000处，钻探显示水泥的平均厚度为23cm.GPR天线发射频率为900MHz，每道采样时间长度为15ns，采集512个数据.图12是GPR实测数据剖面及反演结果，其中介电常数分别为4、6.5和15，电导率分别为0.01和0.025，把时间单位换算成长度单位，水泥层的平均厚度为24cm，基层平均厚度为21cm，反演结果与钻探结果基本一致，但水泥层厚度却略小于原设计厚度，也说明施工的水泥层厚度略微不够.

为了验证本文反演算法对于薄层的效果，我们收集了某沥青公路的GPR数据.雷达天线发射频率900MHz，每道采样时间长度为15ns，采集512个数据.雷达波在沥青中的波长大约为0.167m，沥青层上面层的厚度为0.04m，略小于1/4波长，因此可以被视为薄层.在时域剖面图上，沥青上面层的下底面回波与地表回波叠加在一起，很难对两者做出区分，而本文的频谱反演算法能给出沥青上层面的厚度.图13为第600道数据的上面层回波反射系数幅度谱，由于该层厚度小，难以确定幅度谱凹陷周期，但可以获取幅度谱的第一个峰值频率fp为940MHz，根据式（21）得到T 的初值为0.53ns，约0.041m.通过全局优化反演后得到该段沥青路面的上面层和下面层的深度和厚度如图14所示，其中反演的从上至下各层的平均相对介电常数的分别为0.83、5.3和3.8，平均电导率分别为6×10－6和0.005.

图12 某混凝土公路检测的GPR剖面（a）及反演结果（b）Fig.12 GPR profile from a highway detection and corresponding inversion result

图13 实际反射系数幅度谱和根据反演模型计算的幅度谱Fig.13 GPR amplitude spectra obtained from measured data and calculated using inversion model，respectively

5 结 论

通过对电磁波在层状介质中的传播规律和单天线模式或发射天线与接收天线间距很小时的GPR信号频谱特征的研究，基于地下层状介质GPR回波信号频谱及全局优化算法，提出了一种同时估算介质几何参数和电性参数的频域反演方法.理论模型和实际资料测试结果表明，该方法效果良好，具有较高的分辨率，对厚度为1/8发射脉冲波长左右的薄层介质仍有较好结果，而且对几何参数的反演精度优于对电性参数的反演.

地下层状介质的不同参数对于反射系数序列频谱属性的影响不同.如目标层时间厚度决定了幅度谱凹陷周期的大小，目标层时间厚度和上顶界面时间位置决定相位谱的变化快慢，上顶界面时间位置和各层的介电常数以及中间层电导率共同影响相位谱.据此，我们提出了利用不同的频谱属性分别估算不同未知参量的分步反演方法.以该分步反演结果作为全局优化反演的初值，可以大大提高反演的收敛速度和反演结果的可靠性.

图14 某沥青路面的GPR剖面（a）及反演结果（b）Fig.14 GPR profile of asphalt pavement and corresponding inversion result

（References）

［1］ 冯德山，张彬，戴前伟等.基于速度估计的改进型线性变换有限差分偏移在探地雷达中的应用.地球物理学报，2011，54（5）：1340－1347.Feng D S，Zhang B，Dai Q W，et al.The application of the improved linear transformation of finite difference migration based on the velocity estimation in the GPR date processing.Chinese J.Geophys.（in Chinese），2011，54（5）：1340－1347.

［2］ 冯晅，邹立龙，刘财等.全极化探地雷达正演模拟.地球物理学报，2011，54（2）：349－357.Feng X，Zou L L，Liu C，et al.Forward modeling for fullpolarimetric ground penetrating radar.Chinese J.Geophys.（in Chinese），2011，54（2）：349－357.

［3］ 田钢，林金鑫，王帮兵等.探地雷达地面以上物体反射干扰特征模拟和分析.地球物理学报，2011，54（10）：2639－2651.Tian G，Lin J X，Wang B B，et al.Simulation and analysis reflections interference from above surface objects of ground penetrating radar.Chinese J.Geophys.（in Chinese），2011，54（10）：2639－2651.

［4］ Huang Y Q，Zhang J Z，Liu Q H.Three－dimensional GPR ray tracing based on wavefront expansion with irregular cells.IEEE Transactions on Geoscience and Remote Sensing，2011，49（2）：679－687.

［5］ Huang Y Q，Liu Y H，Liu Q H，et al.Improved 3－D GPR detection by NUFFT combined with MPD method.Progress in Electromagnetics Research，2010，103：185－199.

［6］ Liu C R，Li J，Gan X，et al.New model for estimating the thickness and permittivity of subsurface layers from GPR data.IEE Proc.－Radar Sonar Navig，2002，149（6）：315－319.

［7］ Meles G A，Van der Kruk J，Greenhalgh S A，et al.A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole－to－surface GPR Data.IEEE Transactionson Geoscienceand Remote Sensing，2010，48（9）：3391－3407.

［8］ 王兆磊，周辉，李国发.用地质雷达数据资料反演二维地下介质的方法.地球物理学报，2007，50（3）：897－904.Wang Z L，Zhou H，Li G F.Inversion of ground－penetrating radar data for 2D electric parameters.Chinese J.Geophys.（in Chinese），2007，50（3）：897－904.

［9］ 丁亮，韩波，刘润泽等.基于探地雷达的混凝土无损检测反演成像方法.地球物理学报，2012，55（1）：317－326.Ding L，Han B，Liu R Z，et al.Inversion imaging method for concrete non－destructive testing based on GPR.Chinese J.Geophys.（in Chinese），2012，55（1）：317－326.

［10］ Kao C P，Li J，Wang Y Q，et al.Measurement of layer thickness and permittivity using a new multilayer model from GPR data.IEEE Transactions on Geoscience and Remote Sensing，2007，45（8）：2463－2470.

［11］ 张蓓.路面结构层材料介电特性及其厚度反演分析的系统识别方法：路面雷达关键技术研究［博士论文］.重庆：重庆大学，2003.Zhang B.System identification method for backcalculating the dielectric property and thickness of pavement structures：study on applied technology of ground penetrating radar（in Chinese）［Doctor′s thesis］.Chongqing：Chongqing University，2003.

［12］ Huang C L，Su Y.A new GPR calibration method for high accuracy thickness and permittivity measurement of multilayered pavement.∥Proceedings of the Tenth International Conference on Ground Penetrating Radar.Delft，The Netherlands：IEEE，2004：627－630.

［13］ Qin Y，Chen J，Fang G Y，et al.Research on thin－layer recognition technique based on the spectrum inversion method of ground penetrating radar.Journal of Electronics ＆Information Technology，2010，32（11）：2760－2763.

［14］ Minet J，Patriarca C，Slob E C，et al.Characterization of layered media using full－waveform inversion of proximal GPR data.∥2010URSL International Symposium on Electromagnetic Theory.Berlin：IEEE，2010：1004－1007.

［15］ Minet J，Lambot S，Delaide G，et al.A generalized frequency domain reflectometry modeling technique for soil electrical properties determination.Vadose Zone Journal，2010，9（4）：1063－1072.

［16］ Patriarca C， Lambot S， Mahmoudzadeh M R，et al.Reconstruction of sub－wavelength fractures and physical properties of masonry media using full－waveform inversion of proximal penetrating radar.Journal of Applied Geophysics，2011，74（1）：26－37.

［17］ Puryear C.I，Castagna J P.Layer－thickness determination and stratigraphic interpretation using spectral inversion：Theory and application.Geophysics，2008，73（2）：R37－R48.

［18］ 李大心.探地雷达方法与应用.北京：地质出版社，1994.Li D X.Methods and Applications of Ground Penetrating Radar（in Chinese）.Beijing：Geological Publishing House，1994.

［19］ Balanis C A.Antenna Theory：Analysis and Design.New Jersey：John Wiley，2005.

［20］ 曾昭发，刘四新，王者江等.探地雷达方法原理及应用.北京：科学出版社，2006.Zeng Z F，Liu S X，Wang Z J，et al.Principles and Applications of GPR Method（in Chinese）.Beijing：Science Press，2006.

［21］ Guru B S，Hızırolu H R.Electromagnetic Field Theory Fundamentals.Cambridge： Cambridge University Press，2004.

［22］ Schmelzbach C，Scherbaum F，Tronicke J，et al.Bayesian frequency－domain blind deconvolution of ground－penetrating radar data.Journal of Applied Geophysics，2011，75（4）：615－630.

［23］ Jol H M.Ground Penetrating Radar：Theory and Applications.Singapore：Elsevier，2009.

［24］ 张霖斌，纪晨，姚振兴.叠后地震道反演的随机爬山法.石油地球物理勘探，1997，32（1）：75－80.Zhang L B，Ji C，Yao Z X.Stochastic hill－climbing algorithm for poststack seismic trace inversion.Oil Geophysical Prospecting （in Chinese），1997，32（1）：75－80.

［25］ Giannopoulos A. Modelling ground penetrating radar by GprMax.Construction and Building Materials，2005，19（10）：755－762.