基于探地雷达频谱反演法的薄层识别技术研究

基于探地雷达频谱反演法的薄层识别技术研究

第32卷第11期 电 子 与 信 息 学 报 Vol.32No.11 2010年11月 Journal of Electronics & Information Technology Nov. 2010

基于探地雷达频谱反演法的薄层识别技术研究

秦 瑶

①

①②

陈 洁 方广有 阴和俊

①①③

(中国科学院电子学研究所 北京 100190)

②

(中国科学院研究生院 北京 100039)

③

(中国科学院 北京 100864)

摘 要：衡量探地雷达对薄层的识别能力通常有两个指标：最小可识别层厚和反射系数的识别精度。该文提出的频谱反演法，通过对电磁波在多层介质中频谱传输函数的推导，得出回波频谱的正演计算模型，然后采用阻尼最小二乘法进行二参数反演。通过实验，证实了这种方法能够对层厚小于八分之一波长的薄层进行层厚和反射系数的有效识别。

关键词：探地雷达；薄层；层厚；反射系数；频谱反演

中图分类号：TN959.71 文献标识码： A 文章编号：1009-5896(2010)11-2760-04 DOI: 10.3724/SP.J.1146.2009.01511

Research on Thin-layer Recognition Technique Based on the

Spectrum Inversion Method of Ground Penetrating Radar

Qin Yao

①②

①②

Chen Jie Fang Guang-you Yin He-jun

①①

③

(Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China) (Graduate University, Chinese Academy of Sciences, Beijing 100039, China)

③

(Chinese Academy of Sciences, Beijing 100864, China)

Abstract: Two parameters are normally employed to evaluate recognition ability of the Ground Penetrating Radar (GPR) for the thin-layer: the least recognition thickness and the reflectance recognition precision. In this paper, the spectrum inversion method is used to improve the recognition ability of the GPR. The spectrum transmission function of the multilayer media is firstly deduced, and the forward model of the reflecting wave spectrum is obtained. Then the two parameters are inversed by using the damp least-squares algorithm. The experiment demonstrates that the method can identify the layer effectively with a thickness less then one-eighth wavelength. Key words: Ground Penetrating Radar (GPR);Thin-layer; Thickness; Reflectance coefficient; Spectrum inversion

1 引言

层状介质是探地雷达(GPR)探测中常见的介质结构，如冰川、公路路基等[1 3]。薄层是指层厚小于四分之一波长的层。衡量GPR对薄层的识别能力通常有两个指标：最小可识别层厚和反射系数的识别精度。GPR对薄层的传统识别方法来源于地震勘探中广泛应用的振幅法[4,5]，这种方法对于反射系数已知的介质层的探测十分有效。然而，GPR探测中更普遍的情况是层厚及反射系数均未知，这时需要针对层状介质建立严谨的雷达波正演模型，进而采用多参数反演法进行薄层的识别。文献[3]针对路面结构层介质模型，利用雷达时域回波实现了多参数2009-11-27收到，2010-07-07改回

国家科技支撑计划(2008BADA8B02)资助课题 通信作者：秦瑶 eqinyao@http://www.77cn.com.cn

的反演，但没有涉及最小层厚的分辨率问题。

本文提出的频谱反演法，首先通过对电磁波在多层介质中传输函数的推导，得出严谨的正演模型，然后利用回波频谱进行反演，不仅能够准确得到薄层反射系数，而且提高了层厚的分辨率，从而增强了GPR对薄层的识别能力。

2 时域回波分析

2.1 层厚的影响

对于层状介质模型，如果目标层上下界面跟上下地层有较明显的电性差异，则在该目标层的顶面、底面上都能形成反射波。GPR对薄层厚度的识别能力就是区分上下界面反射波的能力，也就是垂直分辨率[6]。目前，对于GPR垂直分辨极限的定义，通

常参照Rayleigh标准，即薄层的分辨率极限为λ/4。

设雷达发射信号为高斯脉冲，经收发天线，接收机采集到的信号为二阶高斯微分信号。图1为反

基于探地雷达频谱反演法的薄层识别技术研究

第11期 秦 瑶等：基于探地雷达频谱反演法的薄层识别技术研究 2761

射系数R=0.33时，3种不同层厚对应的时域回波。可以看出，当薄层厚度l较大时，介质上、下层反射的信号在回波中能够明显区分开，可以直接通过回

随着薄层厚度波信号波峰与波谷的时间差值计算[5]。

的减少，薄层上下界面的回波在时间轴上逐渐重叠，但波形的幅度逐渐增大，这时可以利用回波幅度求取目标层厚，这就是振幅法的探测原理。 2.2 反射系数的影响

