朱自強(qiáng)，彭凌星，魯光銀，密士文

中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083

0 引言

經(jīng)濟(jì)的快速發(fā)展和資源的日益匱乏，使得人們對勘察技術(shù)提出了更高的要求，尤其是對探測的精度方面要求更高。鉆孔雷達(dá)不僅具有地表地質(zhì)雷達(dá)高分辨率的優(yōu)點(diǎn)，而且探測深度可以達(dá)到數(shù)百米［1］，這是地表地質(zhì)雷達(dá)所不具有的。因此，綜合探測深度和分辨率兩方面考慮，鉆孔雷達(dá)具有較高的分辨率和較大的探測范圍，這是探測深部地下構(gòu)造較為理想的方法［2］。

鉆孔雷達(dá)有3種探測方式：單孔反射探測（single－h(huán)ole reflection survey）、跨孔探測（crosshole survey）、地面－孔中探測（surface－borehole survey，又稱為垂直雷達(dá)剖面 VRP）［3］。單孔反射探測數(shù)據(jù)解釋和地表地質(zhì)雷達(dá)的數(shù)據(jù)解釋方法基本相同，主要根據(jù)同相軸的連續(xù)性以及振幅的變化等來對異常位置進(jìn)行判定［4］?？缈啄Ｊ街械臄?shù)據(jù)解釋就需要對其進(jìn)行成像。目前對跨孔數(shù)據(jù)進(jìn)行成像的方法大都是進(jìn)行層析成像和全波場反演。2007年，Jacques等［5］采用有限差分法對跨孔雷達(dá)數(shù)據(jù)進(jìn)行了全波場反演。2011年，Kim等［6］通過在頻率中計(jì)算共發(fā)射源產(chǎn)生的電場和磁場值，并在頻率域中計(jì)算出這互相正交的兩者的比值，克服了全波場反演中信號源不夠精確等問題，對跨孔雷達(dá)數(shù)據(jù)進(jìn)行全波場反演，并得到了較好的結(jié)果。2012年，F(xiàn)lorian等［7］在時(shí)間域中進(jìn)行迭代反褶積計(jì)算發(fā)射源的子波，并據(jù)此進(jìn)行全波場反演。層析成像主要分為2種，速度層析成像和衰減層析成像?？紤]到天線輻射的影響，Bernard等［8］認(rèn)為衰減層析成像不如旅行時(shí)層析成像的結(jié)果精確。2008年，F(xiàn)ernandez等［9］在MATLAB環(huán)境下，給出了速度層析成像中旅行時(shí)精確度分析的程序包。2009年，Madeleine等［10］在把孔間介質(zhì)視為多層狀介質(zhì)的前提下，通過速度層析成像反演出每層的波速從而對異常體進(jìn)行定位。2009年，Bernard等［11］提出了一種改進(jìn)的旅行時(shí)提取方法，基于AIC（Alaike信息準(zhǔn)則）和CWT（連續(xù)小波變換）的旅行時(shí)提取方法，改進(jìn)速度層析成像的精度。2010年，Gokhan等［12］通過計(jì)算波前的旅行時(shí)（不是基于射線追蹤）的方法來對跨孔雷達(dá)進(jìn)行速度層析成像。

在國內(nèi)，對鉆孔雷達(dá)的研究目前大部分還只是限于應(yīng)用方面：1998年，鉆孔雷達(dá)在國內(nèi)首次得到應(yīng)用，對范各莊礦的灰?guī)r進(jìn)行了單孔反射探測［13］；2005年，劉四新等［14］通過對跨孔雷達(dá)數(shù)據(jù)進(jìn)行速度層析成像探測地下含水裂縫；2010年，劉四新等［15］通過有限差分對跨孔雷達(dá)數(shù)據(jù)進(jìn)行數(shù)值模擬，得到了跨孔模式下直達(dá)波和反射波的形態(tài)變化，證明了跨孔雷達(dá)在探測金屬礦時(shí)同樣具有較好的效果；2010年，李華等［16］通過對鉆孔雷達(dá)數(shù)據(jù)進(jìn)行速度和衰減層析成像，成功探測了地下巖溶的發(fā)育情況；2010年，李玉喜［17］通過對 ART、SIRT、LSQR等方法進(jìn)行分析研究，對鉆孔雷達(dá)的層析成像進(jìn)行了較為詳細(xì)的分析。

