探地雷达薄层信号的谱反演算法＊

黄忠来 张建中 黄吉林

(1)厦门大学信息科学与技术学院,厦门 361005 2)中国海洋大学海洋地球科学学院,青岛 266100)

探地雷达薄层信号的谱反演算法＊

黄忠来1)张建中2)黄吉林1)

(1)厦门大学信息科学与技术学院,厦门 361005 2)中国海洋大学海洋地球科学学院,青岛 266100)

为了从多层介质的探地雷达回波数据中提取薄层的位置和厚度信息,推导和验证了在频率域内实现的谱反演算法。该算法给出了各反射面的广义反射系数,并建立了用于联系反射系数序列频谱和薄层参数的代价函数并进行反演。仿真结果表明,对于厚度小于调谐厚度的薄层,该算法也可以进行准确反演,从而提高了对于薄层的分辨能力。

探地雷达;薄层;谱反演算法;广义反射系数;代价函数

1 引言

探地雷达 (Ground Penetrating Radar,GPR)通过向地下发射电磁脉冲,并接收地下回波信号来探测地下目标,它具有分辨率高、无损、快速等优点。本文研究的对象是含有薄层的层状结构地下介质,目的是通过反演算法得到地下各反射面的位置以及薄层的厚度,这里的薄层是指厚度小于电磁波在其中波长 1/4的层。对于层状介质参数的反演,文献[1]针对地震勘探共中心点方法 (CMP)中的模型进行改进,并在时域反演,以求取地下两层介质的厚度和介电常数,但是反演时需要在时域得到反射面回波的到达时间。文献[2]研究路面层状结构,利用雷达波传播规律进行时域反演,得到层的复介电常数和厚度,但是同样没有考虑目标层是薄层的情况。本文在文献[3]提出的地震勘探谱反演算法的基础上,结合电磁波在层状介质中传播的特点,提出用广义反射系数代替原算法中的反射系数,推导并验证了对地下薄层 GPR信号进行反演的频域算法,提高了对地下薄层的分辨力。

2 谱反演算法过程

设 re=(r1+r2)/2,ro=(r1-r2)/2,则：

其中 Re[g(f)]=2recos(πfT),Im[g(f)]=2rosin(πfT)分别是实部和虚部。从定义式可以看出 re是函数 f(r1,r2)=r1δ(t-t1)+r2δ(t-t2)做奇偶分解后的偶部分,ro则是奇部分(图 2)。

令 G(f)表示 g(f)的幅度谱,则：

图1 反射系数对Fig.1 Reflection coefficients pairs

图2 反射系数对奇偶分解Fig.2 Even-odd decomposition of coefficients pairs

对 G(f)取导数得：

将 (3)、(4)两式相乘得：

找到使得代价函数值为零的参数和后,就可以根据

求出上、下反射系数的奇偶分量,进而得到反射系数值。

由于事先并不知道反射系数谱 g(f),所以在计算 G(f)时要利用接收到的雷达回波和雷达发射子波的频谱。此外,算法建立代价函数的前提是分析窗口时间零点位置位于两反射系数的中点。下面,将在此算法的基础上推导适用于多反射系数的代价函数,本文中所推导出的函数表达式和文献[3]中的有些不同。

以包含 4个反射系数的模型为例,图 3中反射系数 r1～r4分别为 0.1、0.3、-0.2、0.4,t1～t4分别是 4个反射系数的距离地面的时间位置,T1、T2是前后两对反射系数的时间厚度,则反射系数序列g(t)及其频谱 g(f)可以表示为 g(t)=r1(t)+r2(t) =[r1δ(t-t1)+r2δ(t-t1-T1)]+[r3δ(t-t3)+r4δ (t-t3-T2)],以及 g(f)=r1(f)+r2(f)。

图3 反射系数对及其奇偶分解Fig.3 Coefficients pairs and their even-odd decompositions

由前面两反射系数的推导可知：

将 4个反射系数 (两个反射系数对)推广到 n个反射系数对,则代价函数可表示为：