反射系数是表示层界面上反射强度的物理量，电磁波在不同介质的界面上的反射系数取决于两种介质的相对介电常数。两种介质的介电常数差别越大，反射系数的模值R也就越大，Rmax=1(全反射)；相反，如果两介质的介电常数非常接近或一样，反射系数模值R就很小，Rmin=0(无反射)。

反射系数对接收回波的影响，可由图2的仿真波形直观表述。可以看出，回波幅度随反射系数R的增大而增大。

由图2结合反射系数与介电常数的关系，我们能够得出这样的结论：界面两边介质的介电常数的差别越大，回波幅度越大；相反，如果两介质的介电常数非常接近或一样，回波幅度就很小，很难确定薄层的存在。因此，运用GPR识别薄层要求薄层的介电特性必须与上下地层有较明显的差别。

其中kiz=kicosθi，ki+1,z=ki+1cosθi+1；ki，ki+1分 别为电磁波在介质i和介质i+1中的复波数，电磁波垂直入射时，θi=θi+1=0；εi，εi+1分别为介质i和介质i+1的复介电常数。

区域1，区域2，区域3中的波可分别表示为

e2ik1zd1+ik1zz (3) e1y=A1 e ik1zz+R12

2ik2zd2+ik2zz

e2y=A2 e ik2zz+R23e (4)

e3y=A3e ik3zz (5)

为广义反射系数，它表示第1界面z= d其中R121

处所有上行波幅值和下行波幅值之比。由于3区延伸至无穷远，故式(5)中只存在下行波。

可由界面边界条件获未知数A1，A2，A3和R12

得。区域2中的下行波为区域1中下行波的透射波和区域2中上行波的反射波的迭加；在上交界面z= d1处，该边界条件可写为

ik1zd1

A2eik2zd1=AeT12+R21A2R23e2ik2zd2 ik2zd1 (6) 1

同时，区域1中的上行波是由区域1下行波的反射和区域2上行波的透射形成的。于是，在界面z= d1处，有如下约束条件：

eik1zd1=RAeik1zd1+TARe2ik2zd2 ik2zd1 (7) A1R1212121223由式(6)，A2可借助于A1表达，即

i(k1z k2z)d1

T12Ae1

(8) A2=

1 R21R23e2ik2zd2 d1将式(8)代入式(7)，得到

T12R23T21e2ik2z(d2 d1) (9) R12=R12+

1 R21R23e2ik2zd2 d1式(9)即为3层介质时的广义反射系数，它包含了第1层和第2层的反射影响。

由此递推, 并利用关系式Ti,i+1=1+Ri.i+1，可

[8] 以得到N层介质结构的广义反射系数Ri,i+1为

2iki+1,z(di+1 di) R+Rei,i+1i+1,i+2 R (10) i,i+1=2iki+1,zdi+1 di 1+Ri.i+1Ri+1,i+2e对探地雷达测量环境而言，收发天线均置于空

气中，因此最简单的测试模型应为4层，自上而下依次为：空气层、沙(土)层、异常层、沙(土)层。因

3 正演模型的建立

3.1 层状介质的传输函数

实现薄层参数的反演，首先需要推导层状介质的系统传输函数。设计简单的3层反射模型如图3所示[7]，源点和场点均在第1层中。以下推导过程中不考虑收发天线间的直达波，忽略层之间的多次散射回波。因此，图3只给出了电磁波在层状介质之间的主波路径示意图。

以TM波为例，设电磁波的入射方向与z方向夹角为θ1，波在界面 d1， d2上均会发生反射和折射，电场反射系数R和透射系数T分别为

εiki+1,z εi+1kiz

Ri,i+1=, i=1,2 (1)

εiki+1,z+εi+1kiz

2εiki+1,z

Ti,i+1=, i=1,2 (2)

εiki+1,z+εi+1kiz

图1 不同层厚l时的时域回波(R=0.33) 图2 不同反射系数时的时域回波(l=λ/2) 图3 电磁波在3层介质中的主波路径示意图

基于探地雷达频谱反演法的薄层识别技术研究

2762 电 子 与 信 息 学 报 第32卷

此需要在3层模型系统函数式(9)的基础上，将R23用

代替。令区域4与区域3的交界面广义反射系数R23

的表达式为 为z= d3，由式(10)得到R23

e2ik3z(d3 d2)R23+R34 R23= (11) 2ik3zd3 d21+RRe

23

34

频谱A。因此，4层介质模型回波频谱B可表示为

