基于互相关函数对钻孔雷达层析成像的改进

朱自强，彭凌星，鲁光银，密士文

中南大学地球科学与信息物理学院，长沙 410083

0 引言

经济的快速发展和资源的日益匮乏，使得人们对勘察技术提出了更高的要求，尤其是对探测的精度方面要求更高。钻孔雷达不仅具有地表地质雷达高分辨率的优点，而且探测深度可以达到数百米［1］，这是地表地质雷达所不具有的。因此，综合探测深度和分辨率两方面考虑，钻孔雷达具有较高的分辨率和较大的探测范围，这是探测深部地下构造较为理想的方法［2］。

钻孔雷达有3种探测方式：单孔反射探测（single－hole reflection survey）、跨孔探测（crosshole survey）、地面－孔中探测（surface－borehole survey，又称为垂直雷达剖面 VRP）［3］。单孔反射探测数据解释和地表地质雷达的数据解释方法基本相同，主要根据同相轴的连续性以及振幅的变化等来对异常位置进行判定［4］。跨孔模式中的数据解释就需要对其进行成像。目前对跨孔数据进行成像的方法大都是进行层析成像和全波场反演。2007年，Jacques等［5］采用有限差分法对跨孔雷达数据进行了全波场反演。2011年，Kim等［6］通过在频率中计算共发射源产生的电场和磁场值，并在频率域中计算出这互相正交的两者的比值，克服了全波场反演中信号源不够精确等问题，对跨孔雷达数据进行全波场反演，并得到了较好的结果。2012年，Florian等［7］在时间域中进行迭代反褶积计算发射源的子波，并据此进行全波场反演。层析成像主要分为2种，速度层析成像和衰减层析成像。考虑到天线辐射的影响，Bernard等［8］认为衰减层析成像不如旅行时层析成像的结果精确。2008年，Fernandez等［9］在MATLAB环境下，给出了速度层析成像中旅行时精确度分析的程序包。2009年，Madeleine等［10］在把孔间介质视为多层状介质的前提下，通过速度层析成像反演出每层的波速从而对异常体进行定位。2009年，Bernard等［11］提出了一种改进的旅行时提取方法，基于AIC（Alaike信息准则）和CWT（连续小波变换）的旅行时提取方法，改进速度层析成像的精度。2010年，Gokhan等［12］通过计算波前的旅行时（不是基于射线追踪）的方法来对跨孔雷达进行速度层析成像。

在国内，对钻孔雷达的研究目前大部分还只是限于应用方面：1998年，钻孔雷达在国内首次得到应用，对范各庄矿的灰岩进行了单孔反射探测［13］；2005年，刘四新等［14］通过对跨孔雷达数据进行速度层析成像探测地下含水裂缝；2010年，刘四新等［15］通过有限差分对跨孔雷达数据进行数值模拟，得到了跨孔模式下直达波和反射波的形态变化，证明了跨孔雷达在探测金属矿时同样具有较好的效果；2010年，李华等［16］通过对钻孔雷达数据进行速度和衰减层析成像，成功探测了地下岩溶的发育情况；2010年，李玉喜［17］通过对 ART、SIRT、LSQR等方法进行分析研究，对钻孔雷达的层析成像进行了较为详细的分析。

在速度层析成像中，最关键的是初至波的提取，也就是旅行时的提取。只有时间精确，才能得出更为精确的速度层析图。而在目前的研究现状中，绝大多数层析成像的旅行时提取都是采用信噪比最大法和能量比最大法。这2种方法都属于常规的方法，效果并不是很好，尤其是在收发角度较大时，很难提取出初至波的到达时间［18］。因此笔者通过互相关函数对初至波的提取做出改进。

1 钻孔雷达层析成像原理

钻孔雷达的层析成像方法有2种：一种是基于射线追踪的射线层析成像；另一种是基于波动方程的波形层析成像［19］。相对来说，前者受到的外界干扰较少，计算结果更为精确；而后者计算起来复杂许多，普通电脑的内存无法满足其计算要求。因此，实际计算中基本都采用基于射线追踪的射线层析成像，同样文中也采用的是前者。

要对跨孔雷达数据进行成像，假设雷达波在某段区域内是沿着直线传播的，那么就可以通过该段路径的长度除以初至波的到达时间来得到雷达波在此区域内的速度，从而得到所需要的速度层析成像图。图1为层析成像原理示意图，可通过图1来简单描述下旅行时层析成像的原理。