在速度層析成像中，最關(guān)鍵的是初至波的提取，也就是旅行時(shí)的提取。只有時(shí)間精確，才能得出更為精確的速度層析圖。而在目前的研究現(xiàn)狀中，絕大多數(shù)層析成像的旅行時(shí)提取都是采用信噪比最大法和能量比最大法。這2種方法都屬于常規(guī)的方法，效果并不是很好，尤其是在收發(fā)角度較大時(shí)，很難提取出初至波的到達(dá)時(shí)間［18］。因此筆者通過互相關(guān)函數(shù)對初至波的提取做出改進(jìn)。

1 鉆孔雷達(dá)層析成像原理

鉆孔雷達(dá)的層析成像方法有2種：一種是基于射線追蹤的射線層析成像；另一種是基于波動方程的波形層析成像［19］。相對來說，前者受到的外界干擾較少，計(jì)算結(jié)果更為精確；而后者計(jì)算起來復(fù)雜許多，普通電腦的內(nèi)存無法滿足其計(jì)算要求。因此，實(shí)際計(jì)算中基本都采用基于射線追蹤的射線層析成像，同樣文中也采用的是前者。

要對跨孔雷達(dá)數(shù)據(jù)進(jìn)行成像，假設(shè)雷達(dá)波在某段區(qū)域內(nèi)是沿著直線傳播的，那么就可以通過該段路徑的長度除以初至波的到達(dá)時(shí)間來得到雷達(dá)波在此區(qū)域內(nèi)的速度，從而得到所需要的速度層析成像圖。圖1為層析成像原理示意圖，可通過圖1來簡單描述下旅行時(shí)層析成像的原理。

把空間介質(zhì)劃分為塊狀，分別命名為S1、S2….S25。從圖1可以看出，遵循旅行時(shí)最短原理，電磁波的傳播途徑如圖中發(fā)射天線與接收天線兩者之間的連線所示。那么根據(jù)光學(xué)幾何原理，可以得出

圖1 層析成像原理示意圖Fig.1 Sketch map of tomography principles

式中：t為初至波到達(dá)時(shí)間；s為慢速（速度的倒數(shù)）；l為電磁波傳播路徑。層析成像的目的就是算出孔間介質(zhì)的慢度場分布。為了實(shí)現(xiàn)在計(jì)算機(jī)上計(jì)算慢度場的分布，假設(shè)將慢度離散化為M個(gè)細(xì)胞（即圖1中的S1、S2等），那么就有M個(gè)未知的慢度值。那么旅行時(shí)就可以寫為電磁波在各細(xì)胞當(dāng)中傳播的時(shí)間總和：

把所有道的旅行時(shí)和傳播距離等寫到一起，就可以得到方程：

式中：L為大量li組成的矩陣；s為大量si組成的矩陣；T為大量ti組成的矩陣。解開這個(gè)方程就可以得到人所需要的s分布圖。這個(gè)方程是一個(gè)非線性的方程，根據(jù)Fermat定律，電磁波在兩點(diǎn)之間傳播的路徑是沿著耗時(shí)最小的路徑進(jìn)行傳播的，那么電磁波傳播的路徑將會更多地集中在高速區(qū)域。

2 初至波的提取

在鉆孔雷達(dá)的實(shí)際應(yīng)用中，當(dāng)收發(fā)角度較大時(shí)，該道數(shù)據(jù)的信噪比將會很低，可能影響到對初至波到達(dá)時(shí)間的提取。2013年，葉佩等［20］對旅行時(shí)線性插值射線追蹤的方法提出了改進(jìn)，取得了較好的效果，但并沒有對旅行時(shí)的提取進(jìn)行改進(jìn)。這里，采用一種新的方法對大收發(fā)角度數(shù)據(jù)進(jìn)行初至波旅行時(shí)的提取，就是互相關(guān)函數(shù)方法。

圖2為一跨孔雷達(dá)數(shù)據(jù)模擬圖，此圖為共發(fā)射雷達(dá)圖，發(fā)射天線位于最左邊地下0.25m處不動，接收天線位于最右邊地下0.25m處并逐漸往下移動。當(dāng)收發(fā)角度達(dá)到最大時(shí)，也即是圖2中收發(fā)距離最遠(yuǎn)時(shí)，從圖2明顯可以看出，隨著收發(fā)角度的增大，數(shù)據(jù)的信噪比越來越小。那么，在沒有經(jīng)過任何處理的情況下，大收發(fā)角度的數(shù)據(jù)初至波的旅行時(shí)是較難提取的，因此，提出一種改進(jìn)的方法對此進(jìn)行優(yōu)化。具體步驟如下：

圖2 發(fā)射天線位于0.25m處的跨孔雷達(dá)數(shù)據(jù)模擬圖（發(fā)射點(diǎn)距地面深度為0.25m）Fig.2 Numerical simulation of cross－h(huán)ole radar（transmitter depth＝0.25m）

