DOI: 10.11817/j.issn.1672-7207.2015.07.037
基于旅行时修正的钻孔雷达层析成像改进
朱自强，彭凌星，密士文
(中南大学 地球科学与信息物理学院，湖南 长沙，410083)
摘要：基于钻孔雷达的波幅采用速度层析成像时，大收发角度雷达波幅因其信噪比低，旅行时提取较困难等问题，钻孔雷达雷达波幅层析成像精度不高。通过计算证明电磁波能量是从发射天线中心点传播至末端，再由接收天线末端传播至其中心。依据此传播路径采用互相关函数对旅行时的提取进行优化，并得到相应的电磁波波速。利用优化后的旅行时与波速得到旅行时的修正值，并根据修正后的旅行时进行速度层析成像，将进行旅行时修正的层析成像图与未进行旅行时修正的层析成像图进行对比。研究结果表明：经过旅行时修正后得到的层析成像结果更加精确，证明对钻孔雷达层析成像的改进方法是可行的。
关键词：钻孔雷达；初至波旅行时；互相关函数；旅行时修正；层析成像
中图分类号：P631.3 文献标志码：A 文章编号：1672-7207(2015)07-2658-07
Improved borehole-GPR tomography based on travel-time correction
ZHU Ziqiang, PENG Lingxing, MI Shiwen
(School of Geoscience and Info-Physics, Central South University, Changsha 410083, China)
Abstract: In the course of borehole GPR tomography processing, the first-arrival time at high transmitter-receive angles is often very difficult to pick by using common tools because of low ratio of signal and noise in the data. The accuracy of borehole radar tomography is not high. It was verified by calculation that the propagation path of electromagnetic wave energy in velocity tomography goes from the center of the transmitter antenna to the end of the transmitter antenna, then to the end of the receiver antenna, and finally to the center of the transmitter antenna. Depending on the true propagation path of electromagnetic wave, the cross-correlation function was used to optimize the extraction of travel time, and then the optimized velocity was obtained. The value of travel-time correction was obtained by calculating the average velocity along each ray and optimized time-picking using cross-correlation, and the optimized velocity tomography was obtained. The results show that compared with the tomography without correction, the optimized tomography is more precise, which verifies the feasibility of the proposed method.
Key words: borehole radar; first-arrival time; cross-correlation function; travel-time correction; tomography
20世纪70年代之后，在地质雷达的基础上，Olhoet[1]为了充分利用钻孔这一已有的通道对地下介质的结构特征等信息进行探测，提出一种新的物探方法即钻孔雷达探测。钻孔雷达在国内的发展起步较晚，到20世纪末国内才出现有关钻孔雷达的研究，但基本上以介绍为主，之后渐渐对其应用和理论进行研究[2]。黄家会等[3]应用钻孔雷达对地下深部灰岩地区的裂隙带和水流通道等进行探测，通过对雷达波幅进行成像对高衰减低速区、低衰减高速区与低衰减低速区等进行区分。刘四新等[4]通过有限差分数值模拟了充水裂缝和矿体在钻孔雷达中的响应特征，证明钻孔雷达在裂缝探测中的可行性。王飞等[5]通过有限差分模拟证明了钻孔雷达单孔反射探测可以有效确定矿体的位置和形态。在钻孔雷达的成像研究方面，李玉喜[6]采用GPR-MAX钻孔雷达进行正演模拟，并通过反演实例讨论了射线覆盖角度、正则化因子等对成像的影响。杨薇等[7]针对解答大型稀疏矩阵的LSQR方法进行了研究，利用钻孔雷达层析成像直观反映了地下介质的速度、吸收系数等参数。在钻孔雷达的成像理论研究方面，Jens等[8-9]结合钻孔雷达直达波和反射波、折射波进行联合成像，提高了成像的精度，能够有效得到分辨地下介质结构。Jacques等[10]在基于时域有限差分法的全波场反演基础上，对2组跨孔雷达实测雷达波幅进行反演，反演结果与实际地质结果基本一致。James等[11]通过对旅行时提取的改进来联合所有雷达波幅进行层析成像，有效降低了层析成像图中大收发角度雷达波幅引起的假象。Bernard等[12]基于Akaike信息准则(AIC)和连续小波变换(CWT)改进了层析成像中旅行时的提取方法，Gokhan等[13]通过程函方程计算雷达波旅行时。还有一些学者在钻孔雷达成像方面采用全波场反演[14-18]。在这些研究中，主要采用层析成像方法进行钻孔雷达成像，这是应用最广泛的方法之一。本文也采用层析成像对钻孔雷达雷达波幅进行处理。针对大收发角度雷达波幅与小角度雷达波幅之间的差异，利用电磁波在天线中传播速度和介质中传播速度之差对旅行时进行修正补偿，从而提高旅行时的准确度以便提高层析成像的效果。
1 旅行时提取优化
钻孔雷达的跨孔探测方式如图1所示。探测方式有ZOP(零偏移距剖面)和MOG(多偏移距雷达探测)共2种，其中MOG方式较常用，对其雷达波幅进行处理的方式也是主要通过层析成像来进行。当所有接收点都接收雷达波幅时，大角度收发雷达波幅因其传播距离较长，会导致电磁波衰减更加强烈。
取1个钻孔雷达实测雷达波幅来对此进行说明。使用瑞典MALA公司的RAMAC地质雷达采集雷达波幅，其天线的主频率为250 MHz。发射天线位于地下11 m处固定不动，接收天线沿接收钻孔从井口匀速下移17 m。对实验雷达波幅进行滤波、归一化振幅以及频谱分析等处理得到的结果如图2所示。随着接收天线深度逐渐增大，收发角度由极大值减小至0°，再逐渐增大。从图2可以看出：随着收发角度的增大，信号的主频也增大，信噪比则慢慢减小。此时，初至波的提取就会变得越来越困难。
图1 跨孔雷达探测示意图
Fig. 1 Sketch map of cross-hole radar detection
图2 钻孔雷达实测雷达波幅处理结果
Fig. 2 Processing results of field data of borehole radar
采用一种新的方法即互相关函数方法对大收发角度雷达波幅提取初至波旅行时。其步骤为：
1) 对雷达波幅进行先期处理，对每道进行归一化处理，然后剔除废道。
2) 将雷达波幅按收发角度进行划分，每5°为1个单元(以收发角度0°~5°为例)。
3) 针对相同收发角度的雷达波幅，对每道雷达波幅与该角度雷达波幅中信噪比最大的一道雷达波幅进行互相关处理，这会使大部分雷达波幅按照合理的方式进行排序，如图3所示。
4) 对这些雷达波幅进行叠加，从而得到收发角为0°~5°时的雷达波幅参考波形，如图4所示。
5) 提取出参考波形的初至波到达时间。
6) 对每道雷达波幅与该角度的参考波形进行互相关处理，从而得到每道雷达波幅的初至波到达时间。
同样，可以得到其他角度的参考波形，与图4类似。这样对每个角度的参考波形与该角度雷达波幅进行互相关处理，提取出初至波到达时间。
图5所示为图2中不同发射点的雷达波幅初至波到达时间提取效果图(取4个雷达波幅为例)，其中黑色粗线对应初至波的旅行时。从图5可以看出：经过互相关处理后的旅行时提取更加精确。这样，就可以得到每道雷达波幅的准确初至波旅行时，从而对雷达波幅进行层析成像。
图3 相同收发角度雷达波幅排序
Fig. 3 Sort sketch map of common-ray-angle data
图4 相同收发角度雷达波幅叠加得到的参考波形
Fig. 4 Reference waveform set by stacking common-ray-angle data
图5 初至波旅行时提取示意图
Fig. 5 Results of automatic picking procedure
2 速度的计算分析
即使可以较准确地提取出大收发角度雷达波幅的初至波旅行时，大收发角度雷达波幅和小收发角度雷达波幅仍存在一定的差异。图6所示为一数值模拟结果，孔间介质的介电常数随机分布在20~30之间，取发射点位于孔中间时接收的雷达波幅成图，这样可以计算出每道雷达波幅的平均波速，如图7所示。
图6 跨孔雷达探测数值模拟结果
Fig. 6 Numerical simulation result of cross-hole radar detection
图7 传播路径平均波速
Fig. 7 Average velocity of propagation path
从图7可以看出：电磁波平均速度随着接收角度的增大而增大然后慢慢减小。大收发角度的雷达波幅明显与小收发角度的雷达波幅波速相差较大，因此，将所有角度的雷达波幅直接进行联合成像会造成成像效果的误差。Peterson[19]阐述了这种类似现象，但是他并没有阐述发生这种现象的原因。James等[20]在2005年通过对数值模拟，证明了不管钻孔内是空气充填还是水充填，电磁波在天线中传播速度通常要比在介质中传播速度快。在此前提下，基于电磁波传播不是在天线中心之间传播而是通过天线和孔间介质传播的假设，本文对旅行时进行计算和分析。
在传统的基于射线理论的速度层析成像中，都是将发射天线和接收天线看成1个点，电磁波沿着2个点之间传播。而在钻孔雷达探测的实际过程中，发射天线和接收天线都是有一定长度的，并且电磁波在天线上传播的速度比在空间介质中传播的速度快，因此，在这个基础上，假设电磁波传播的路径如图8所示。图8中带箭头的线为电磁波传播路径，电磁波首先沿着接收天线传播至端头，然后经由孔间介质传播至接收天线端头，最后到达接收天线中心。
图8 射线传播路径示意图
Fig. 8 Sketch map of propagation path
按照此理论，可以计算出每道雷达波幅平均波速的理论值。因此，设置1个模型，并据此计算其电磁波的波速。仍然以发射点位于中心位置时为例。孔间介质为各向均匀介质，介电常数为25；电磁波在天线中传播的速度为0.11 m/ns；天线长度为0.8 m；2个孔之间的距离为5 m。图9所示为计算的每道雷达波幅的平均波速。从图9可以看出：波速随着角度的增大会逐渐增大然后慢慢减小。但在收发角度较小时，速度基本不变。
但这种曲线的形态过于理想化，是建立在脉冲持续时间相对于电磁波在天线中传播的时间要小得多的基础上。实际上，对于现在所使用的仪器来说不大可能实现。因此，计算时，可以采用10 ns高斯脉冲作为代表进行计算。在不考虑介质的色散情况下，所得计算结果如图10所示。
图9 波速理论值示意图
Fig. 9 Theoretical value of wave velocity
图10 高斯脉冲波速理论值
Fig. 10 Theoretical value of Gauss pulse wave velocity
考虑到实际的钻孔雷达探测过程中，介质会使电磁波发生一定程度色散，这里采用与图9所示模型中一样的参数和脉冲，给予电磁波一定色散，分析色散对不同收发角度雷达波幅波速的影响。Turner等[21]给出色散系数Q的取值范围为2~30。本文取Q为16，所得结果如图11所示。
图11 色散下的高斯脉冲波速理论值
Fig. 11 Theoretical value of Gauss pulse wave velocity in dispersion
从图11可以看出：在基于色散和长脉冲的基础上电磁波波速计算结果与根据模拟雷达波幅计算出的波速分布图非常接近，曲线形状也基本一样。由此可以推断：在钻孔雷达实际探测过程中电磁波的传播路径是按照图8所示的路径进行传播的。层析成像基于射线从天线中心点之间按直线传播的理论进行计算，所以，必须根据实际传播路径对其旅行时进行修正，从而更接近于实际情况。
3 旅行时修正
按照上述计算结果，在对钻孔雷达雷达波幅进行层析成像时，需要对其旅行时进行修正。标准的基于射线理论的速度层析成像方程为
(1)
其中：L为传播路径的长度；s为电磁波在路径中传播的慢度(速度的倒数)；T为旅行时。根据前面的分析，需要对旅行时进行一部分修正。可以将方程改写为
(2)
式中：
为观察到的旅行时与基于射线理论得到的旅行时之间的差值；
为慢度差值，加入到初始模型的速度中以便使旅行时更符合观察到的结果。将式(2)改写为
(3)
其中： M为每道雷达波幅进行修正的权重；
为不同收发角度的旅行时修正值。
在解方程(3)时，因为电磁波的波速较大，所以，旅行时的修正值会很小。故在对修正值曲线进行校正时，不需要对其进行平滑和对称约束处理。需要注意的是：当收发角度为0°(收发天线位于同一水平位置)时，该角度的也为0°。通过解式(3)，可以将所有收发角度的雷达波幅进行速度层析成像。
为了对上述方法进行验证，通过数值模拟对模型层析成像进行对比，图12所示为模型示意图。图12中，异常体形状都较简单，为规则的矩形；背景是介电常数为25的各向均匀介质，左上和右下2个异常体的介电常数为29，其余异常体的介电常数为21；所有介质的电导率均设为1 mS/m；磁导率为真空中的磁导率。尽管图12所示模型较简单，与实际地质情况不太相符，但仅用于对旅行时修正法进行验证，简单的模型可以足够说明问题。
模型的正演雷达波幅采用基于CPML边界条件的有限差分正演[22]。图13(a)所示为没有对旅行时进行修正时根据图12所示模型得出的速度层析成像图；图13(b)所示为对旅行时进行修正之后得到的速度层析成像图。从图13可以看出：进行旅行时修正之后的层析成像图更加精确，能够较清晰地判断出异常体的形状和位置，且能判断出异常体内的波速。据图13(a)能够较清晰地判断异常体的波速、形状和位置，但并不能对模型中所有的异常体进行区分，不能反映出模型左上以及右下2个低速异常体。由此可以判断：在钻孔雷达的速度层析成像中旅行时是需要进行修正的，文中提出的修正方法也可取，得出的层析成像图更加精确，对异常体的反映也更加敏感，能够为判断地下介质的结构特征等提供帮助。
图12 模型示意图
Fig. 12 Sketch map of model
图13 速度层析成像图
Fig. 13 Images of velocity tomographies
4 结论
1) 在针对大收发角度雷达波幅进行初至波提取时，因其具有很低的信噪比，从而较难提取出准确的旅行时。本文对同收发角度的雷达波幅进行排序、叠加得到参考波形，并对参考波形和该收发角度的雷达波幅进行互相关处理，从而得到每道雷达波幅的初至波旅行时。对合成雷达波幅的处理结果显示，该方法对大收发角度雷达波幅的初至波提取具有很好的效果。
2) 以发射点位于钻孔中间位置时为例，本文计算了色散条件下的电磁波波速与收发角度之间的关系，计算结果与模拟结果较吻合，证明了文中提出的电磁波传播路径更符合实际情况，即电磁波首先由发射天线中心点传至发射天线末端，再由发射天线末端传播至接收天线末端最后传播至接收天线中心。因此，在针对钻孔雷达雷达波幅进行层析成像时，需要对旅行时进行一定修正。
3) 加入修正值的成像结果更加精确，更加接近于实际模型。证明了该方法的正确性与可行性。
4) 文中的方法主要针对的是电磁波在孔间介质内波速相差不大的情况，但是，当孔间介质中电磁波的波速产生较大的差异时，无法对旅行时进行修正。此外，层析成像的分辨率受Fresnel带的限制只能达到某一程度。基于这些不足，有必要进行其他反演方法的研究，如全波场反演等。
参考文献：
[1] Olhoet G R. Application of ground penetration radar[C]// Proceedings of the 6th Int’l Conf on Ground Penetration Radar. Sendadi, Japan. 1996: 1-4.
[2] 南生辉. 钻孔雷达定向测量技术与应用[J]. 煤田地质与勘察, 1998, 26(S): 56-59.
NAN Shenghui. Technique and application of directional prospecting of borehole radar system[J]. Coal Geology & Exploration, 1998, 26(S): 56-59.
[3] 黄家会, 宋雷, 崔广心. 应用跨孔雷达层析成像技术研究深部岩层特性[J]. 中国矿业大学学报, 1999, 25(6): 578-581.
HUANG Jiahui, SONG Lei, CUI Guangxin. Application of cross-hole radar tomography in studying characteristics of strata in depths[J]. Journal of China University of Mining & Technology, 1999, 25(6): 578-581.
[4] 刘四新, 佐藤源之. 时间域有限差分法(FDTD)对井中雷达的数值模拟[J]. 吉林大学学报(地球科学版), 2003, 33(4): 545-550.
LIU Sixin, Sato M. Numerical simulation for borehole radar by FDTD[J]. Journal of Jilin University (Erath Science Edition), 2003, 33(4): 545-550.
[5] 王飞, 刘四新, 吴俊军, 等. 时间域有限差分法在钻孔雷达探测金属矿中的应用[J]. 吉林大学学报(地球科学版), 2011, 40(S): 70-72.
WANG Fei, LIU Sixin, WU Junjun, et al. The application of finite difference time domain method in borehole radar for metal mine detection[J]. Journal of Jilin University (Erath Science Edition), 2011, 40(S): 70-72.
[6] 李玉喜. 钻孔雷达层析成像软件系统的研究与开发[D]. 长春: 吉林大学地球探测科学与技术学院, 2010: 36-41.
LI Yuxi. Borehole radar tomography software system study and development[D]. Changchun: University of Jilin. College of Earth Science and Technology, 2010: 36-41.
[7] 杨薇, 刘四新, 冯彦谦. 跨孔层析成像LSQR算法研究[J]. 物探与化探, 2008, 32(2): 199-202.
YANG Wei, LIU Sixin, FENG Yanqian. A study of the LSQR algorithm for cross-hole tomography[J]. Geophysical & Geochemical Exploration, 2008, 32(2): 199-202.
[8] Jens T, Dary R T, Peter D. Improved cross-hole radar tomography by using direct and reflected arrival times[J]. Journal of Applied Geophysics, 2001, 47(2): 97-105.
[9] Dale F R, Ty F. Automated water content reconstruction of zero-offset borehole ground penetrating radar data using simulated annealing[J]. Journal of Hydrology, 2005, 309(4): 1-16.
[10] Jacques R E, Alan G G, Hansruedi M. Application of a new 2D time-domain full-waveform inversion scheme to cross-hole radar data[J]. Geophysics, 2007, 72(5): J53-J64.
[11] James I, Michaael K, Rosemary K. Improving cross-hole radar velocity tomograms: A new approach to incorporating high-angle travel-time data[J]. Geophysics, 2007, 72(4): J31-J41.
[12] Bernard G, Abderrezak B, Michal C. Assisted travel-time picking of cross-hole GPR data[J]. Geophysics, 2009, 74(4): J35-J48.
[13] Gokhan G, Caglayan B. Travel-time tomography of cross-hole radar data without ray tracing[J]. Journal of Applied Geophysics, 2010, 72(4): 213-224.
[14] Evert S, Motoyuki S, Gary O. Surface and borehole ground penetrating radar developments[J]. Geophysics, 2010, 75(5): A103-A120.
[15] Knud S C, Thomas M H. Monte Carlo full-waveform inversion of cross-hole GPR data using multiple-point geo-statistical a priori information[J]. Geophysics, 2012, 77(2): H19-H31.
[16] Julien M, Agung W, Patrick B. Mapping shallow soil moisture profiles at the field scale using full-waveform inversion of ground penetrating radar data[J]. Geoderma, 2011, 167(4): 225-237.
[17] Giovanni M, Stewart G. Taming the non-linearity problem in GPR full-waveform inversion for high contrast media[J]. Journal of Applied Geophysics, 2011, 73(2): 174-186.
[18] Claudio P, Sebastien L, Mahmoudzadeh M R. Reconstruction of sub-wavelength fractures and physical properties of masonry media using full-waveform inversion of proximal penetrating radar[J]. Journal of Applied Geophysics, 2011, 74(1): 26-37.
[19] Peterson J E. Pre-inversion corrections and analysis of radar tomographic data[J]. Journal of Environmental and Engineering Geophysics, 2001, 6(1): 1-18.
[20] James I, Rosemary J K. Effect of antennas on velocity estimates obtained from cross-hole GPR data[J]. Geophysics, 2005, 70(5): K39-K42.
[21] Turner G, Siggins A F. Constant Q attenuation of subsurface radar pulses[J]. Geophysics, 1994, 59(8): 1192-1200.
[22] ZHU Ziqiang, PENG Lingxing, LU guangyin, et al. Borehole-GPR numerical simulation of full wave field based on CPML boundary[J]. Journal of Central South University, 2013, 20(3): 764-769.
(编辑 陈灿华)
收稿日期：2014-07-12；修回日期：2014-09-23
基金项目(Foundation item)：国家自然科学基金资助项目(41174061)；中南大学自由探索计划项目(2011QNZT011) ( Project(41174061) supported by the National Natural Science Foundation of China; Project(2011QNZT011) supported by Exploration Program of Central South University)
通信作者：彭凌星，博士研究生，从事地质和环境灾害探测与评价；E-mail: innocent-eyes@163.com