B=A F(l,R) (15) 式(15)即为雷达回波频谱的正演计算公式。

4 二参数反演

首先给出模型中待反演参数的初始假设值为

假设第4层延伸至无穷远，故该区内只存在下行波，满足：

R 34=R34

(12) 将式(11)，式(12)代入式(10)，再利用关系式

Ri,i+1= Ri+1,i，整理后得

R (1 R2iD12)R23e1+(1 R2i(D12)R34e1+D2)12=R12+1+RRiDReiD1D2 2334e2 12R23eiD1 R12R34(13)

式中D1=2k2z(d2 d1)；D2=2k3z(d3 d2)=2k3z

l，l=d3 d2，即为异常层的厚度。

由于探地雷达大多采用地面耦合的测量方式，

因此，收发天线间的直达波与第1层地表的反射波

R12在时域上无法区分，在滤波处理中会随直达波一

起滤除，因此，对于地面耦合探地雷达的4层介质

谱域传输函数应为

R (1 R212)R23eiD1+(1 R2

i(D12)R34e1+D2)12=1+RiDiD23R34e2 RiD12R23e1 R12R34e1+D2(14) 为简化书写、明确物理意义，令式(14)中的R 12

=F(l,R)。F(l,R)即为异常层层厚为l，异常层反射系数为R时4层介质的谱域传输函数。同理，

可按上述推导方法得到N层介质的谱域传输函数。

需要特别说明的是，式(14)中的参数除了目标

层厚度l和与薄层反射系数相关的目标层介电常数

ε3外，还有Δx，Δz，f，ε1，d1，ε2和d2。其中，

当测量系统一定时，收发天线相对位置Δx，Δz，

发射频率f，

空气介电常数ε1，地面的位置坐标d1均为已知量；对于一般的探测介质，ε2也是可以采用点源反射，TDR(Time Domain Reflectometry)等方法预先测量得到的；d2为异常体上表面的相对位置，可通过对时域波形的分析来确定。 3.2 回波频谱 回波频谱的正演计算还需要用到雷达的入射

波。由于雷达的型号、天线、测试环境都会对雷达

入射波产生很大影响，因此入射波不能简单地用一

个函数来定义。在每次试验前，需要采用金属板反

射法提取雷达入射波[3]。

具体实现方法是：将一块表面平整的金属板铺在介质表面，记录接收到的金属

板全反射波形；再将天线对向天空，记录无目标时

的回波波形；最后将两次记录的波形相减，即可得到雷达入射波。将其进行傅里叶变换，得到入射波