1）對數(shù)據(jù)進(jìn)行先期處理，對每道進(jìn)行歸一化處理，然后剔除廢道；2）把數(shù)據(jù)按收發(fā)角度進(jìn)行劃分，每5°為一個(gè)單元（以收發(fā)角度0°～5°為例）；3）針對相同收發(fā)角度的數(shù)據(jù)，對每道數(shù)據(jù)進(jìn)行與該角度數(shù)據(jù)中信噪比最大的一道數(shù)據(jù)進(jìn)行互相關(guān)處理，這會使大部分?jǐn)?shù)據(jù)按照合理的方式進(jìn)行排序，如圖3所示（僅以40道為例）；4）對這些數(shù)據(jù)進(jìn)行疊加，從而得到收發(fā)角度0°～5°的數(shù)據(jù)參考波形，如圖4所示；5）提取出參考波形的初至波到達(dá)時(shí)間；6）對每道數(shù)據(jù)與該角度的參考波形進(jìn)行互相關(guān)處理，從而得到每道數(shù)據(jù)的初至波到達(dá)時(shí)間。

同樣，可以得到其他角度的參考波形，與圖4所示類似。這樣對每個(gè)角度的參考波形與該角度數(shù)據(jù)進(jìn)行互相關(guān)處理，提取出初至波到達(dá)時(shí)間。圖5即為圖2中數(shù)據(jù)的初至波到達(dá)時(shí)間提取效果圖，粗線即為初至波的旅行時(shí)。這樣，就可以得到每道數(shù)據(jù)準(zhǔn)確的初至波旅行時(shí)，從而對數(shù)據(jù)進(jìn)行層析成像。

圖3 相同收發(fā)角度數(shù)據(jù)排序示意圖Fig.3 Sketch map of common－ray－angle data

圖4 疊加后得到的參考波形示意圖Fig.4 Mean trace set by stacking common－ray－angle data

圖5 初至波旅行時(shí)提取示意圖Fig.5 Results of our automatic picking procedure

3 層析成像實(shí)例

3.1 合成數(shù)據(jù)的對比

為了驗(yàn)證筆者所提出的旅行時(shí)提取方法的準(zhǔn)確性，通過數(shù)值模擬對其進(jìn)行檢驗(yàn)，即對合成數(shù)據(jù)進(jìn)行層析成像來對該方法進(jìn)行檢驗(yàn)。因僅僅是驗(yàn)證，所以此次模型較為簡單，模型中異常體的形狀都較為規(guī)則，為矩形，異常體的磁導(dǎo)率和電導(dǎo)率與背景場相同。圖6為模型的示意圖。

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

圖7為根據(jù)2種初至波旅行時(shí)提取方法得出的合成數(shù)據(jù)速度層析成像圖。兩圖相比較而言，圖7b更為清晰，能夠分辨出模型中所有異常體的位置，并且能夠?qū)ζ涠ㄐ?。采用常?guī)方法提取初至波旅行時(shí)得出的速度層析圖能夠分辨異常體的位置，但是并沒有對模型中所有異常體進(jìn)行定位，其中的左上和右下2個(gè)位置的異常體并沒有在圖中得到反映。

圖7 合成數(shù)據(jù)的速度層析成像圖Fig.7 Synthetic data results

3.2 實(shí)際數(shù)據(jù)的對比

為了驗(yàn)證筆者提出的方法在針對實(shí)際數(shù)據(jù)的處理是否有效，筆者在湖南某高速公路的路基探溶中做了實(shí)驗(yàn)。該路段受地形、地層巖性及地下水等因素影響，沿線不良地質(zhì)現(xiàn)象及特殊性巖土主要為巖溶、軟土和紅黏土。

圖8 鉆孔雷達(dá)實(shí)際探測數(shù)據(jù)圖（發(fā)射點(diǎn)距地面深度為7m）Fig.8 Field data of borehole radar（transmitter depth＝7m）

此次數(shù)據(jù)采集應(yīng)用的是瑞典MALA公司的RAMAC地質(zhì)雷達(dá)。圖8為其中的一組數(shù)據(jù)（為便于觀看，將圖旋轉(zhuǎn)了90°），發(fā)射天線位于距離地面7 m處的鉆孔內(nèi)固定不動，接收天線從另一鉆孔內(nèi)勻速往下運(yùn)動，采樣間隔為0.2m。那么收發(fā)角度就是從大到小，到接收天線也位于距離地面7m處收發(fā)角度為0°，然后又逐漸增大。所以從圖8中可以看出，隨著收發(fā)角度的增大，信噪比在降低，初至波提取也變得越來越難。

