三维介质中探地雷达(GPR )波传播逐段迭代射线追踪方法研究和应用
曾昭发1
,高尔根
2
(1.吉林大学地球探测科学与技术学院,吉林长春 130026;2.中国科技大学地球与空间科学系,安徽合肥 230026)
摘要:基于Fer mat 原理和Snell 定理对三维任意界面情况下的两点间GP R 波传播射线追踪问题进行了研究,采用了一个适用于任意界面情况下计算反(折)射点的一阶近似公式,从任意给定的初始路径出发进行逐段迭代计算,给出了三维逐段迭代射线追踪算法的计算格式,并进行了三维介质下射线追踪计算。结果表明计算速度相当快,且其计算精度可以根据需要达到任意要求。在此基础上,根据雷达波传播规律,考虑每一段的能量损失。根据地质雷达的发射子波,最终获得三维介质中雷达波传播和接收记录。并将结果应用于实际剖面拟合中,提高了解释的精度。
关键词:探地雷达;三维结构;逐段迭代;射线追踪中图分类号:V 443.2 文献标识码:A
基金项目:国家自然科学基金项目(40474042)
作者简介:曾昭发(1966-),男,博士,教授,从事工程与环境地球物理教学和研究工作,T el:0431-*******。
The 3-D Modeling by the Segment -by -Segment Iterative
Ray -Tracing Method Study of GPR Wave
ZEN G Zhao -fa 1,GAO Er-gen 2
(1.Colleg e o f GeoExp loration S cience and T echnology ,J ilin Unive rsity ,Chang chun 130026,China ;2.Dep t .of E arth and Sp ace S cience ,S cienc e and T echnology Unive rsity of China ,H ef ei 230026,China )
Abstract :A new kind o f ray -tracing method -seg ment by segm ent iterative ray -tracing method in random interface of GPR(Gr ound Penetrating Radar )w ave is studied in this paper.First,it starts from the fact that the ray satisfies the sam e r ay param eter in the w hole ray -path,reflection or refraction point
in the interface isobtained thr oug h a one -order approx imate form ula w hich fits to random interface .Finally,the segment by seg ment iterative ray-tracing is completed.Num er ical tests based on this me-thod show that this alg orithm is faster than traditio nal shooting method in r andom interface velocity distr ibution and the results are accurate enoug h.At the m ean tim e,travel tim e can be given in the pr ocedure of ray -tracing .It offers a new kind o f efficient r ay -tracing metho d for GPR tomo graphy .Base on the study ,we calculated the ener gy loss of GPR'w ave acco rding tothe GPR w ave's law of pr opag https://www.doczj.com/doc/3613343036.html,ing the w av elet,we got the reco rding of GPR w ave in 3-D m edium.
Key words :GPR ;tracing ;ray path ;to mog raphy
0 引 言
探地雷达已成为工程和环境地球物理领域中的一种重要的方法,解决诸如公路路面层和基底结构、高层建筑基底形态、近地表土壤层结构、地下水位面
和岩石分层、地下水污染评价等问题,而这些问题大都可以归结为地下三维曲面的探测问题。显然,研究
雷达波在这类问题的介质中传播,将为提高探测的
效果和解释的准确性具有重要的意义。对于三维介质中雷达波的传播问题研究方法很多,如有限元和
第35卷 专辑吉林大学学报(地球科学版)
V ol.35 Sup.
2005年7月
Jo ur nal o f Jilin U niver sity (Ear th Science Edition)July 2005
时间域有限差分(FDTD )方法,但这些方法计算量较大,运算速度较慢。探地雷达的探测通常需要实时
地得到解释的结果和尽可能快地对地下做出评价。因而计算速度快的射线追踪方法成为探地雷达快速解释的最具潜力方法。
射线追踪方法已广泛应用于波传播问题(正演问题)研究一类主要方法。射线追踪方法以计算速度快、结果直观,在层析成像技术、偏移、反演及模型试算中都占有很重要的地位。射线追踪方法在地震勘探和研究中发展很快,并取得了实际的进展(Julian 等,1977;Chander ,1977;Pereyr a 等,1980;Keller
等,1983;Um ,J 等,1987;Hanyg a ,1988;Sambridg e 等,1990;Farra ,1992;马争鸣等,1991;黄联捷等,1992)。由于地质雷达理论和地震勘探理论的相似性(Ursin ,1983),近年也有学者将射线追踪方法用于地质雷达波的传播和正演研究,如Cai(1995)应用射线追踪方法进行了二维介质中雷达波的传播与模拟研究,并进行了二维剖面反演研究(Cai ,1999)。但三维介质中的雷达波正反演研究还是一个重要的研究方向。近年来,探地雷达的多道测量仪器技术和三维测量的发展,进行三维介质的快速拟合计算和反演已成为迫切需要。到目前为此,尚未见快速而准确的三维地质介质中雷达波传播和正演模拟研究成果。
本文借鉴了地震波的三维逐段迭代射线追踪算法,考虑了雷达波的动力学特征,实现了探地雷达波在三维介质中的正演模拟,获得了三维逐段迭代射线追踪算法的计算格式,并进行了雷达波在三维介质下射线追踪计算结果。模型计算表明:它的计算速度相当快,且其计算精度可以根据需要达到任意要求。在此基础上,根据雷达波传播规律,考虑每一段和界面附近的能量损失。根据地质雷达的发射子波,获得三维介质中雷达波传播和接收记录。开发了计算的可视化计算软件,并在资料的解释中进行应用,取得很好的效果。
1 方法原理
地质雷达方法是地球物理方法中的一种高分辨率、高效率、实时的探测方法。地质雷达方法的基本原理是利用高频电磁波(1M ～1GHz),以脉冲形式通过发射天线被定向地送入地下。雷达波在地下介质中传播时,当遇到存在电性差异的地下地层或目
标体时,电磁波便发生反射,返回地面后由接收天线所接收。在对接收天线所接收到的雷达波进行处理
和分析的基础上,根据接收到的雷达波波形、强度、时间等推断地下介质空间位置、结构、电性质及几何
形态,从而达到对地下地层或目标体的探测。可见,雷达波在介质中传播的运动学特征类似于地震勘探理论中的地震波。
根据射线追踪理论,同一条射线路径上满足相同射线参数,实际隐含着任意连续三点满足此参数,而三点间的具体形式为Snell 定律。这样,可以从射线的任一端出发,依据Snell 定律,求出中间点,然后以一点移动为步长,顺次进行下去,直到另一端点。当整条路径上的总校正量满足精度要求时,认为最
后一次的追踪结果即为射线路径(图1)。图中v 和l 分别为雷达波传播的速度和传播的路径。图1 雷达波在三维介质中透射示意图
Fig .1 The transmission sketch map of GPR wave in 3-D
120 吉林大学学报(地球科学版) 第35卷
1.1 一阶近似公式的推导
以透射波为例,对于反射波情况,原理相同。如图1(b)所示:
t =l v +l
-v -=
(x -x )2+(y -y )2+(z -z )2
v +
(x --x )2+(y --y )2+(z --z )2
v 。
(1)
对上式求导,根据Ferm at 原理可得
(x -x )+(z -z )z ′x
v l
+(x -x -)+(z -z -)z ′y
v -l
-=0,(y -y )+(z -z )z ′y
v l
+(y -y -)+(z -z -)z ′y
v -l
-=0。(2) 对(2)式在初始路径点(x *,y *,z *)进行一阶
Tay lor 展开,令
x *-x
=a ,y *-y =b ,z *-z =c ,x *-x -= ,y *-y -=v ,z *
-z
-= ,p =a +cz ′x l ,q =b +cz ′y
l ,s =
+ z ′x l -,t =v + z ′y l
-, 1=1+cz ″x +z ′2x , 1=1+ z ″x +z ′2
x , 2=1+cz ″y +z ′2y , 2=1+ z ″y +z ′2y ,
3=cz ″x y +z ′x z ′y , 3= z ″xy +z ′x z ′y , 可推得
AX =B 。
(3)
这里A 为2×2矩阵,X 为校正量向量,B 为一与
初始路径有关的向量,其具体形式为:
A 11=p sv -l +v -l - 1+p sv l -+v l 1A 12=p tv -l +v -l - 3+qsv l -+v l 3A 21=qsv -l +v -l - 3+p tv l -+v l 3A 22=qtv -l +v -l - 2+qtv l -+v l
2B 1=-(p v -+sv )l l -,
B 2=-(qv -+tv )l l -X =( x , y )T
。
(4)
通过求解(3)式,解得 x 和 y 值后,利用(x *+ x ,y *+ y ,z (x *+ x ,y *+ y ))代替原来的(x *,y *,z *),继续重复上述过程,直至满足射线追踪精度为止。
上述推导是对于透射波情况下的结果,在反射波情况下,可按上述办法推导出射线路径的矫正公式如下:
A ~X =
B ~
。(5)
这里A ~为2×2矩阵,X 为校正量向量,B ~为一与初始路径有关的向量,其具体形式为:
A ~11=v (p sl +l - 1+p sl -+l
1)A ~12=v (p tl +l - 3+qsl -+l 3)A ~21=v (qsl +l - 3+p tl -+l 3)A ~22=v (qtl +l - 2+qtl -+l
2)B ~1=-(p +s )vl l -,B ~2=-(q +t )vl l -X =( x , y )
T
。
(6)
1.2 介质特殊分布情况
以上给出的一阶近似公式适用于任意界面情况下,对于介质特殊分布时,可进一步简化。
(1)当介质为水平层状时,z ′x =z ′y =0,假设分界面为z =z 0,对于透射波,此时(4)式可简化为
A 11=a (v -l -+v l )+v -l -+v l A 12=av v -l -+b v l
A 21=b v -l -+av v
l
A 22=bv (v -l -+v l )+v -l -+v l
B 1=-(av -l -+ v l ),B 2=-(bv -l -+v v l ) 。
(7) 对于反射波,只须将上式中的所有v 和v 替换成
v 即可。
当介质分界面为倾斜层时,z ′x =M ,z ′y =N ,M 与N 为常数。此时对于透射波情况,式(4)可简化为
A 11=(a +cM )( + M )(v -l -+v
l
)+(1+
M 2)(v -l -+v
l ),A 12=(a +cM )(v + N )v
-l
-+(b +cN )( +
v M )v
l
+M N (v -l -+v l ),
A 21=(b +cN )( + M )v
-l
-+(a +cM )(v +
N )v
l
+M N (v -l -+v l ),
A 22=(b +cN )(v + N )(v -l -+v
l
)+(1+
N 2)(v -l -+v
l ),121专 辑 曾昭发,高尔根:三维介质中探地雷达(GP R)波传播逐段迭代射线追踪方法研究和应用
B 1=-(a +cM )v -l --( + M )v l ,
B 2=-(b +cN )v -l --(v + N )v l 。而对于反射波,同上。
2 收敛性问题
上述公式推导中,曾采用一阶Tay lor 展开办法进行逐段迭代求取射线路径,这种办法对于初始路
径状态不太好时,有可能造成迭代失败,当然,对于探地雷达来说,由于发射天线与接收天线,按照一定规律排列,采用上次射线路径作为下一次初始路径可以克服迭代失败,但就研究方法本身来说,研究稳定算法具有重要意义。
研究中发现,当采用不完全一阶Tay lor 展开时,可获得一正定方程组,其形式与(3)式一致,各分量具体形式如下:
A 11=(v -l -+v l )(1+z ′2x )A 12=(v -l -+v l )z ′x z ′y A 21=A 12
A 22=(v -l -+v l )(1+z ′2y )
B 1=-(p v -+sv )l l -B 2=-(qv -+tv )l l
-
。
(8)
利用(4)式各元素进行射线追踪计算时,总能计算出
射线路径。
3 雷达波的衰减问题
雷达波在地下介质的传播中,根据雷达波传播
的动力学特性,雷达波能量具有较强的衰减。根据雷达波传播的功率原理,雷达波在通过j 层介质时的振幅为
A p =S D s D R R
G i G 0
j T i e - J l J 。可见,反射雷达波的振幅(或能量)与如下因素有关(Cai,1995)。
(1)雷达波的有效振幅(S );(2)雷达波发射天线的定向因子D s ;(3)雷达波入射平面的几何扩散因子G i ;(4)雷达波入射平面的几何扩散因子G 0;(5)雷达波在界面的反射系数R ;(6)雷达波在界面的透射系数T ;(7)雷达波在介质中的衰减因子 ;(8)雷达波接收天线的定向因子D R ;(9)雷达波在介质中传播的传播距离l j 。
其中(1)～(8)的参数,在分层均匀介质中,当初始条件确定后,这些参数均为常数。而雷达波的传播
距离l j 可以按以上射线追踪的方法获得。
4 数值模拟计算
本文在典型数据模型之上利用逐段迭代射线追踪方法进行数字模拟计算,并可使我们从中体会到
该方法的计算速度,精度与可行性。
图2是一个三层任意界面的地质雷达的透射波射线追踪介质模型和波形图,发射天线位于T(200,200,0),接收天线位于R(600,600,0)。
第一层底界的界面函数为z 1=200-30sin (x /100)+40cos(y /200)cm,速度为5cm /ns;
第二层底界的界面函数为z 2=500-[(x -500)/50]2
+[(y -400)/50]2
(cm ),速度为10cm/ns ;
第三层底界的界面函数为z 3=700-30cos (x /100)+40sin(y /200)(cm ),速度为12cm /ns;下部介质速度为15cm /ns 。
图2(a )为射线和模型图。图2(b )为波形图,图2(b)中1为射线在无耗介质中传播,2为考虑射线在传播过程中由于电导率的影响而产生电磁波能量的衰减,而模型中介质的电阻率自上而下分别为100 ?m 、
200 ?m 、500 ?m 、1000 ?m 。5 应用实例
利用探地雷达对吉林省长春市某大厦场地进行探测,目的是查明该场地的基底分布,以便合理地确定桩基础的类型和确定场地30m 以内物性层的分布,及是否存在物性层的明显变化,并以此推断断层的存在及其分布规律。测量采用美国GSSI 公司生产
的SIR —2型地质雷达。完成了9条测线中南北向分布的NS2测线,长度为120m 。如图3。图中ZK2为
钻孔,钻探深度为30m 左右。进行三维正演模拟的初始资料根据钻孔资料获得。
根据三维正演模拟,首先根据钻孔资料获得了地下各层介质中雷达波的传播速度,以该资料为基础,拟合层面不同部位的深度资料,解释获得6个较清楚的地质层和一个地下水位面。根据地质雷达测量资料和层面分布情况,可以推断该场地不存在大的、延伸较长断层或破碎带。
由于测量过程中仪器的增益系数未知和多次波的干扰,正演拟合中没有获得雷达波传播的真振幅的波形。
122 吉林大学学报(地球科学版) 第35卷
图2 三层三维任意界面的介质模型的透射波射线追踪(a)和波形图(b)
Fig.2 The transmission ray tracing in3layer3-D medium(a)and wave-form(b
)
图3 NS2剖面地质雷达测量结果图(ZK2-为钻孔) Fig.3 GPR prof ile map at a building f ield in Changchun
6 结论与讨论
模型的数值模拟计算和实际应用表明,本文提出的三维结构下雷达波传播的逐段迭代射线追踪算法,其基本思想明确,增量以显式给出,计算速度快,追踪计算的射线精度高,对于给定模型,可以追踪到真实路径的任意精度下的逼近值。
研究表明,该射线追踪方法不仅能模拟三维曲面中雷达波传播问题,而且将为地质雷达资料的三维反演解释提供了基础。由于二维问题为三维问题的一个特例,因而该软件不仅能对三维资料的解释,也能对地质雷达的二维剖面资料进行解释。
本文提出的算法可模拟任意地质雷达的天线组合,对地质雷达方法应用的可行性研究和地质雷达的数据采集参数的确定具有重要的作用。
这里需要说明的是,对于介质分布比较特殊,或两点之间存在多条射线时,利用本方法总可以追踪到初始射线附近的一条射线路径。如果要追踪一条时间最短射线路径,就需要结合其它方法或通过人机对话形式对射线初始路径进行全局搜索。
还有,本方法的计算速度主要取决于中间点的计算速度快慢。本文中的中间点求取是通过一阶T ay ler展开来逐段迭代计算的,其计算速度比用优化发等方法求解方程(3)快得多,并且计算精度高。本方法计算速度快的主要原因可能是:它仅用加,减,乘,除等四则运算来进行射线追踪计算,从而避免了像打靶法等方法那样一方面需要反复搜索射线路径与界面的交点,另一方面需要多次计算三角函数值或求解方程组。
参考文献:
[1] Chander R.On tr acing seismic ray s w it h specified end
points in la yer s of constant v elocity and pla ne inter face
[J].G eophys.Pr osp.,1997,25:120-124.
[2] Far r a V.Bending method r evisited:a Hamilto nian
appr oach[J].G eophys.J.I nt.,1992,109:138-150. [3] Fa rr a V.R ay tr acing in complex media[J].J.A ppl.
Geo phy s.,1993,30:55-73.
[4] 高尔根,徐果明,赵炎.一种任意界面的逐段迭代射线
追踪方法[J].石油地球物理勘探,1998,33:54-60.
(下转136页)
123
专 辑 曾昭发,高尔根:三维介质中探地雷达(GP R)波传播逐段迭代射线追踪方法研究和应用
相一致,200～370m处为煤柱,其形态更好反映了煤柱体的存在。370～600m为低阻反映,为采空区的位置及范围。
以上三种方法对于采空区,塌陷区及残留煤柱位置、范围的判断是一致的。
2 结 论
通过以上分析,对于营城煤矿采空区、塌陷区的地球物理特征可归为以下几点:
1)采空区的电阻率为低阻,氡浓度为高正异常。
2)塌陷区的电阻率上部为高阻,下部为低阻;氡浓度曲线为锯齿状,波动较大,浓度值为中等,50～100Bq/m3。
3)残留煤柱电阻率为高阻,氡浓度为低值。
由于地质资料的缺乏,所以解释推断没有已知资料作为验证,但综合地球物理勘探方法可以互相佐证,综合判断,使解释与推断更趋合理。
在野外数据测量过程中得到了徐波等研究生及2001级本科生的帮助,在此表示感谢!
参考文献:
[1] 韩许恒,郁春霞.氡射气探测在采空区勘察中的应用
[J].工程勘察,1996,(5):63-72.
[2] 寇绳武,李克祥,郭舜,等.高密度电阻率法探测洞穴、
采空区的效果分析[J].工程勘察,1994,(6):61-65. [3] 任基林,邓金灿,苏亚汝,等.瞬变电磁法在探明某矿
区采空区方面的应用[J].矿产与地质,2001,15(3):
39-61.
[4] 朱叶.综合物探方法在煤矿采空区的应用[J].河北煤
炭,2001,(3):39-40.
[5] 郭崇光,刘君.综合电法在采空区探测中的应用[J].
科技情报开发与经济,2004,14(2):142-143.
[6] 段鸿杰,唐岱茂,曹为民.测氡技术圈定采空区影响边
界的应用[J].华北地质矿产杂志,1999,14(1):71-76.
[7] 郭崇光,李敬宇,刘君.氡气测量在山西采空区探测中
的应用[J].科技情报开发与经济,2004,14(1):180-
181.
(上接123页)
[5] 高尔根,徐果明,蒋先艺,等.三维结构下逐段迭代射
线追踪方法[J].石油地球物理勘探,2002,37(1):11-
16.
[6] 高尔根,席道瑛,曾昭发.三维结构下全路径迭代射线
追踪方法[J].岩石力学与工程学报,2002,21(11):
1610-1614.
[7] Sambr idge M S.Kennett,B.L.N.,Boundary value
r ay tr acing in het ero geneous medium:a simple and
v ersatile algo rit hm[J].G eophys.J.Int.,1991,101:
157-168.
[8] U m J,T hurber.C H.A fast alg or ithm fo r tw o-point
seismic r ay tr acing[J].Bull.Seis.Soc.A m.,1987,
77:972-986.
[9] 李大心.探地雷达方法及应用[M].北京:地质出版
社,1994.[10] Cai J,M cM echan G A.R ay-based synthesis of
bistatic gr ound-penetr ating r adar pr ofile[J].
Geo phy sics,1995,60:87-96.
[11] Cai J,M cM echan G A.2-D R ay-based to mog ra phy
fo r velocit y,layer shape,a nd at tenuat ion fro m G PR
data[J].G eophysics,1999,64:1579-1593.
[12] U r sin B.R eview o f elastic a nd electr o magnetic w av e
pr o pag atio n in hor izontally layer ed media[J].
Geo phy sics,1983,48:1063-1081.
[13] 马争鸣,李衍达.二步法射线追踪[J].地球物理学
报,1991,34:501-508.
[14] 黄联捷,李幼铭,吴如山.用于图像重建的波前法射
线追踪[J].地球物理学报,1992,35(2):223-233.
136 吉林大学学报(地球科学版) 第35卷