此代价函数含有以下未知量：反射系数对中,上层反射系数的时间位置 t,两反射系数之间时间厚度T,反射系数对的偶分量 re以及反射系数对的奇分量 ro。以上就是谱分解算法的大致过程。通过最小化代价函数的值,可以得到地下各层的时间位置、厚度以及反射系数。但是,算法中讨论的是简单的反射系数褶积模型,未讨论波在地下传播时的衰减。对于探地雷达波来说,反演出的反射系数并不等于反射面真正的反射系数,因为地下介质实际上是非理想的有耗介质,导致电磁波在传播时会逐渐衰减,影响其传播的除了反射系数、透射系数外还有衰减常数。因此原代价函数里的反射系数奇偶分量需要重新定义。

3 广义反射系数

同理,第二个反射界面的广义反射系数为：

以此类推,第 i个广义反射系数为：

式中 zi为波在第 i层中单向传播的距离,c为光速,并认为α0=0。

图4 电磁波在地下传播路径示意图Fig.4 Electromagnetic wave propagation path under ground

4 反演算法

对于包含多参数的代价函数的反演,本文采用在模拟退火算法上改进得到的随机爬山法[4]。其算法步骤为：1)产生初始参数向量 X={x1,x2,…, xi,…,xN},每个参数的值可以在规定的范围内随机生成,也可以根据实际情况下得到的数据进行设定。这时,向量对应的代价函数为Obj。

2)在{x1,x2,…,xi,…,xN}中随机选取一个参数进行修改,例如选取的是第 i个变量 xi,则修改后xi变为且=xi+sign＊d＊ rang,sign是在 -1和 1两者间随机选取的符号,d为修改的步长,用以控制反演的速度和精度,rand为 0～1之间的随机数。

重复 2)到 4)的操作,直到最后计算出的代价函数值满足收敛条件,或者反演的次数超过设定值为止。每次修改时的步长 d不是固定的,开始时为了加快收敛速度,步长可以取得较大,以后逐渐减小,从而提高每次修改的精度,尽量使算法收敛到全局最优值。

5 数据仿真

5.1 楔形模型仿真

为讨论目标层厚度对算法结果的影响,利用GPRMax软件[5]进行正演仿真,产生作为实际参考值的数据。建立上、下层面反射系数同号和异号两种情况下的楔形模型A和B(图 5和图 6分别是A和B的时域回波图)。模型 A中,令楔形上界面反射系数为 -0.25,下界面为 0.25。模型 B中的楔形,令上界面反射系数 -0.14,下界面为 -0.11。楔形厚度最厚处为发射波在楔形中波长的 1/2,雷达发射波中心频率设定为 300 MHz(图 5～8)。

图7和图 8分别是楔形模型 A和 B的反射系数反演结果。从图7和图8可以看到,在两种反射系数组合下的模型中,1/4波长调谐厚度均位于记录道第 25道,当楔形厚度逐渐减小到 1/8波长时,算法仍然可以正确反演出反射面的位置和厚度,这说明了算法对于薄层的分辨效果非常好。

图6 模型B回波Fig.6 Received waves ofmodelB

图7 模型A反射系数反演结果Fig.7 Inversion results of coefficients in modelA

图8 模型B反射系数反演结果Fig.8 Inversion results of coefficients in modelB

5.2 多反射系数层状模型

设计模型在地下有 4个反射界面,反射系数分别是 r1=-0.250 0、r2=0.111 1、r3=0.200 0、r4= 0.500 0。其相应的广义反射系数分别是 gr1= -0.113 1、gr2=0.043 7、gr3= -0.030 6、gr4= 0.069 9。地下包含了 2个层,上方层的厚度设计成电磁波在层中波长的 1/10,下方层的厚度设计成波长的 1/2,即上方是个相对薄层,而下方是相对厚层。

表1 实际值与反演结果比较Tab.1 Comparison between i nversion results and true values

图9 回波、发射波及反射系数序列频谱Fig.9 Spectrum of reflected wave,incident wave and reflection coefficients

图10 反演结果与实际值比较Fig.10 Comparison between inversion results and true values