采集數(shù)據(jù)的鉆孔分布為：ZK1和ZK2分別為雷達(dá)探測孔和接收孔，兩孔相隔為4m；ZK3為驗(yàn)證孔，位于ZK1和ZK2之間，與ZK2的距離約為1 m。圖9為用2種方法對實(shí)際采集的跨孔數(shù)據(jù)進(jìn)行的速度層析成像圖。

圖9 實(shí)際采集數(shù)據(jù)的速度層析成像圖Fig.9 Velocity tomography of field data

通過對比圖9a和9b可以看出，互相關(guān)處理旅行時(shí)得出的層析成像圖在效果上比普通方法得出的結(jié)果好。結(jié)合最后對探測結(jié)果進(jìn)行的鉆孔驗(yàn)證可知：巖心顯示9～10m處為空洞，說明此處電磁波速度應(yīng)該為高速異常；采用互相關(guān)處理旅行時(shí)得出的圖9b壓制了圖9a中深度為4m處的假異常，并且突出了深度為9～10m處的高速異常體的特征。由此可以得知，通過互相關(guān)函數(shù)對旅行時(shí)進(jìn)行處理確實(shí)能夠提高速度層析成像的準(zhǔn)確度。因此，對于鉆孔雷達(dá)探測隱伏巖溶，可以更準(zhǔn)確地探測出隱伏巖溶體的特征和規(guī)模。

4 結(jié)論

1）在針對大收發(fā)角度數(shù)據(jù)進(jìn)行初至波提取時(shí)，因其具有很低的信噪比，從而較難提取出準(zhǔn)確的旅行時(shí)。筆者對同收發(fā)角度的數(shù)據(jù)進(jìn)行排序、疊加得到參考波形，并對參考波形和該收發(fā)角度的數(shù)據(jù)進(jìn)行互相關(guān)處理，從而得到每道數(shù)據(jù)的初至波旅行時(shí)。根據(jù)對合成數(shù)據(jù)的處理來看，該方法對大收發(fā)角度數(shù)據(jù)的初至波提取具有很好的效果。

2）通過對不同初至波提取方法得出的層析成像進(jìn)行對比可以看出，采用互相關(guān)函數(shù)提取初至波旅行時(shí)的方法得到的層析成像圖更為精確，對異常的分辨也更為準(zhǔn)確。這也證明互相關(guān)函數(shù)在改進(jìn)初至波旅行時(shí)提取中的有效性，為鉆孔雷達(dá)在對隱伏巖溶的探測提供了更大的幫助。

3）在基于射線的射線層析成像當(dāng)中，計(jì)算區(qū)域的4個(gè)角落處沒有足夠的射線覆蓋，從而使得該處的層析成像圖沒有計(jì)算區(qū)域中部的層析圖那樣精確。因此，下一步的研究工作是聯(lián)合反射波和直達(dá)波進(jìn)行層析成像，這樣就能保證射線可以覆蓋到更多的區(qū)域，從而使得層析成像圖更為精確，這樣對隱伏巖溶的判斷更為準(zhǔn)確和直接。

（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］陳建勝，陳從新.鉆孔雷達(dá)技術(shù)的發(fā)展和現(xiàn)狀［J］.地球物理學(xué)進(jìn)展，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á)正演模擬及應(yīng)用研究［D］.長沙：中南大學(xué)，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，F(xiàn)erná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］黃家會，宋雷，崔廣心，等.應(yīng)用跨孔雷達(dá)層析成像技術(shù)研究深部巖層特性［J］.中國礦業(yè)大學(xué)學(xué)報(bào)，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］劉四新，曾昭發(fā)，徐波.鉆孔雷達(dá)探測地下裂縫［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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］劉四新，周峻峰，吳俊軍.金屬礦鉆孔雷達(dá)探測的數(shù)值模擬［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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］李華，焦彥杰，李富，等.鉆孔雷達(dá)層析成像技術(shù)的研究與應(yīng)用［J］.地球物理學(xué)進(jìn)展，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á)層析成像軟件系統(tǒng)的研究與開發(fā)［D］.長春：吉林大學(xué)，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］袁志亮.井間聲波電磁波層析成像技術(shù)應(yīng)用研究與軟件開發(fā)［D］.北京：中國地質(zhì)大學(xué)，2007.Yuan Zhiliang.Cross－Well Acoustic Wave and Electromagnetic Wave Computerized Tomography Application and Software Opened［D］.Beijing：China University of Geosciences，2007.

［20］葉佩，李慶春.旅行時(shí)線性插值射線追蹤提高計(jì)算精度和效率的改進(jìn)方法［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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.