基于互相關(guān)函數(shù)對鉆孔雷達層析成像的改進

朱自強，彭凌星，魯光銀，密士文

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

0 引言

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

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

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

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

1 鉆孔雷達層析成像原理

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

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

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

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

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

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

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

2 初至波的提取

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

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

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

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

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

圖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 初至波旅行時提取示意圖Fig.5 Results of our automatic picking procedure

3 層析成像實例

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

（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］陳建勝，陳從新.鉆孔雷達技術(shù)的發(fā)展和現(xiàn)狀［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，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］黃家會，宋雷，崔廣心，等.應用跨孔雷達層析成像技術(shù)研究深部巖層特性［J］.中國礦業(yè)大學學報，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ā)，徐波.鉆孔雷達探測地下裂縫［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］劉四新，周峻峰，吳俊軍.金屬礦鉆孔雷達探測的數(shù)值模擬［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］李華，焦彥杰，李富，等.鉆孔雷達層析成像技術(shù)的研究與應用［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］李玉喜.鉆孔雷達層析成像軟件系統(tǒng)的研究與開發(fā)［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］袁志亮.井間聲波電磁波層析成像技術(shù)應用研究與軟件開發(fā)［D］.北京：中國地質(zhì)大學，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.