把空间介质划分为块状，分别命名为S1、S2….S25。从图1可以看出，遵循旅行时最短原理，电磁波的传播途径如图中发射天线与接收天线两者之间的连线所示。那么根据光学几何原理，可以得出

图1 层析成像原理示意图Fig.1 Sketch map of tomography principles

式中：t为初至波到达时间；s为慢速（速度的倒数）；l为电磁波传播路径。层析成像的目的就是算出孔间介质的慢度场分布。为了实现在计算机上计算慢度场的分布，假设将慢度离散化为M个细胞（即图1中的S1、S2等），那么就有M个未知的慢度值。那么旅行时就可以写为电磁波在各细胞当中传播的时间总和：

把所有道的旅行时和传播距离等写到一起，就可以得到方程：

式中：L为大量li组成的矩阵；s为大量si组成的矩阵；T为大量ti组成的矩阵。解开这个方程就可以得到人所需要的s分布图。这个方程是一个非线性的方程，根据Fermat定律，电磁波在两点之间传播的路径是沿着耗时最小的路径进行传播的，那么电磁波传播的路径将会更多地集中在高速区域。

2 初至波的提取

在钻孔雷达的实际应用中，当收发角度较大时，该道数据的信噪比将会很低，可能影响到对初至波到达时间的提取。2013年，叶佩等［20］对旅行时线性插值射线追踪的方法提出了改进，取得了较好的效果，但并没有对旅行时的提取进行改进。这里，采用一种新的方法对大收发角度数据进行初至波旅行时的提取，就是互相关函数方法。

图2为一跨孔雷达数据模拟图，此图为共发射雷达图，发射天线位于最左边地下0.25m处不动，接收天线位于最右边地下0.25m处并逐渐往下移动。当收发角度达到最大时，也即是图2中收发距离最远时，从图2明显可以看出，随着收发角度的增大，数据的信噪比越来越小。那么，在没有经过任何处理的情况下，大收发角度的数据初至波的旅行时是较难提取的，因此，提出一种改进的方法对此进行优化。具体步骤如下：

图2 发射天线位于0.25m处的跨孔雷达数据模拟图（发射点距地面深度为0.25m）Fig.2 Numerical simulation of cross－hole radar（transmitter depth＝0.25m）

1）对数据进行先期处理，对每道进行归一化处理，然后剔除废道；2）把数据按收发角度进行划分，每5°为一个单元（以收发角度0°～5°为例）；3）针对相同收发角度的数据，对每道数据进行与该角度数据中信噪比最大的一道数据进行互相关处理，这会使大部分数据按照合理的方式进行排序，如图3所示（仅以40道为例）；4）对这些数据进行叠加，从而得到收发角度0°～5°的数据参考波形，如图4所示；5）提取出参考波形的初至波到达时间；6）对每道数据与该角度的参考波形进行互相关处理，从而得到每道数据的初至波到达时间。

同样，可以得到其他角度的参考波形，与图4所示类似。这样对每个角度的参考波形与该角度数据进行互相关处理，提取出初至波到达时间。图5即为图2中数据的初至波到达时间提取效果图，粗线即为初至波的旅行时。这样，就可以得到每道数据准确的初至波旅行时，从而对数据进行层析成像。

图3 相同收发角度数据排序示意图Fig.3 Sketch map of common－ray－angle data

图4 叠加后得到的参考波形示意图Fig.4 Mean trace set by stacking common－ray－angle data

图5 初至波旅行时提取示意图Fig.5 Results of our automatic picking procedure

3 层析成像实例

3.1 合成数据的对比

为了验证笔者所提出的旅行时提取方法的准确性，通过数值模拟对其进行检验，即对合成数据进行层析成像来对该方法进行检验。因仅仅是验证，所以此次模型较为简单，模型中异常体的形状都较为规则，为矩形，异常体的磁导率和电导率与背景场相同。图6为模型的示意图。

图6 模型示意图Fig.6 Sketch map of model

图7为根据2种初至波旅行时提取方法得出的合成数据速度层析成像图。两图相比较而言，图7b更为清晰，能够分辨出模型中所有异常体的位置，并且能够对其定性。采用常规方法提取初至波旅行时得出的速度层析图能够分辨异常体的位置，但是并没有对模型中所有异常体进行定位，其中的左上和右下2个位置的异常体并没有在图中得到反映。

