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)

其中，为衰减常数，单位为Np/m，表示电磁波在传播1 m后单位振幅衰减至e^-1；为相位常数，单位为rad/m，表示波传播1 m后的相移；ε为介电常数，σ为电导率，ω为角频率.对于理想介质，σ=0，则α=0，β=k=ω

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

(2)

(3)

且τ_i，i+1=1+r_i，i+1.

定义广义反射系数_i，i+1为各反射回波与雷达发射波的电场强度之比，则地下第一个反射界面的广义反射系数₁₂为

(4)

第二个反射界面的广义反射系数₂₃为

(5)

以此类推，第i个广义反射系数_i，i+1为

(6)

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

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

(7)

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

(8)

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

(9)

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

(10)

令r_e=（r₁ +r₂）/2，r_o=（r₁ -r₂）/2，则^[17]：

(11)

其中，Re[（f）]=2r_ecos（πfT），Im[（f）]=2r_osin（πfT）分别是实部和虚部.从定义式可以看出r_e和r_o分别是对（t）进行奇偶分解后所得的偶分量和奇分量^[17]，如图 2（b，c）所示.实际上，分析窗口的时间零点很难正好位于两反射系数的中点，因此，在式（11）的基础上导出窗口零点位置不在反射系数对中点的代价函数.当时间零点位置在地面时，反射系数序列r（t）为

(12)

容易看出，，其频谱为

(13)

其实部和虚部分别为

(14)

(15)

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

(16)

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

(17)

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

3 反演问题的解法 3.1 全局优化算法

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

（1）产生初始未知参数向量X={x₁，x₂，…，x_i，…，x_N}.每个参数的值可以在规定的范围内随机生成，也可以根据先验信息给定.记初始参数向量对应的代价函数为O.

（2）在{x₁，x₂，…，x_i，…，x_N}中随机选取一个参数进行修改.例如选取的是第i个变量xi，则修改后x_i变为，且=x_i+sign·d·rand.sign是随机选取的正负符号，-1或1.rand为0~1之间的随机数，d为修改步长，用以控制修改参数的速度和精度.对于本文的几何和电性参数，需给定不同的d值.

（3）用修改后的参数代替x_i，产生新的参数向量={x₁，x₂，…，，…，x_N }，它对应的代价函数为.

（4）若 < O，则接受本次对参数进行的修改；否则不接受修改，继续随机选择下一个变量.重复（2）到（4）的操作，直到最后计算出的代价函数值满足收敛条件，或者迭代次数超过设定值为止.每次修改时的步长d不是固定的，开始时为了加快收敛速度，步长可以取得较大，在算法的结果不再明显变化后，逐渐减小d，从而提高每次修改的精度.

3.2 分步反演方法

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

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

(18)

其幅度谱为

(19)

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

(20)

解之得

(21)

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

(22)

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

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

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

(23)

其中，

(24)

(25)

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

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

(26)

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

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

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

4 实验例子

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

4.1 理论模型实验

使用GPRMax^[25]建立了楔形目标层模型并合成了理论GPR数据.发射雷克子波中心频率为300 MHz；上中下三层的介电常数ε₁、ε₂和ε₃分别为9、16和9，电导率σ₁、σ₂和σ₃分别为0.02、0.01和0.02；楔形层上界面的反射系数为-0.25，下界面反射系数为0.25；楔形层上界面距离地表0.38m，楔形层的最大厚度为0.125 m，是雷达发射波在楔形中波长的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的反演较准确，这是因为幅度谱凹陷周期大小主要与它相关，而与其它参数的联系不密切，分步反演给出了非常接近真值的初始值.同样，使用与反射面位置参数t₀相关而与其它参数不相关的相位谱属性，使t₀的反演结果也比较准确.介电常数既影响相位谱，又影响幅度谱.当需要反演的介电常数个数较多时，反演结果就容易陷入局部最优，这也是介电常数的反演结果没有T和t₀结果准确的主要原因.相比之下，电导率反演结果不够准确，因为它决定雷达波的衰减程度，主要体现在雷达信号的幅值变化上，但幅值不仅与电导率有关，还与反射面上、下介电常数的差值以及层厚度有关系.

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

4.2 实际资料应用

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

5 结论

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

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

致谢

感谢重庆交通科研设计院为本文提供实际公路检测资料.