朱自強，彭凌星，密士文

?

基于旅行時修正的鉆孔雷達層析成像改進

朱自強，彭凌星，密士文

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

基于鉆孔雷達的波幅采用速度層析成像時，大收發(fā)角度雷達波幅因其信噪比低，旅行時提取較困難等問題，鉆孔雷達雷達波幅層析成像精度不高。通過計算證明電磁波能量是從發(fā)射天線中心點傳播至末端，再由接收天線末端傳播至其中心。依據此傳播路徑采用互相關函數(shù)對旅行時的提取進行優(yōu)化，并得到相應的電磁波波速。利用優(yōu)化后的旅行時與波速得到旅行時的修正值，并根據修正后的旅行時進行速度層析成像，將進行旅行時修正的層析成像圖與未進行旅行時修正的層析成像圖進行對比。研究結果表明：經過旅行時修正后得到的層析成像結果更加精確，證明對鉆孔雷達層析成像的改進方法是可行的。

鉆孔雷達；初至波旅行時；互相關函數(shù)；旅行時修正；層析成像

20世紀70年代之后，在地質雷達的基礎上，Olhoet[1]為了充分利用鉆孔這一已有的通道對地下介質的結構特征等信息進行探測，提出一種新的物探方法即鉆孔雷達探測。鉆孔雷達在國內的發(fā)展起步較晚，到20世紀末國內才出現(xiàn)有關鉆孔雷達的研究，但基本上以介紹為主，之后漸漸對其應用和理論進行研究[2]。黃家會等[3]應用鉆孔雷達對地下深部灰?guī)r地區(qū)的裂隙帶和水流通道等進行探測，通過對雷達波幅進行成像對高衰減低速區(qū)、低衰減高速區(qū)與低衰減低速區(qū)等進行區(qū)分。劉四新等[4]通過有限差分數(shù)值模擬了充水裂縫和礦體在鉆孔雷達中的響應特征，證明鉆孔雷達在裂縫探測中的可行性。王飛等[5]通過有限差分模擬證明了鉆孔雷達單孔反射探測可以有效確定礦體的位置和形態(tài)。在鉆孔雷達的成像研究方面，李玉喜[6]采用GPR-MAX鉆孔雷達進行正演模擬，并通過反演實例討論了射線覆蓋角度、正則化因子等對成像的影響。楊薇等[7]針對解答大型稀疏矩陣的LSQR方法進行了研究，利用鉆孔雷達層析成像直觀反映了地下介質的速度、吸收系數(shù)等參數(shù)。在鉆孔雷達的成像理論研究方面，Jens等[8?9]結合鉆孔雷達直達波和反射波、折射波進行聯(lián)合成像，提高了成像的精度，能夠有效得到分辨地下介質結構。Jacques等[10]在基于時域有限差分法的全波場反演基礎上，對2組跨孔雷達實測雷達波幅進行反演，反演結果與實際地質結果基本一致。James等[11]通過對旅行時提取的改進來聯(lián)合所有雷達波幅進行層析成像，有效降低了層析成像圖中大收發(fā)角度雷達波幅引起的假象。Bernard等[12]基于Akaike信息準則(AIC)和連續(xù)小波變換(CWT)改進了層析成像中旅行時的提取方法，Gokhan等[13]通過程函方程計算雷達波旅行時。還有一些學者在鉆孔雷達成像方面采用全波場反演[14?18]。在這些研究中，主要采用層析成像方法進行鉆孔雷達成像，這是應用最廣泛的方法之一。本文也采用層析成像對鉆孔雷達雷達波幅進行處理。針對大收發(fā)角度雷達波幅與小角度雷達波幅之間的差異，利用電磁波在天線中傳播速度和介質中傳播速度之差對旅行時進行修正補償，從而提高旅行時的準確度以便提高層析成像的效果。

1 旅行時提取優(yōu)化

鉆孔雷達的跨孔探測方式如圖1所示。探測方式有ZOP(零偏移距剖面)和MOG(多偏移距雷達探測)共2種，其中MOG方式較常用，對其雷達波幅進行處理的方式也是主要通過層析成像來進行。當所有接收點都接收雷達波幅時，大角度收發(fā)雷達波幅因其傳播距離較長，會導致電磁波衰減更加強烈。

圖1 跨孔雷達探測示意圖

取1個鉆孔雷達實測雷達波幅來對此進行說明。使用瑞典MALA公司的RAMAC地質雷達采集雷達波幅，其天線的主頻率為250 MHz。發(fā)射天線位于地下11 m處固定不動，接收天線沿接收鉆孔從井口勻速下移17 m。對實驗雷達波幅進行濾波、歸一化振幅以及頻譜分析等處理得到的結果如圖2所示。隨著接收天線深度逐漸增大，收發(fā)角度由極大值減小至0°，再逐漸增大。從圖2可以看出：隨著收發(fā)角度的增大，信號的主頻也增大，信噪比則慢慢減小。此時，初至波的提取就會變得越來越困難。

(a) 深度與時間的關系；(b) 深度與頻率的關系；(c) 信噪比與深度的關系