PT0=[p10,p20,",pM0] (16)

其中M为待反演参数个数。

设Y为测量回波的频谱值，用泰勒级数将Y展

开，略去二阶以上项，得

M yn=bnP bn 0+∑ i=1 pi Δpi, pi=pi0

n=1,2,",N; i=1,2,",M (17)

其中N为测量点数。将式(17)写成矩阵形式：

ε=Y B=JΔP (18)

式中Y=[y1,y2,",yTN]为测量数据矩阵；B=[b1,b2, ",bN]T

为假设模型的计算值；ε=[y1 b1,y2 b2, ",yTN bN]为测量值与计算值之差列阵；ΔP= [p1,p2,",pTM]为待反演参数增量阵列。J为Jacobi矩阵，本系统针对层厚和反射系数反演，Jacobi矩阵是二参数表达式如下： b1 b1 p 1

p2

J= ## b b (19) nn p1 p2

式(18)是关于待反演参数修正步长ΔP的N×M线性方程组。采用阻尼最小二乘法计算参数增量ΔP(阻尼最小二乘法的实现过程可参见文献[9,10])，由此得到新的参数列阵Pi=P0+ΔP，再将这组新的参数代入正演模型，可求出新的预测数据y

，新的差值列阵ε ，以及新的Jacobi矩阵；再用这些结果得出新的参数增量ΔP

i，由此构成迭代过程直到满足收敛条件为止。

反演迭代误差e的定义为：测量频谱值与反演频谱值之间相对误差绝对值的和，即

e(i)=∑(Y A F(li,Ri))i=1,2,",I (20)

其中I为迭代允许次数。当反演误差e(i)小于给定值，此时的li,Ri即为反演得到的薄层厚度和反射系数。 由于反演结果具有不唯一性，因此在反演过程中需要加上l>0的约束条件。在具体应用中，还可以利用待测介质的时域回波相位等先验信息来排除赝解。

基于探地雷达频谱反演法的薄层识别技术研究

第11期 秦 瑶等：基于探地雷达频谱反演法的薄层识别技术研究 2763

5 实验验证

实验采用拉脱维亚Zond 12c雷达进行探测。在沙坑中0.3 m深处埋设厚度为1.2 cm、ε木板=2.8的木板。用TDR测得ε沙=4，计算知反射系数R=0.0889。雷达天线中心频率为1.5 GHz，电磁波在沙中的波长λ=0.1 m。因此，目标层厚可以用波长表示为l=0.12λ。图4为进行滤波后的B-scan探测图，由图可以较清楚判断目标层的位置。同时，由图4提取含层目标的第150道回波数据进行反演计算。首先判断出目标上表面距地表距离约为0.3 m；再给出一组反演初始值，令两个待反演参数的

步长均为0.001，按式(19)构造Jacobi矩阵；之后代入式(18)得到反演方程；接着用阻尼最小二乘法(阻尼因子α=0.01)，迭代10次后输出反演结果，并计算层厚和反射系数的反演误差。

图5为金属板反射法得到的发射波形频谱。图6为反演初始值l0=0.2,R0=0.07时的反演频谱与测量频谱的对比图，由于实验噪声等影响，二者出现一定偏差。表1为不同反演初始值迭代10次的反演结果。可以看出，薄层厚度和反射系数的反演值与真实值的误差均小于6%，可以满足工程探测需求。

图4 B-scan测量图 图5 发射波形频谱 图6 计算频谱与测量频谱对比图

表1 不同初始值迭代10次的反演结果 初始值 层厚

反演值

误差(%)

University, 2008. [3] [4] [5] [6] [7]

张蓓. 路面结构层材料介电特性及其厚度反演分析的系统识别方法[D]. [博士论文], 重庆: 重庆大学, 2003.

序号

(单位:λ)

R

层厚

(单位: λ)

R 层厚 R

De Voogd N and Den Rooijen H. Thin-layer response and spectral bandwidth [J]. Geophysics, 1983, 48(1): 12-18. 陈先红. 基于叠前资料道集上的远近偏移距属性差异的薄层识别[J]. 油气地球物理, 2008, 6(4): 51-55.

袁明德. 浅析探地雷达的分辨率[J]. 物探与化探, 2003, 27(1):

1 0.2 0.07 0.119 0.0884 0.8 0.6 2 0.3 0.33 0.118 0.0898 1.7 1 3 0.5 -0.05 0.125 0.0876 4.2 1.5 4 0.6 -0.1

0.127 0.0870 5.8 2.1

28-32.

Lin Z W, Zhang X J, and Fang G Y. The second-order SPM solution for scattering from multi-layer dielectric media with slightly rough surface[C]. PIERS Proceedings, Moscow, Russia, Aug. 18-21, 2009: 1387-1390. [8] [9]

Weng Cho-Chew. Waves and Fields in Inhomogeneous Media [M]. NewYork: Van Nostrand Reinhold, 1990: 45-55. 张建华，刘振华，仵杰. 电法测井原理与应用[M]. 西安: 西北大学出版社, 2002: 189-193.

6 结论

本文提出采用频谱反演法提高探地雷达对薄层的识别能力，通过对电磁波在分层媒质中频谱传输函数的推导，实现雷达回波频谱的正演计算，进而构造反演模型，采用阻尼最小二乘法实现了二参数反演。通过对实验数据的处理，证实了频谱反演法能够对层厚小于传统Rayleigh标准的薄层进行有效识别，使探地雷达对薄层的识别能力得到了提高。

参 考 文 献

[1]

黄饶, 陈小宏, 李景叶. 层间多次波模拟与影响分析[J]. 地球物理学进展, 2009, 24(3): 981-987.

[10] 刘国林. 非线性最小二乘与测量平差[M]. 北京: 测绘出版社,

2002: 12-13.

秦 瑶： 女，1981年生，博士生，研究方向为雷达信号处理及成

像技术研究.

陈 洁： 男，1978年生，博士，助研，研究方向为雷达信号处理

及成像技术研究.

方广有： 男，1963年生，研究员，博士生导师，中科院“百人计

划”入选者，研究方向为超宽带电磁学及其工程应用、损耗介质中隐蔽目标的电磁探测、微波成像新技术和新方法和电磁场的数值计算方法等.

[2] Seyed Alireza Tabatabaeenejad. Forward and inverse models of electromagnetic scattering from layered media with rough interfaces[D]. [Ph.D. dissertation], Ann Arbor: Michigan

本文来源：https://www.bwwdw.com/article/8w3i.html