图9是雷达回波、发射波的频谱以及由前两者计算出的模型反射系数序列频谱,图 10是由反演出的参数计算得到的频谱和真实反射系数序列频谱之间的对比,上方是两个频谱的实部,下方是虚部,其中红色星划线是真实值,蓝色圈划线是反演值。从表 1的多个反射面模型反演结果可以看到,薄层的位置以及厚度的反演结果准确,表明当地下为多层结构时,算法仍然正确有效。其中薄层的广义反射系数反演值与真实值存在一定误差。

6 结语

通过研究电磁波在地下层状介质中的传播规律,重新定义了反射系数表达式,推导并验证了适用于 GPR信号的谱反演算法。正是因为有了新的广义反射系数定义式,才可以根据实际模型参数来验证算法反演结果的正确性。利用谱反演算法可以得到厚度小于 1/8波长的薄层的反射面位置,层厚以及每个反射面广义反射系数的大小。仿真结果表明该算法对于提高薄层的分辨有很好的效果,但是在优化算法的结果上存在多解性问题,有进一步改进的余地。

1 Kao Chienping,et al.Measurement of layer thickness and permittivity using a new multilayer model from GPR data [J].IEEE Transactions on geoscience and remote sensing, 2007,45(8)：2 463-2 470.

2 张蓓.路面结构层材料介电特性及其厚度反演分析的系统识别方法——路面雷达关键技术研究 [D].重庆大学,2003. (Zhang Bei.System identification method for back-calculating the dielectric property and thickness of pavement structures_study on applied technology of ground penetrating radar[D].ChongqingUniversity,2003)

3 Charles I Puryear and John P Castagna.Layer-thickness deter mination and stratigraphic interpretation using spectral inversion：Theory and application[J].Geophysics,2008, 73(2)：37-48.

4 张霖斌,纪晨,姚振兴.叠后地震道反演的随机爬山法[J].石油地球物理勘探,1997,32(1)：75-80.(Zhang Linbin,JiChen and Yao Zhenxing.Stochastic hill-climbing algorithm of post-stack seismic inversion[J].Oil Geophysical Prospecting,1997,32(1)：75-80)

5 Giannopoulos A.Modelling ground penetrating radar by Gpr Max[J].Construction and Building Materials,2005, 19：755-762.

6 秦德文.基于谱反演的薄层预测与反演方法研究[D].中国石油大学,2009.(Qin Dewen.Thin bed prediction using spectral inversion and inversion method analysis[D].China University of Petroleum,2009)

7 Liu C R,et al.New model for estimating the thickness and per mittivity of subsurface layers from GPR data[J]. IEE Proe.-Rodur SonorNaviz,2002,149(6)：315-319.

A SPECTRAL INVERSI ON ALGORITHM FOR GPR SIGNALS OF THIN LAYERS

Huang Zhonglai1),Zhang Jianzhong2)and Huang Jilin1)

(1)School of Infor m ation Science and Technology,X iam en University,X iam en 361005 2)College of M arine Geo-science,Ocean University of China,Q ingdao 266100)

In order to extract infor mation including thin layer’s position and thickness,from reflected wave data of Ground Penetrating Radar(GPR)in multi-layered situation,a spectral inversion algorithm realized in frequency domain is presented.Through analysis of electromagnetic wave propagation in multi-layermedia,generalized reflection coefficients are derived.After dividing the coefficients into pairs followed by decomposing the pairs into even and odd components,a cost function is established to connect reflection coefficients’spectrum and layer’sparameters before inversion is carried out.Some simulations show that accurate inversion results can still be expected when thickness of the layer is less than tuning thickness,which therefore i mproves resolution of thin layers.

Ground Penetrating Radar;thin layer(GPR);spectral inversion algorithm;general reflection coefficients;cost function

waves ofmodelA

1671-5942(2011)04-0154-06

2011-03-16

国家自然科学基金(40774065)

黄忠来,博士研究生,主要研究方向为探地雷达及随机信号处理.E-mail：huangzhonglai@163.com

P631.2

A