采用一種新的方法即互相關函數(shù)方法對大收發(fā)角度雷達波幅提取初至波旅行時。其步驟為：

1)對雷達波幅進行先期處理，對每道進行歸一化處理，然后剔除廢道。

2) 將雷達波幅按收發(fā)角度進行劃分，每5°為1個單元(以收發(fā)角度0°~5°為例)。

3) 針對相同收發(fā)角度的雷達波幅，對每道雷達波幅與該角度雷達波幅中信噪比最大的一道雷達波幅進行互相關處理，這會使大部分雷達波幅按照合理的方式進行排序，如圖3所示。

圖3 相同收發(fā)角度雷達波幅排序

4) 對這些雷達波幅進行疊加，從而得到收發(fā)角為0°~5°時的雷達波幅參考波形，如圖4所示。

圖4 相同收發(fā)角度雷達波幅疊加得到的參考波形

5) 提取出參考波形的初至波到達時間。

6) 對每道雷達波幅與該角度的參考波形進行互相關處理，從而得到每道雷達波幅的初至波到達時間。

同樣，可以得到其他角度的參考波形，與圖4類似。這樣對每個角度的參考波形與該角度雷達波幅進行互相關處理，提取出初至波到達時間。

圖5所示為圖2中不同發(fā)射點的雷達波幅初至波到達時間提取效果圖(取4個雷達波幅為例)，其中黑色粗線對應初至波的旅行時。從圖5可以看出：經過互相關處理后的旅行時提取更加精確。這樣，就可以得到每道雷達波幅的準確初至波旅行時，從而對雷達波幅進行層析成像。

發(fā)射天線距離/m：(a) 0.25；(b) 2.50；(c) 5.00；(d) 10.00

2 速度的計算分析

即使可以較準確地提取出大收發(fā)角度雷達波幅的初至波旅行時，大收發(fā)角度雷達波幅和小收發(fā)角度雷達波幅仍存在一定的差異。圖6所示為一數(shù)值模擬結果，孔間介質的介電常數(shù)隨機分布在20~30之間，取發(fā)射點位于孔中間時接收的雷達波幅成圖，這樣可以計算出每道雷達波幅的平均波速，如圖7所示。

圖6 跨孔雷達探測數(shù)值模擬結果

圖7 傳播路徑平均波速

從圖7可以看出：電磁波平均速度隨著接收角度的增大而增大然后慢慢減小。大收發(fā)角度的雷達波幅明顯與小收發(fā)角度的雷達波幅波速相差較大，因此，將所有角度的雷達波幅直接進行聯(lián)合成像會造成成像效果的誤差。Peterson[19]闡述了這種類似現(xiàn)象，但是他并沒有闡述發(fā)生這種現(xiàn)象的原因。James等[20]在2005年通過對數(shù)值模擬，證明了不管鉆孔內是空氣充填還是水充填，電磁波在天線中傳播速度通常要比在介質中傳播速度快。在此前提下，基于電磁波傳播不是在天線中心之間傳播而是通過天線和孔間介質傳播的假設，本文對旅行時進行計算和分析。

在傳統(tǒng)的基于射線理論的速度層析成像中，都是將發(fā)射天線和接收天線看成1個點，電磁波沿著2個點之間傳播。而在鉆孔雷達探測的實際過程中，發(fā)射天線和接收天線都是有一定長度的，并且電磁波在天線上傳播的速度比在空間介質中傳播的速度快，因此，在這個基礎上，假設電磁波傳播的路徑如圖8所示。圖8中帶箭頭的線為電磁波傳播路徑，電磁波首先沿著接收天線傳播至端頭，然后經由孔間介質傳播至接收天線端頭，最后到達接收天線中心。

圖8 射線傳播路徑示意圖

按照此理論，可以計算出每道雷達波幅平均波速的理論值。因此，設置1個模型，并據此計算其電磁波的波速。仍然以發(fā)射點位于中心位置時為例?？组g介質為各向均勻介質，介電常數(shù)為25；電磁波在天線中傳播的速度為0.11 m/ns；天線長度為0.8 m；2個孔之間的距離為5 m。圖9所示為計算的每道雷達波幅的平均波速。從圖9可以看出：波速隨著角度的增大會逐漸增大然后慢慢減小。但在收發(fā)角度較小時，速度基本不變。

圖9 波速理論值示意圖

但這種曲線的形態(tài)過于理想化，是建立在脈沖持續(xù)時間相對于電磁波在天線中傳播的時間要小得多的基礎上。實際上，對于現(xiàn)在所使用的儀器來說不大可能實現(xiàn)。因此，計算時，可以采用10 ns高斯脈沖作為代表進行計算。在不考慮介質的色散情況下，所得計算結果如圖10所示。

圖10 高斯脈沖波速理論值

考慮到實際的鉆孔雷達探測過程中，介質會使電磁波發(fā)生一定程度色散，這里采用與圖9所示模型中一樣的參數(shù)和脈沖，給予電磁波一定色散，分析色散對不同收發(fā)角度雷達波幅波速的影響。Turner等[21]給出色散系數(shù)的取值范圍為2~30。本文取為16，所得結果如圖11所示。