图7 合成数据的速度层析成像图Fig.7 Synthetic data results

3.2 实际数据的对比

为了验证笔者提出的方法在针对实际数据的处理是否有效，笔者在湖南某高速公路的路基探溶中做了实验。该路段受地形、地层岩性及地下水等因素影响，沿线不良地质现象及特殊性岩土主要为岩溶、软土和红黏土。

图8 钻孔雷达实际探测数据图（发射点距地面深度为7m）Fig.8 Field data of borehole radar（transmitter depth＝7m）

此次数据采集应用的是瑞典MALA公司的RAMAC地质雷达。图8为其中的一组数据（为便于观看，将图旋转了90°），发射天线位于距离地面7 m处的钻孔内固定不动，接收天线从另一钻孔内匀速往下运动，采样间隔为0.2m。那么收发角度就是从大到小，到接收天线也位于距离地面7m处收发角度为0°，然后又逐渐增大。所以从图8中可以看出，随着收发角度的增大，信噪比在降低，初至波提取也变得越来越难。

采集数据的钻孔分布为：ZK1和ZK2分别为雷达探测孔和接收孔，两孔相隔为4m；ZK3为验证孔，位于ZK1和ZK2之间，与ZK2的距离约为1 m。图9为用2种方法对实际采集的跨孔数据进行的速度层析成像图。

图9 实际采集数据的速度层析成像图Fig.9 Velocity tomography of field data

通过对比图9a和9b可以看出，互相关处理旅行时得出的层析成像图在效果上比普通方法得出的结果好。结合最后对探测结果进行的钻孔验证可知：岩心显示9～10m处为空洞，说明此处电磁波速度应该为高速异常；采用互相关处理旅行时得出的图9b压制了图9a中深度为4m处的假异常，并且突出了深度为9～10m处的高速异常体的特征。由此可以得知，通过互相关函数对旅行时进行处理确实能够提高速度层析成像的准确度。因此，对于钻孔雷达探测隐伏岩溶，可以更准确地探测出隐伏岩溶体的特征和规模。

4 结论

1）在针对大收发角度数据进行初至波提取时，因其具有很低的信噪比，从而较难提取出准确的旅行时。笔者对同收发角度的数据进行排序、叠加得到参考波形，并对参考波形和该收发角度的数据进行互相关处理，从而得到每道数据的初至波旅行时。根据对合成数据的处理来看，该方法对大收发角度数据的初至波提取具有很好的效果。

2）通过对不同初至波提取方法得出的层析成像进行对比可以看出，采用互相关函数提取初至波旅行时的方法得到的层析成像图更为精确，对异常的分辨也更为准确。这也证明互相关函数在改进初至波旅行时提取中的有效性，为钻孔雷达在对隐伏岩溶的探测提供了更大的帮助。

3）在基于射线的射线层析成像当中，计算区域的4个角落处没有足够的射线覆盖，从而使得该处的层析成像图没有计算区域中部的层析图那样精确。因此，下一步的研究工作是联合反射波和直达波进行层析成像，这样就能保证射线可以覆盖到更多的区域，从而使得层析成像图更为精确，这样对隐伏岩溶的判断更为准确和直接。

（References）：

［1］Gilles B，Michel C.Massive Sulphide Delineation U－sing Borehole Radar：Tests at the McConnel Nickel Deposit［J］.Journal of Applied Geophysics，2001，47（1）：45－61.

［2］陈建胜，陈从新.钻孔雷达技术的发展和现状［J］.地球物理学进展，2008，23（5）：1634－1640.Chen Jiansheng，Chen Congxin.The Review of Borehole Radar Technology［J］.Progress in Geophysics，2008，23（5）：1634－1640.

［3］彭凌星.钻孔雷达正演模拟及应用研究［D］.长沙：中南大学，2010.Peng Lingxing.Study on Simulation and Application of Borehole－GPR［D］.Changsha：Central South Uni－versity，2010.

［4］Harry M J.Ground Penetrating Radar：Theory and Applications［M］.Oxford：Elsevier Science ＆ Technology Rights Department，2009.

［5］Jacques R E，Alan G G，Hansruedi M，et al.Application of a New 2DTime－Domain Full－Waveform Inversion Scheme to Cross－Hole Radar Data［J］.Geophysics，2007，72（5）：53－64.

