探地雷达各向异性介质有限差分偏移线性变换
戴前伟,张彬,冯德山
(中南大学 地球科学与信息物理工程学院,湖南 长沙,410083)
摘要:将地震偏移技术中的线性变换差分偏移应用于探地雷达的数据偏移,在讨论各向同性介质线性变换差分偏移算法的基础上,研究各向异性介质偏移归位情况来解决该算法对速度的高敏感性问题,编写探地雷达线性变换差分偏移算法的Matlat程序;通过对线性变换差分方程推导,得到一系列三对角线性方程组;在求解这些特殊方程组过程中,通过调用Matlat程序中定义的Pursue函数来循环求解这一系列大型的三对角线性方程组,并将计算结果进行储存供后续迭代使用,最后将满足终止条件的偏移成像矩阵输出并可视化,运用该算法分别对各向异性的正演雷达数据和实测雷达数据进行处理。研究结果表明:该算法在对各向异性介质中存在多次反射质点和主要反射层的地质问题进行偏移归位时具有明显优势。
关键词:探地雷达;有限差分偏移;线性变换;各向异性介质
中图分类号:P631 文献标志码:A 文章编号:1672-7207(2012)05-1814-07
Linear transformation of finite difference migration of GPR in anisotropic medium
DAI Qian-wei, ZHANG Bin, FENG De-shan
(School of Info-Physics and Geomatics Engineering, Central South University, Changsha 410083, China)
Abstract: The linear transformation of finite difference migration method in reflection seismology was applied to improve the GPR (ground penetrating radar) migration technology. Based on the fact that the migration effect is highly sensitive to the medium velocity in isotropic condition, the anisotropic case was extended when the theory of the traditional finite difference migration was established. In order to solve the velocity restrictions of the traditional algorithm, the detailed migration process was presented as well as a part of numerical calculation process using Matlab. Through solving circular tri-diagonal equation systems, the entire solve process was implemented by the functions of Pursue, which is defined in the Matlat program, and then this results were compared to those of the traditional migration. The results indicate that the improved migration algorithm of GPR data migration processing technology is more flexible to solve multiple reflection particles and a main reflector layers when dealing with anisotropic geological medium.
Key words: ground penetrating radar (GPR); finite difference migration; linearly transformation; anisotropic medium
在速度和构造均较复杂的地区,利用有限差分波动方程进行探地雷达波场的偏移归位[1-2]。以往研究有限差分偏移方法是沿深度方向逐层延拓进行,为了消除波场对深度的高阶偏导,采用延时坐标变换或浮动坐标变换[3],在小倾角地层反射前提下略去波场沿深度方向的高阶偏导[4-7],求解零时刻各深度上波场值来解决该边值问题。Li等[8-9]提出的线性变换有限差分偏移方法突破了传统有限差分偏移倾角限制,实现了波场的精确归位。李广场[10]运用线性变换的方法通过精确推导得到均匀介质下的雷达波场外推差分方程,实现了探地雷达各向同性介质的全倾角偏移。在此,本文作者对探地雷达各向同性介质下线性变换差分偏移算法进行优化,讨论各向异性介质条件下该算法的运用情况,利用Matlab软件[11-15]编制相应的迭代及偏移程序,给出详细的求解过程,并通过对比可视化的偏移结果,体现优化的算法在解决探地雷达各向异性介质偏移问题时的优势。
1 各向异性介质线性变换差分偏移原理
1.1 线性变换
探地雷达高频电磁波在二维介质传播的麦克斯韦方程可以写成:
(1)
对于均匀介质,将线性变换代入方程(1)中得:
![](/web/fileinfo/upload/magazine/12230/299905/image004.gif)
没有略去任何项时即得与15°偏移方程形式相同的精确方程:
(2)
研究精确方程(2)与传统的15°近似偏移方程可知:在寻求电磁波遵循的上行波方程中,实现波场精确归位的关键并不是线性变换,而是该变换采用的波场延拓方向,即15°偏移方程是沿深度方向延拓;而线性变换差分偏移则是取上行波波前方向为延拓方向,最终实现反射界面的精确偏移归位。
对于各向异性变速介质,在原算法基础上进行推广,将速度参数
其代入方程(1)得:
(3)
定义
(其中,
为参考
速度,nx为雷达记录道编号),代入波动方程得:
(4)
定义
(其中,
为伪时间深度),有
,则在
中,方程(4)变为:
![](/web/fileinfo/upload/magazine/12230/299905/image026.gif)
(5)
同样引入线性变换:
(6)
则方程(5)变为:
![](/web/fileinfo/upload/magazine/12230/299905/image032.gif)
![](/web/fileinfo/upload/magazine/12230/299905/image034.gif)
![](/web/fileinfo/upload/magazine/12230/299905/image036.gif)
(7)
式中:
;
。这样,在新的坐标系
中,方程(7)为各向异性介质中线性变换差分偏移方程的表达式。这里重点研究式(7)右端第1项对雷达偏移效果带来的影响,将其称为混合因子项e。分析该项可知:它与速度的纵向梯度除以速度再乘以波场的纵向梯度的结果成正比。
1.2 差分格式
在传统有限差分偏移中,
为雷达地表数据,
为雷达偏移剖面,上行波边界条件为
,如图1所示。而经线性变换后的新坐标系
中,雷达地表数据矩阵为
,雷达偏移剖面数据矩阵为
,边界条件则为
,即在
平面上,线性变换差分格式如图2所示。
![](/web/fileinfo/upload/magazine/12230/299905/image062.jpg)
图1 传统有限差分偏移的差分网格
Fig.1 Differential grid of traditional finite difference migration method
![](/web/fileinfo/upload/magazine/12230/299905/image064.jpg)
图2 线性变换差分偏移的差分网格
Fig.2 Differential grid of linear transformation of finite difference migration method
在线性变换后的
坐标下,地表数据矩阵沿
分布。在上行波边界条件约束下,向上或向左进行,将
平面上的雷达地表数据矩阵从t轴向z轴外推直至成像,从数值稳定性出发,采用隐式差分格式,在
坐标系中,方程(2)的线性变换Crank-Nicolson差分格式为:
![](/web/fileinfo/upload/magazine/12230/299905/image071.gif)
(8)
其中:
;P为沿
轴分布的波场向量;I和T分别为以(0,1,0)和(-1, 2, -1)为主元素的三对角矩阵;
为网格j时间、l+1深度处的波场向量,该项由向量
,
和
给出,即与15°方程差分格式不同,经线性变换的差分格式为四点差分格式,如图3和图4所示。
![](/web/fileinfo/upload/magazine/12230/299905/image087.jpg)
图3 (x, z, t)中的五点差分格式
Fig.3 Five-point difference format in (x, z, t)
![](/web/fileinfo/upload/magazine/12230/299905/image089.jpg)
图4
中的四点差分格式
Fig.4 Four-point difference format in ![](/web/fileinfo/upload/magazine/12230/299905/image090.gif)
方程(8)写成矩阵向量形式为:
(9)
未知量S的分量
,已知量B的分量
,系数矩阵为:
![](/web/fileinfo/upload/magazine/12230/299905/image098.gif)
![](/web/fileinfo/upload/magazine/12230/299905/image100.gif)
矩阵在Matlab程序迭代初始阶段,置l=0,j=T,
由雷达地表数据得知,
和
为0。在重复的求解过程中,设置j为外循环,l为内循环,单元格每向上或向左移动1步,就要调用1次计算三对角线性方程组[16-18]的pursue(a, b, c, f1)函数(其中:a,b和c分别为系数矩阵的三对角元素;f1为右端常数项)。求得沿
轴分布的波场向量,循环迭代后求得vt平面上的偏移剖面数据矩阵。
2 Matlab实现流程
2.1 算法步骤
在各向同性介质中,线性变换差分偏移的计算程序实现过程与15°有限差分偏移程序实现过程相似,只是延拓方向不同。15°差分方程偏移按深度延拓,而线性变换算法则是沿上行波波前方向进行延拓,该算法的实现步骤如下。
(1) 设置网格参数迭代步长:
,
和
及网格节点数。
(2) 输入探地雷达时间剖面aryin。
(3) 输入与网格行数量相等的速度矩阵v。
(4) 计算
,并计算线性方程组的系数矩阵参数a,b和c。
(5) 保持j=T不变,设置l内循环,根据上行波第一类边界条件
,从探地雷达地表数据矩阵
开始,先置l=0,j=T即得差分方程:
![](/web/fileinfo/upload/magazine/12230/299905/image121.gif)
调用M程序中定义的pursue(a, b, c, f1)函数求解在l=0,j=T时的线性方程组,求得
。
(1) 在
方向上逐步增加
,即分别置l=
,l=2
…l=j,重复步骤(5)的所有计算,求得各深度对于的波场值
(l=
, 2
, …, j)。
(2) 进入j外循环,在
方向上逐步减少
,即分别置j=T-
,T-2
,…,重复上述步骤(5)和(6)的所有计算,求得各时刻、各深度的波场
,
,…,
。
(3) 在l内循环迭代过程中,即当差分格式分别向上进行外推时,在给定的时间长度内进行,将第j时间层各深度的波场矩阵储存于1个三维数组X的第j层,将三维数组X中l=j的向量组成数据矩阵输出即得偏移剖面
。
2.2 各向异性介质线性变换有限差分偏移算法的Matlab实现流程
如前所述,在(t1,t2)平面内采用四点差分格式,不需要另外增加边界条件。对于输入的速度参数矩阵V,首先需将其在
方向上分布的各个向量转化为
,然后,对
在
坐标中求其差分格式值。具体Matlab程序流程见图5。
![](/web/fileinfo/upload/magazine/12230/299905/image144.jpg)
图5 线性变换有限差分偏移算法Matlab程序流程图
Fig.5 Matlab program flow chart of linear transformation of finite difference migration algorithm
函数Pursue(a,b,c,f)的部分代码如下:
function x=Pursue(a,b,c,f)
…
for i=2:n-1
b(i)=b(i)-a(i)*c(i-1);
c(i)=c(i)/b(i);
f(i)=(f(i)-a(i)*f(i-1))/b(i);
end
f(n)=(f(n)-a(n)*f(n-1))/(b(n)-a(n)*c(n-1));
for i=n-1: -1:1
f(i)=f(i)-c(i)*f(i+1);
end
x=f;
…
在设置j外循环、l内循环步骤中,部分数值的计算过程如下:
------j外循环、l内循环------
…
for j=T:-dt:1
if j==T
%迭代初始时将地表数据第一行作为右端常数项
f1=(aryin(1,:))';
for l=0:dz:j
%快速追赶法计算T时间层、各深度波场值p(i, l, T-1) ,并将其存于三维数组的第T层
x_T(:,l+1,T)= pursue(a,b,c,f1)
end
else
x2=(aryin(T-j,:))';f1=(aryin(T+1-j,:))';
for l=0:dz:j
f2=x_T(:,l+1,j+1);
x_T(:,l+1,j)=pursue(a,b,c,(f1));
f1=x_T(:,l+1,j);
x2=x_T(:,l+1,j+1);
end
end
end
…
3 偏移实例
3.1 各向异性介质地电模型
选取复杂的各向异性地电模型,模型为3层变速介质:第1层设置为干黏土,第2层为混凝土,第3层为湿砂层。在混泥土层内设置球体,球体下方为1个具有不同倾角界面的主体反射层,球体半径为0.05 m,埋深为0.20 m,各介质层速度平缓变化。粗略估计混凝土层与球体交界处速度由0.12 m/ns突变为0.27 m/ns,各层的相对介电常数ε分别为:ε1=12,ε2=6.0,ε3=25;球体的相对介电常数为ε4=2.0,各介质层的电导率分别为:σ1=2.0×10-4 S/m, σ2=0.04 S/m,σ3=0.006 9 S/m,球体的电导率σ4=10-6 S/m。
![](/web/fileinfo/upload/magazine/12230/299905/image146.jpg)
图6 各向异性介质层状地电模型
Fig.6 Layer geoelectric model with anisotropic medium
图7(a)所示为该模型的正演结果。从图7(a)可见:球体处产生强烈绕射双曲线,不同倾角的斜面交界处产生强烈的绕射波,球体形状无法辨明,不同倾角的界面也只能分辨出波浪形轮廓,且斜面的倾角大幅度变小,正演模拟结果发生了严重的畸变。图7(b)所示为15°有限差分法偏移的结果。从图7(b)可见:球体处绕射波的收敛归位能力较强,而不同倾角分界面的偏移归位效果较弱,表明在陡倾角界面的各向异性介质中,15°算法对倾斜反射界面的偏移归位能力较差。图7(c)所示为利用Kirchhoff积分法进行偏移处理的效果。从图7(c)可见:球体的归位能力比15°算法的强,倾斜界面的偏移结果也较差。同样的情况也出现在频率波数域偏移上,其结果如图7(d)所示,倾斜界面的轮廓完全畸变,只是该算法处理噪声的能力比前者强。图7(e)所示为带混合因子项的线性变换差分偏移算法处理的结果。从图7(e)可见:不论缓倾角还是陡倾角成像都清晰,虽然球体界面处产生的绕射波收敛效果一般,但较易辨明其形态及位置,整体偏移的效果比较理想。图7(f)所示为不带混合因子项的算法偏移的结果,可见倾斜分界面的连续性很好且偏移成像清晰,球体周围无反射波产生。这是由于略去了混合因子项,当波场外推至速度变化剧烈的球体界面处时无反射波产生。
3.2 实测探地雷达数据
图8所示为武广高速铁路某隧道管线探测的雷达剖面。为了查明地下管道的分布及定位,采用中心频率为900 MHz的屏蔽天线,每道扫描采样点数设为512,采样间隔为0.15 ns,采集时窗为15 ns。从图8可以看出:该段存在诸多双曲线绕射波,还有2处较明显的反射带,其中绕射双曲线与倾斜的反射带交错分布,给管道的定位及解释带来很大困难。
![](/web/fileinfo/upload/magazine/12230/299905/image147.jpg)
图7 不同偏移算法对探地雷达各向异性介质层状地电模型处理的结果
Fig.7 Migrated GPR results of layer geoelectric model using different migration algorithm
![](/web/fileinfo/upload/magazine/12230/299905/image149.jpg)
图8 探地雷达管道实测剖面图
Fig.8 Real GPR section of pipe detection
图9(a)和9(b)所示分别为带混合因子项和不带混合因子项的各向异性介质线性变换差分偏移算法处理后的结果。对比图9(a)和9(b)可见:带混合因子项的算法在介质速度突变处产生较强烈的反射波,而略去了混合因子项的算法在该处则减弱了这种效应,图示居中矩形框的一道双曲线绕射波处表现尤其明显,二者均可显示2条几乎平行的反射带,推断为管道的上下管壁;图9(a)中诸多收敛的强反射质点推断为钢筋网,经证实该推断符合实际情况。
![](/web/fileinfo/upload/magazine/12230/299905/image150.jpg)
图9 线性变换差分偏移算法处理后的探地雷达实测剖面结果
Fig.9 Migrated real GPR sections using linear transformed finite difference migration method
4 结论
(1) 将各向异性介质线性变换有限差分偏移技术应用于探地雷达数据处理,编制了求解大型三对角矩阵方程组的Matlab程序和探地雷达各向异性介质线性变换有限差分Matlab偏移程序。
(2) 当横向介质速度变化剧烈时,混合因子项可提高对速度矩阵中横向速度剧变处(多次反射质点)的偏移精度;而在对主要反射体进行偏移归位时,为了减弱这种在速度突变处产生的透射和反射,必须忽略该混合因子项。
参考文献:
[1] Bednar J B. A brief history of seismic migration[J]. Geophysics, 2005, 70: 3MJ-20MJ.
[2] 周彬忠, 包吉山. 复杂地区地震波场的正演模拟与反时偏移[J]. 石油地球物理勘探, 1988, 23(3): 257-262.
ZHOU Bing-zhong, BAO Ji-shan. The forward modeling and reverse time migration of seismic wave field in complex medium[J]. Oil Geophysical Prospecting, 1988, 23(3): 257-262.
[3] Claerbout J F. Coarse grid calculations of waves in inhomogeneous media with application to delineation of complicated seismic structure[J]. Geophysics, 1970, 35(3): 407-418.
[4] Stolt R H. Migration by Fourier transforms[J]. Geophysics, 1978, 43(1): 23-48.
[5] 马在田. 高阶方程偏移的分裂算法[J]. 地球物理学报, 1983, 26(4): 377-388.
MA Zai-tian. High order equation migration of splitting algorithm[J]. Chinese Journal of Geophysics, 1983, 26(4): 377-388.
[6] 张关泉. 利用低阶偏微分方程组的大倾角差分偏移[J]. 地球物理学报, 1986, 29(3): 273-282.
ZHANG Guan-quan. Steep dip finite-difference migration using the system of lower-order partial differential equations[J]. Chinese Journal of Geophysics, 1986, 29(3): 273-282.
[7] 董良国, 马在田, 曹景忠, 等. 一阶弹性波方程交错网格高阶差分解法[J]. 地球物理学报, 2000, 43(3): 411-419.
DONG Liang-guo, MA Zai-tian, CAO Jin-zhong, et al. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(3): 411-419.
[8] LI Zhi-ming. Wave field extrapolation by the linearly transformed wave equation operator[J]. Geophysics, 1986, 51(8): 38-41.
[9] LI Zhi-ming. Imaging steep-dip reflection by the linearly transformed wave equation method[C]//Proceeding of SEG 56th Annual Meeting. Tulsa, 1986: 47-49.
[10] 李广场. 有限差分法探地雷达波动方程偏移[D]. 天津: 河海大学土木学院, 2004: 33-36.
LI Guan-chang. The wave equation migration of ground penetrating radar by the method of finite difference[D]. Hohai University. College of Civil Engineering, 2004: 33-36.
[11] 张志涌, 杨祖樱. 掌握和精通MATLAB[M]. 北京: 北京航空航天大学出版社, 2000: 162-183.
ZHANG Zhi-yong, YANG Zu-yin. Master and master in MATLAB[M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 2000: 162-183.
[12] 熊彬, 阮百尧. MATLAB在有限差分法中的应用[J]. 桂林工学院学报, 2001, 21(2): 104-109.
XIONG Bing, RUAN Bai-yao. The application of MATLAB to finite difference method[J]. Journal of Guilin Institute of Technology, 2001, 21(2): 104-109.
[13] Chapman S J. MATLAB programming for engineers[M]. 2卷. 北京: 科学出版社, 2002: 81-179.
Chapman S J. MATLAB programming for engineers[M]. 2nd ed. Beijing: Science Press, 2002: 81-179.
[14] 王飞, 裴永祥. 有限差分方法的MATLAB编程[J]. 新疆师范大学学报: 自然科学版, 2003, 22(4): 22-27.
WANG Fei, PEI Yong-xiang. The edit program by MATLAB for finite difference method[J]. Journal of Xiangjiang Normal University: Natural Science Edition, 2003, 22(4): 22-27.
[15] Irving J, Knight R. Numerical modeling of ground-penetrating radar in 2-D using MATLAB[J]. Computers Geosciences, 2006(32): 1246-1258.
[16] 刘晓, 李文强. 追赶法在求解循环和拟循环三对角方程组中的一种推广[J]. 河南师范大学学报: 自然科学版, 2009, 37(1): 13-16.
LIU Xiao, LI Wen-qiang. An extension of chasing method for solving circular and quasi-circular tri-diagonal systems[J]. Journal of Henan Normal University: Natural Science, 2009, 37(1): 13-16.
[17] 迟利华, 刘杰, 李晓梅. 三对角线性方程组的一种有效并行算法[J]. 计算机学报, 1999, 22(2): 218-221.
CHI Li-hua, LIU Jie, LI Xiao-mei. An effective parallel algorithm for tri-diagonal linear equation[J]. Chinese Journal of Computers, 1999, 22(2): 218-221.
[18] 李青, 王能超. 解循环三对角线性方程组的追赶法[J]. 小型微型计算机系统, 2002, 23(11): 1393-1395.
LI Qing, WANG Neng-chao. An algorithm for solving circular tri-diagonal systems[J]. Mini-micro Systems, 2002, 23(11): 1393-1395.
(编辑 陈灿华)
收稿日期:2011-07-15;修回日期:2011-09-26
基金项目:国家自然科学基金资助项目(40804027,41074085);教育部博士点新教师专项科研基金资助项目(200805331082);湖南省自然科学基金重点资助项目(09JJ3084)
通信作者:冯德山(1978-),男,湖南祁阳人,博士,副教授,从事探地雷达正、反演及数据处理研究;电话:0731-88836145;E-mail: fengdeshan@126.com