圖11 色散下的高斯脈沖波速理論值

從圖11可以看出：在基于色散和長脈沖的基礎上電磁波波速計算結果與根據模擬雷達波幅計算出的波速分布圖非常接近，曲線形狀也基本一樣。由此可以推斷：在鉆孔雷達實際探測過程中電磁波的傳播路徑是按照圖8所示的路徑進行傳播的。層析成像基于射線從天線中心點之間按直線傳播的理論進行計算，所以，必須根據實際傳播路徑對其旅行時進行修正，從而更接近于實際情況。

3 旅行時修正

按照上述計算結果，在對鉆孔雷達雷達波幅進行層析成像時，需要對其旅行時進行修正。標準的基于射線理論的速度層析成像方程為

其中：為傳播路徑的長度；為電磁波在路徑中傳播的慢度(速度的倒數(shù))；為旅行時。根據前面的分析，需要對旅行時進行一部分修正?？梢詫⒎匠谈膶憺?/p>

(2)

在解方程(3)時，因為電磁波的波速較大，所以，旅行時的修正值會很小。故在對修正值曲線進行校正時，不需要對其進行平滑和對稱約束處理。需要注意的是：當收發(fā)角度為0°(收發(fā)天線位于同一水平位置)時，該角度的也為0°。通過解式(3)，可以將所有收發(fā)角度的雷達波幅進行速度層析成像。

為了對上述方法進行驗證，通過數(shù)值模擬對模型層析成像進行對比，圖12所示為模型示意圖。圖12中，異常體形狀都較簡單，為規(guī)則的矩形；背景是介電常數(shù)為25的各向均勻介質，左上和右下2個異常體的介電常數(shù)為29，其余異常體的介電常數(shù)為21；所有介質的電導率均設為1 mS/m；磁導率為真空中的磁導率。盡管圖12所示模型較簡單，與實際地質情況不太相符，但僅用于對旅行時修正法進行驗證，簡單的模型可以足夠說明問題。

圖12 模型示意圖

模型的正演雷達波幅采用基于CPML邊界條件的有限差分正演[22]。圖13(a)所示為沒有對旅行時進行修正時根據圖12所示模型得出的速度層析成像圖；圖13(b)所示為對旅行時進行修正之后得到的速度層析成像圖。從圖13可以看出：進行旅行時修正之后的層析成像圖更加精確，能夠較清晰地判斷出異常體的形狀和位置，且能判斷出異常體內的波速。據圖13(a)能夠較清晰地判斷異常體的波速、形狀和位置，但并不能對模型中所有的異常體進行區(qū)分，不能反映出模型左上以及右下2個低速異常體。由此可以判斷：在鉆孔雷達的速度層析成像中旅行時是需要進行修正的，文中提出的修正方法也可取，得出的層析成像圖更加精確，對異常體的反映也更加敏感，能夠為判斷地下介質的結構特征等提供幫助。

(a) 未對旅行時進行修正；(b) 對旅行時進行修正

4 結論

1) 在針對大收發(fā)角度雷達波幅進行初至波提取時，因其具有很低的信噪比，從而較難提取出準確的旅行時。本文對同收發(fā)角度的雷達波幅進行排序、疊加得到參考波形，并對參考波形和該收發(fā)角度的雷達波幅進行互相關處理，從而得到每道雷達波幅的初至波旅行時。對合成雷達波幅的處理結果顯示，該方法對大收發(fā)角度雷達波幅的初至波提取具有很好的效果。

2) 以發(fā)射點位于鉆孔中間位置時為例，本文計算了色散條件下的電磁波波速與收發(fā)角度之間的關系，計算結果與模擬結果較吻合，證明了文中提出的電磁波傳播路徑更符合實際情況，即電磁波首先由發(fā)射天線中心點傳至發(fā)射天線末端，再由發(fā)射天線末端傳播至接收天線末端最后傳播至接收天線中心。因此，在針對鉆孔雷達雷達波幅進行層析成像時，需要對旅行時進行一定修正。

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]. 中國礦業(yè)大學學報, 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)對井中雷達的數(shù)值模擬[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] 李玉喜. 鉆孔雷達層析成像軟件系統(tǒng)的研究與開發(fā)[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. Constantattenuation 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.

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)

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.

borehole radar; first-arrival time; cross-correlation function; travel-time correction; tomography

10.11817/j.issn.1672-7207.2015.07.037

P631.3

A

1672?7207(2015)07?2658?07

2014?07?12；

2014?09?23

國家自然科學基金資助項目(41174061)；中南大學自由探索計劃項目(2011QNZT011) ( Project(41174061) supported by the National Natural Science Foundation of China; Project(2011QNZT011) supported by Exploration Program of Central South University)

彭凌星，博士研究生，從事地質和環(huán)境災害探測與評價；E-mail: innocent-eyes@163.com

(編輯 陳燦華)