［6］Kim J H，Kobayashi T，Lee S K.Admittance Inversion of GPR Transmission for Cross－Hole Tomography［J］.Journal of Applied Geophysics，2012，87（2）：57－67.

［7］Florian B，James I，Jacques E，et al.Analysis of an Iterative Deconvolution Approach for Estimating the Source Wavelet During Waveform Inversion of Cross－Hole Geo－Radar Data［J］.Journal of Applied Geophysics，2012，78（3）：20－30.

［8］Bernard G，Erwan G，Michel C.Bh＿Tomo：A Matlab Borehole Geo－Radar 2DTomography Package［J］.Computers ＆ Geosciences，2007，33（1）：126－137.

［9］Fernández－Martínez J L，Fernández－Alvarez J P，Pedruelo－González L M.MTCLAB：A Matlab－Based Program for Travel－Time Quality Analysis and Pre－Inversion Velocity Tuning in 2DTransmission Tomography［J］.Computers ＆ Geosciences，2008，34（3）：213－225.

［10］Madeleine M，Perroud H，Dominique R.First－Arrival Travel－Times Inversion Based on a Minimal Number of Parameters in Shallow Cross－Well GPR Tomography［J］.Journal of Applied Geophysics，2009，67（4）：278－287.

［11］Bernard G，Abderrezak B，Michel C.Assisted Travel－Time Picking of Cross－Hole GPR Data［J］.Geophysics，2009，72（4）：35－48.

［12］Gökhan G，Çaˇglayan B.Travel－Time Tomography of Cross－Hole Radar Data Without Ray Tracing［J］.Journal of Applied Geophysics，2010，72（4）：213－224.

［13］黄家会，宋雷，崔广心，等.应用跨孔雷达层析成像技术研究深部岩层特性［J］.中国矿业大学学报，1999，25（6）：578－581.Huang Jiahui，Song Lei，Cui Guangxin，et al.Application of Crosshole Radar Tomography in Studying Characteristics of Strata in Depths［J］.Journal of China Mining ＆Technology，1999，25（6）：578－581.

［14］刘四新，曾昭发，徐波.钻孔雷达探测地下裂缝［J］.吉林大学学报：地球科学版，2005，35（增刊）：128－131.Liu Sixin，Zeng Zhaofa，Xu Bo.Subsurface Water－Filled Fractures Detection by Borehole Radar［J］.Journal of Jilin University：Earth Science Edition，2005，35（Sup.）：128－131.

［15］刘四新，周峻峰，吴俊军.金属矿钻孔雷达探测的数值模拟［J］.吉林大学学报：地球科学版，2010，40（6）：1479－1484.Liu Sixin，Zhou Junfeng，Wu Junjun.Numerical Simulation of Borehole Radar for Metal Ore Detection［J］.Journal of Jilin University：Earth Science Edition，2010，40（6）：1479－1484.

［16］李华，焦彦杰，李富，等.钻孔雷达层析成像技术的研究与应用［J］.地球物理学进展，2010，25（5）：1863－1869.Li Hua，Jiao Yanjie，Li Fu，et al.The Study and Application of Borehole Radar Tomography Technology［J］.Progress in Geophysics，2010，25（5）：1863－1869.

［17］李玉喜.钻孔雷达层析成像软件系统的研究与开发［D］.长春：吉林大学，2010.Li Yuxi.Borehole Radar Tomography Software System Study and Development［D］.Changchun：Jilin University，2010.

［18］James D I，Michael D K，Rosemary J K.Improving Cross－Hole Radar Velocity Tomograms：A New Approach to Incorporating High－Angle Travel－Time Data［J］.2007，72（4）：31－41.

［19］袁志亮.井间声波电磁波层析成像技术应用研究与软件开发［D］.北京：中国地质大学，2007.Yuan Zhiliang.Cross－Well Acoustic Wave and Electromagnetic Wave Computerized Tomography Application and Software Opened［D］.Beijing：China University of Geosciences，2007.

［20］叶佩，李庆春.旅行时线性插值射线追踪提高计算精度和效率的改进方法［J］.吉林大学学报：地球科学版，2013，43（1）：291－298.Ye Pei，Li Qingchun.Improvements of Linear Traveltime Interpolation Ray Tracing for the Accuracy and Efficiency［J］.Journal of Jilin University：Earth Science Edition，2013，43（1）：291－298.