探地雷達衰減補償逆時偏移成像方法

朱尉強， 黃清華

北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系, 北京　100871

?

探地雷達衰減補償逆時偏移成像方法

朱尉強， 黃清華*

北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系, 北京100871

探地雷達信號在地下介質(zhì)中傳播時易受到電導(dǎo)率所產(chǎn)生的衰減影響，從而使得傳統(tǒng)偏移成像結(jié)果在高衰減區(qū)域變得模糊.本文提出了衰減補償?shù)哪鏁r偏移方法來消除電導(dǎo)率的影響.該方法基于麥克斯韋方程組實現(xiàn)電磁波的正演模擬和逆時傳播.通過改變衰減項的正負號，保證了逆時傳播的時間對稱性，從而能夠重構(gòu)出原始波場，實現(xiàn)衰減補償.數(shù)值實驗比較了傳統(tǒng)逆時偏移方法和衰減補償逆時偏移方法在存在高導(dǎo)異常區(qū)域情況下的成像效果，結(jié)果證明了衰減補償逆時偏移方法能夠很好地恢復(fù)由電導(dǎo)率造成的信號衰減，從而提高探地雷達剖面的分辨率.關(guān)鍵詞探地雷達; 衰減補償; 逆時偏移

1　引言

探地雷達是基于高頻電磁波對淺層地下介質(zhì)信息進行探測的有效手段.由于電磁波與地震波間的相似性，地震數(shù)據(jù)處理方法被廣泛地運用于探地雷達數(shù)據(jù)的處理中，例如偏移成像技術(shù)(Fisher et al., 1992; Leuschen et al., 2001; Radzevicius 2008; Liu et al., 2014).但是這些探地雷達偏移方法都未考慮衰減的影響.對于電磁波而言，電導(dǎo)率的存在會造成振幅的衰減，從而導(dǎo)致在相對高導(dǎo)區(qū)域無法得到清晰的偏移成像結(jié)果，例如污水泄露區(qū)域，海水侵入環(huán)境等(Heteren et al., 1998; Sauck et al., 1998).電磁波的這一特點使得電導(dǎo)率成為探地雷達數(shù)據(jù)處理中必需考慮的因素.

對于探地雷達數(shù)據(jù)中的波形衰減效應(yīng)的研究目前主要有兩個方面：衰減的估計和衰減的補償.Turner et al.,(1994)與 Bradford (2007)提出了用Q值函數(shù)來描述探地雷達衰減效應(yīng).Turner (1994)與 Irving et al.(2003)借鑒地震學(xué)中的反Q濾波技術(shù)討論了對于探地雷達數(shù)據(jù)的Q值估計和反Q濾波方法，經(jīng)過反Q濾波后的數(shù)據(jù)被進一步用于偏移成像.張先武等(2014)通過地下介質(zhì)等效濾波器的振幅譜來求取反濾波器并對探地雷達數(shù)據(jù)進行反濾波處理消除衰減效應(yīng).但是反濾波方法使用了一維的衰減模型，因而無法考慮復(fù)雜的地質(zhì)條件.在偏移成像過程中考慮衰減補償則可適用于更為復(fù)雜的模型.電磁波的衰減發(fā)生于波場傳播過程，而偏移成像基于波場對時間的反向延拓，在此過程中同時考慮衰減的補償可以增強高衰減區(qū)域的信號強度，改善成像效果.現(xiàn)有的考慮衰減補償?shù)奶降乩走_偏移方法都為頻率域方法(Bano 1996; Bitri et al., 1998; Sena et al., 2006; Oden et al., 2007)，此類方法將衰減效應(yīng)作為復(fù)數(shù)波速的虛部進行補償.底青云等(2000)提出了基于2D有限元方法的GPR 偏移方法，考慮了電導(dǎo)率對于速度的影響，但未修正振幅的衰減.

逆時偏移方法是勘探地震學(xué)中針對復(fù)雜地下結(jié)構(gòu)的非常有效的偏移成像方法(Yoon et al., 2003; Etgen et al., 2009).該方法也被有效地應(yīng)用于探地雷達的數(shù)據(jù)處理中(Fisher et al., 1992; Leuschen et al., 2001; 雷林林等, 2015)，但之前的工作中都未考慮電導(dǎo)率的影響.借鑒勘探地震學(xué)中Q-RTM(Q-compensated reverse-time migration)方法(Zhu et al., 2014; Zhu 2014)，本文提出了針對探地雷達的衰減補償逆時偏移方法，通過改變衰減項的正負號保持了逆時傳播的時間對稱性，實現(xiàn)在逆時偏移過程中的衰減補償.通過數(shù)值實驗分析了不同逆時傳播方式與重構(gòu)波形之間的關(guān)系，驗證了本文方法相比傳統(tǒng)不考慮電導(dǎo)率的逆時偏移方法的優(yōu)勢.

2　方法

地下介質(zhì)的電導(dǎo)率會對電磁波的振幅與速度都產(chǎn)生影響.考慮平面波垂直入射均勻介質(zhì)的情況，設(shè)地表為xy平面，z軸垂直向下，根據(jù)麥克斯韋方程：

(1)

可以求出導(dǎo)體內(nèi)部電場為

(2)

(3)

其中，E為電場分量，ω為角頻率，μ為磁導(dǎo)率，ε為介電常數(shù)，σ為電導(dǎo)率.從方程(2)可以看出在電導(dǎo)率對速度(v=ω/β)與衰減項e-αz都有一定影響.圖 1展示了不同電導(dǎo)率下，電磁波速度與振幅的變化，可以看到電導(dǎo)率主要造成電磁波振幅的衰減，而對速度的影響很小.當σ=0.01S·m-1，在深度1m處振幅已經(jīng)衰減為地表振幅的54%，而速度相比沒有電導(dǎo)率情況下僅減小約1%.由于電導(dǎo)率對于振幅的影響，對于存在高衰減區(qū)域的探地雷達剖面，衰減補償成為必要的數(shù)據(jù)處理方法.

圖1　電導(dǎo)率對電磁波速度和振幅的影響(a) 速度與電導(dǎo)率關(guān)系； (b) 相對振幅與電導(dǎo)率關(guān)系.地面電磁波振幅設(shè)為1.Fig.1　The effect of conductivity on velocity and amplitude(a) Relationship between velocity and conductivity; (b) relationship between relative amplitude and conductivity. The amplitude at surface is set as 1.

如果不考慮電導(dǎo)率項，電磁波方程與聲波方程則具有相同的形式：

(4)

(5)

其中，p為壓強，c為聲速.因而，傳統(tǒng)的探地雷達的逆時偏移方法都直接借鑒了勘探地震學(xué)中的逆時偏移方法.由于方程(4)與方程(5)僅僅與時間的二次導(dǎo)數(shù)相關(guān)，根據(jù)時間反轉(zhuǎn)不變原則(Fink1992),如果p(x,t)為聲波方程(4)的一個解，則p(x,-t)是同一問題的另外一個解.逆時偏移成像通過將接收信號時間反轉(zhuǎn) (t→-t)后作為邊界條件重新注入計算區(qū)域，從而實現(xiàn)波場的重構(gòu)，再通過成像條件提取出目標體的形態(tài)與位置(Kaelinetal., 2006).

但是對于包含電導(dǎo)率項的完整電磁波方程(方程(1))，存在對于時間的一階導(dǎo)數(shù)，失去了時間反轉(zhuǎn)不變性：

(6)

(7)

(8)

此時方程(8)與方程(1)在形式上完全一致，人為保持了時間的對稱性.從而基于方程(8)的波場逆時傳播可以完全重構(gòu)原始波場，實現(xiàn)對于衰減的補償.

3　數(shù)值算例

3.1逆時波場重構(gòu)

為了驗證逆時波場重構(gòu)的效果，我們采用了圖 2a所示均勻模型(εr=9,σ=0.01 S·m-1)進行驗證.我們在圖 2a中四個方向都放置了接收點(紅色)，圖中黃色星號為發(fā)射源.本文采用時間域有限差分方法(FDTD)來實現(xiàn)基于麥克斯韋方程的電磁波的模擬和逆時偏移.空間步長選為dx=dz=0.02 m，時間步長為dt=0.04 ns.激發(fā)信號采用了中心頻率為100 MHz的ricker子波.圖 2b為在模型中心點接收到的波形，展示了基于方程(6)與方程(8)下的逆時重構(gòu)的波場間差別.圖中藍線(No.1)為正演時波形，即在黃色星號點處輸入激發(fā)信號后記錄到的波場.粉線(No.2)，綠線(No.3)，紅色叉號(No.4)為逆時偏移時記錄到的波形，即將圖 2a中紅色接收點接收到的信號經(jīng)時間反轉(zhuǎn)后作為邊界條件注入計算區(qū)域后重構(gòu)的波場.其中，粉線為不考慮電導(dǎo)率時(εr=9,σ=0 S·m-1)的逆時偏移結(jié)果，綠線為基于方程(6)的逆時偏移結(jié)果，紅色叉號為基于方程(8)的逆時偏移結(jié)果.從圖 2b可以清楚地看到，傳統(tǒng)不考慮電導(dǎo)率的逆時偏移無法恢復(fù)電導(dǎo)率帶來的衰減，基于方程(6)的逆時偏移方法引入了二次衰減，相比傳統(tǒng)不考慮電導(dǎo)率方法在振幅上更加惡化了成像效果，而本文的衰減補償逆時偏移方法則可以精確地進行波場重構(gòu)，消除電導(dǎo)率帶來的振幅衰減.

圖2　(a) 測試電磁模型； (b) 不同逆時偏移方法重構(gòu)的電磁波振幅.圖2a中測試模型為均勻電磁模型(εr=9,σ=0.01 S·m-1)，黃色星號為發(fā)射源位置，紅色叉號為接收點位置.波場的正演模擬和逆時傳播都基于FDTD方法.圖 2b中藍線(No.1)為模型中心點記錄到的正演模擬波形.粉線(No.2)，綠線(No.3)，紅色叉號(No.4)為逆時傳播時記錄波形.其中，粉線為采用沒有電導(dǎo)率模型(εr=9,σ=0 S·m-1)的逆時傳播結(jié)果，綠線為基于方程(6)和圖2a中模型的逆時偏移結(jié)果，紅色叉號為基于方程(8)和圖2a中模型的逆時偏移結(jié)果.Fig.2　(a) Test model; (b) Reconstructed waveform by different reverse time modeling methods. The test model in Fig.2a is a homogeneous model with εr=9,σ=0.01 S·m-1. The yellow star is the location of the source and the red crosses are the locations of the receivers. FDTD method is used for the forward and time-reverse modeling of electromagnetic wavefields. The blue line (No.1) in Fig.2a is the waveform recorded at the center of the model during forward modeling. The pink line (No.2), green line (No.3), red crosses (No.4) are the recorded waveform during time-reverse modeling. The pink line is the reconstructed waveform produced by a model without conductivity (εr=9,σ=0 S·m-1). The green line is the reconstructed waveform using eq. (4) and the model in Fig.2a. The red crosses are the reconstructed waveform using eq. (6) and the model in Fig.2a.

3.2污水滲漏模型

為了驗證衰減補償逆時偏移方法相對于傳統(tǒng)不考慮衰減的逆時偏移方法的優(yōu)勢，我們首先設(shè)計了一個簡單的三層模型(圖3a):空氣層(εr=1)，泥土層(εr=9)和含水層(εr=80).在均勻電導(dǎo)率背景模型(σ=0.001 S·m-1)中間存在一高導(dǎo)區(qū)域(σ=0.01 S·m-1)，模擬生活中存在污水滲漏情況(Chang et al., 2004).為了僅驗證電導(dǎo)率的影響，我們只在電導(dǎo)率模型中放置了這一異常.高電導(dǎo)區(qū)域邊界通過高斯平滑處理以避免在邊界產(chǎn)生強烈反射.發(fā)射源間距為0.5 m,接受點間距為0.1 m.空間采樣，時間采樣，激發(fā)信號都與之前相同.一個均勻的相對介電常數(shù)模型(εr=9)和圖 3b的電導(dǎo)率模型被用于逆時偏移成像.

為了比較衰減補償逆時偏移方法的效果，我們首先對于不存在高電導(dǎo)率區(qū)域的模型進行了逆時偏移成像，結(jié)果如圖4a所示.在沒有高導(dǎo)異常影響情況下，界面可以被清晰地成像.然后我們針對存在高導(dǎo)異常的模型進行了傳統(tǒng)不考慮電導(dǎo)率的逆時偏移(圖4b)和衰減補償逆時偏移(圖4c).比較圖4中成像結(jié)果，由于高電導(dǎo)率導(dǎo)致的波形衰減，高導(dǎo)異常區(qū)域內(nèi)的界面信號被強烈衰減.通過衰減補償逆時偏移后，電導(dǎo)率導(dǎo)致的波形衰減被有效恢復(fù).

圖3　污水滲漏模型(a) 相對介電常數(shù)(εr)模型； (b) 電導(dǎo)率(σ)模型.相對介電常數(shù)模型包括三層:空氣層(εr=1)，泥土層(εr=9)和含水層(εr=80). 電導(dǎo)率模型背景為σ=0.001 S·m-1，中間具有存在一高導(dǎo)區(qū)域(σ=0.01 S·m-1)，模擬生活中存在污水滲漏情況.發(fā)射源與接收點都位于地面.圖中星號為發(fā)射源位置，間距為0.5 m,接收點間距為0.1 m.Fig.3　Wastewater infiltration model(a) model of dielectric permittivity; (b) model of conductivity. The model of dielectric permittivity consists of three layers: air, unsaturated soil, saturated soil. The background of the conductivity model is σ=0.001 S·m-1. There is a high conductive anomaly zone in the center, which corresponds to an area contaminated by waste water. The sources and receivers are distributed on the surface. The asterisks are locations of sources with a spacing 0.5 m; and the spacing of receivers is 0.1 m.

圖4　污水滲漏模型逆時偏移結(jié)果比較(a) 原始模型不包含電導(dǎo)率異常下的偏移結(jié)果； (b) 高電導(dǎo)異常下傳統(tǒng)不考慮電導(dǎo)率逆時偏移結(jié)果； (c) 高電導(dǎo)異常下衰減補償逆時偏移結(jié)果.圖(a)中逆時偏移結(jié)果是基于原始不含有電導(dǎo)率及相應(yīng)波形衰減影響的模擬數(shù)據(jù)得到，作為參照結(jié)果.圖(b)與圖(c)中逆時偏移結(jié)果都是基于圖 3模型的模擬數(shù)據(jù)得到.圖(b)與圖(c)中虛線為高導(dǎo)異常區(qū)域.Fig.4　Results of reverse time migration(a) Result using data without the conductivity anomaly; (b) result of conventional reverse time migration; (c) result of attenuation compensated reverse time migration. The result in Fig.4a which is for comparison is based on GPR data generated by a model without the conductivity anomaly. The results in Fig.4b and Fig.4c are based on GPR data generated by the model in Fig.3. The dashed lines in Fig.4b and Fig.4c show the locations of the high conductive anomaly zones.

3.3滲水斷層模型

為了進一步說明衰減補償逆時偏移方法在探地雷達中的應(yīng)用前景，我們設(shè)計了一個滲水的斷層模型(圖5)，斷層區(qū)域存在高電導(dǎo)率異常(σ=0.01 S·m-1，背景為σ=0.001 S·m-1).用于逆時偏移的模型通過對原始模型進行高斯平滑后得到(圖6).發(fā)射源間距為0.3 m，接收點間距為0.1 m.空間采樣，時間采樣，激發(fā)信號都與之前相同.圖7展示了傳統(tǒng)不考慮電導(dǎo)率的逆時偏移和衰減補償逆時偏移的結(jié)果.電磁波在經(jīng)過斷層內(nèi)部時，由于電導(dǎo)率的存在會使波形受到衰減.如果在偏移中不考慮電導(dǎo)率的影響，斷層內(nèi)部及其下部界面的成像結(jié)果(圖7中虛線方框區(qū)域)變得模糊.通過衰減補償逆時偏移后，斷層內(nèi)部及其下部界面的成像結(jié)果都被顯著改善.

3.4結(jié)果分析與討論

以上數(shù)值模擬實驗初步證明了本文采用的衰減補償逆時偏移成像方法在存在高電導(dǎo)率影響的探地雷達數(shù)據(jù)處理中所起作用.

在理想的采集方式下(圖2a)，即在源的四周都進行信號的采集記錄，通過衰減補償?shù)哪鏁r重構(gòu)可以完全恢復(fù)原始的波場(圖2b)，電導(dǎo)率造成的波形衰減被完全補償.這與第二節(jié)中理論分析相符.對于實際的測量，理想的采集方式很難實現(xiàn)，我們往往只能在探測區(qū)域的一側(cè)進行測量.但這并不影響逆時傳播中對波形的補償效果.之后的兩個數(shù)值模型結(jié)果(圖4與圖7)也證實了這一點.

在污水滲漏模型中，傳統(tǒng)的逆時偏移成像結(jié)果(圖4b)會使高導(dǎo)異常區(qū)域內(nèi)的信號變得較弱.在實際的數(shù)據(jù)中，如果存在一定噪聲，并且存在更多數(shù)量、形態(tài)更加復(fù)雜的地下構(gòu)造時就會對解釋造成很大的困難.經(jīng)過衰減逆時偏移成像后的界面信號強度與原始模型不包含電導(dǎo)率異常下的偏移結(jié)果相近，有效地消除電導(dǎo)率的影響，有利于數(shù)據(jù)的解釋分析.

滲水斷層模型雖然結(jié)構(gòu)較為復(fù)雜，但是通過使用衰減補償逆時偏移方法后仍然得到了很好的補償效果.圖7(b)與圖7(a)都為施加互相關(guān)成像條件后的直接結(jié)果，未添加增益.圖7(b)整體信號強度更大是衰減補償?shù)男Ч?圖7(b)與圖7(a)相比，被衰減補償?shù)男盘柤扔袛鄬觾?nèi)部的錯斷面，也有背景的層面.對于斷層下部的層面，電磁波在向下傳播或者遇到界面向上反射的過程中，會經(jīng)過斷層高導(dǎo)高衰減區(qū)域，從而造成衰減.相應(yīng)的，在逆時偏移的過程中，正傳波場或者逆時反傳波場在經(jīng)過斷層區(qū)域的時候就可以得到衰減補償，恢復(fù)真實的振幅.

相比傳統(tǒng)的探地雷達逆時偏移方法，衰減補償逆時偏移方法不需增加額外的步驟和計算量.但在實際數(shù)據(jù)處理中需要注意噪聲的影響.因為噪聲也會在衰減補償過程中被放大，特別是高頻噪聲.在強衰減的情況下，如果地面本身接收到的電磁信號已經(jīng)被衰減到遠低于噪聲水平，噪聲在逆時傳播后期的衰減補償過程中會被指數(shù)級放大.因而在處理實際數(shù)據(jù)時，必要的數(shù)據(jù)預(yù)處理是很有必要的，例如濾波處理等.對衰減補償程度加以限制，防止過度補償增強穩(wěn)定性也是另外一個改進方向.

本文采用的衰減補償逆時偏移方法理論上適用于任意復(fù)雜的電導(dǎo)率模型.但是在實際探地雷達數(shù)據(jù)處理中獲取較為精確的電導(dǎo)率模型是主要的限制因素.前人已經(jīng)在電法勘探中電阻率成像方法(Tripp et al., 1984; Sasaki 1994; Zhang et al., 1995; Auken et al., 2004)做了大量的工作，這些方法可被借鑒用于獲取衰減補償逆時偏移方法所需的電導(dǎo)率模型.井間探地雷達層析成像方法(Holliger et al., 2001; Chang et al., 2004; Giroux et al., 2007)與全波形反演方法(Ernst et al., 2007; Meles et al., 2010)也提供另外的獲取精確電導(dǎo)率模型的方法.

圖5　滲水斷層模型(a) 相對介電常數(shù)(εr)模型； (b) 電導(dǎo)率(σ)模型.模型模擬一個正斷層，地面為陡坎形狀，上部為空氣層，中間高導(dǎo)區(qū)域(σ=0.01 S·m-1)為滲水后的斷層破碎帶.同時我們在相對介電常數(shù)模型中放置了兩個高介電常數(shù)界面模擬內(nèi)部的錯斷面.Fig.5　Fluid-infiltrated fault zone model(a) model of dielectric permittivity; (b) model of conductivity. This model has the shape of a normal fault with scarp. The top layer is air and the high conductive area (σ=0.01 S·m-1) corresponds to the fault-fracture zone which is permeated with water. Two interface of high dielectric permittivity are placed inside the fault-fracture zone to simulate fault-planes.

圖6　用于逆時偏移的模型(a) 相對介電常數(shù)(εr)模型； (b) 電導(dǎo)率(σ)模型.該模型由圖 5中模型經(jīng)過高斯平滑后得到.Fig.6　Models for reverse time migration(a) model of dielectric permittivity; (b) model of conductivity. These models are generated by a Gaussian smoothing on the models in Fig.5.

圖7　滲水斷層逆時偏移結(jié)果比較(a) 傳統(tǒng)不考慮電導(dǎo)率逆時偏移結(jié)果； (b) 衰減補償逆時偏移結(jié)果.Fig.7　Results of reverse time migration(a) Result of conventional reverse time migration; (b) Result of attenuation compensated reverse time migration.

4　結(jié)論

衰減補償逆時偏移成像方法基于完整的電磁波方程，通過改變衰減項的正負號，保持了波動方程在逆時傳播時的時間對稱性，從而實現(xiàn)了對電導(dǎo)率造成的振幅衰減的補償.通過與傳統(tǒng)不考慮電導(dǎo)率影響的逆時偏移方法的比較，證明了衰減補償逆時偏移成像能夠很好地重構(gòu)原始波場，恢復(fù)高衰減區(qū)域的信號強度，提高了成像的精度.

致謝感謝審稿專家和編輯部老師的支持和幫助.

Bradford J H. 2007. Frequency-dependent attenuation analysis of ground-penetrating radar data.Geophysics, 72(3):J7-J16.

Chang P Y, Alumbaugh D, Brainard J, et al. 2004. The application of ground penetrating radar attenuation tomography in a vadose zone infiltration experiment.JournalofContaminantHydrology, 71(1-4):67-87.Di Q Y, Xu K, Wang M Y. 2000. The attenuated radar wave migration with finite element method.ChineseJournalofGeophysics(in Chinese), 43(2):257-263.

Ernst J R, Maurer H, Green A G, et al. 2007. Full-Waveform Inversion of Crosshole Radar Data Based on 2-D Finite-Difference Time-Domain Solutions of Maxwell′s Equations.IEEETransactionsonGeoscienceandRemoteSensing, 45(9):2807-2828.

Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics.Geophysics, 74(6):WCA5-WCA17.Fink M. 1992. Time reversal of ultrasonic fields. I. Basic principles. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 39(5):555-566.Fisher E, McMechan G A, Annan A P, et al. 1992. Examples of reverse-time migration of single-channel, ground-penetrating radar profiles.Geophysics, 57(4):577-586.

Giroux B, Gloaguen E, Chouteau M. 2007. bh_tomo—a Matlab borehole georadar 2D tomography package.Computers&Geosciences, 33(1):126-137.

Heteren S V, Fitzgerald D M, Mckinlay P A, et al. 1998. Radar facies of paraglacial barrier systems: coastal New England, USA.Sedimentology, 45(1):181-200.

Holliger K, Musil M, Maurer H R. 2001. Ray-based amplitude tomography for crosshole georadar data: a numerical assessment.JournalofAppliedGeophysics, 47(3):285-298.Irving J D, Knight R J. 2003. Removal of wavelet dispersion from ground-penetrating radar data.Geophysics, 68(3):960-970.

Kaelin B, Guitton A. 2006. Imaging condition for reverse time migration. In 2006 SEG Annual Meeting. Society of Exploration Geophysicists.Lei L L, Liu S X, Fu L, et al. 2015. Reverse time migration applied to GPR data based on full wave inversion.ChineseJournalofGeophysics(in Chinese), 58(9):3346-3355.

Leuschen C J, Plumb R G. 2001. A matched-filter-based reverse-time migration algorithm for ground-penetrating radar data.IEEETransactionsonGeoscienceandRemoteSensing, 39(5):929-936.

Liu S, Lei L, Fu L, et al. 2014. Application of pre-stack reverse time migration based on FWI velocity estimation to ground penetrating radar data.JournalofAppliedGeophysics, 107:1-7.

Meles G A, Van Der Kruk J, Greenhalgh S A, et al. 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to- surface GPR data.IEEETransactionsonGeoscienceandRemoteSensing, 48(9):3391-3407.

Oden C P, Powers M H, Wright D L, et al. 2007. Improving GPR image resolution in lossy ground using dispersive migration.IEEETransactionsonGeoscienceandRemoteSensing, 45(8):2492-2499.

Radzevicius S. 2008. Practical 3-d migration and visualization for accurate imaging of complex geometries with gpr.JournalofEnvironmental&EngineeringGeophysics, 13(2):99-112.

Sasaki Y. 1994. 3-D resistivity inversion using the finite-element method.Geophysics, 59(12):1839-1848.

Sauck W A, Atekwana E A, Nash M S. 1998. High Conductivities Associated With an LNAPL Plume Imaged by Integrated Geophysical Techniques.JournalofEnvironmentalandEngineeringGeophysics, 2:203-212.

Sena A R, Stoffa P L, Sen M K. 2006. Split-step Fourier migration of GPR data in lossy media.Geophysics, 71(4):K77-K91.

Tripp A C, Hohmann G W, Swift C M. 1984. Two-dimensional resistivity inversion.Geophysics, 49(10):1708-1717.Turner G. 1994. Subsurface radar propagation deconvolution.Geophysics, 59(2):215-223.

Turner G, Siggins A F. 1994. Constant Q attenuation of subsurface radar pulses.Geophysics, 59(8):1192-1200.

Yoon K, Shin C, Suh S, et al. 2003. 3D reverse-time migration using the acoustic wave equation: An experience with the SEG/EAGE data set.TheLeadingEdge, 22(1):38-41.

Zhang J, Mackie R L, Madden T R. 1995. 3-D resistivity forward modeling and inversion using conjugate gradients.Geophysics, 60(5):1313-1325.

Zhang X W, Gao Y Z, Fang G Y. 2014. An inverse filtering method for removing the wavelet attenuation and dispersion of Ground Penetrating Radar data.ChineseJournalofGeophysics(in Chinese), 57(3):932-938.

Zhu T. 2014. Time-reverse modelling of acoustic wave propagation in attenuating media.GeophysicalJournalInternational, 197(1):483-494.

Zhu T, Harris J, Biondi B. 2014. Q-compensated reverse-time migration.Geophysics, 79(3):S77-S87.

附中文參考文獻

底青云，許琨，王妙月. 2000. 衰減雷達波有限元偏移. 地球物理學(xué)報, 43(2):257-263.

雷林林，劉四新，傅磊等. 2015. 基于全波形反演的探地雷達數(shù)據(jù)逆時偏移成像. 地球物理學(xué)報, 58(9):3346-3355.

張先武，高云澤，方廣有. 2014. 消除探地雷達數(shù)據(jù)的子波衰減和頻散的反濾波方法. 地球物理學(xué)報, 57(3):932-938.

(本文編輯劉少華)

Attenuation compensated reverse time migration method of ground penetrating radar signals

ZHU Wei-Qiang, HUANG Qing-Hua*

DepartmentofGeophysics,SchoolofEarthandSpaceSciences,PekingUniversity,Beijing100871,China

Conductivity of subsurface media causes attenuation of ground penetrating radar (GPR) signals, so that the imaging results of conventional migration methods are blurred at strong attenuation zones. We proposed an attenuation compensated reverse time migration method for GPR signals. The forward and time-reverse propagation of electromagnetic waves was based on the Maxwell′s equations. The sign of the conductive term was reversed to keep equations′ temporal symmetry, so that the original wavefield could be reconstructed and the attenuation effect was compensated. The numerical experiments further compared the results of conventional reverse time migration and attenuation compensated reverse time migration. The imaging results showed that attenuation compensated method recovered the weak signals within and beneath high conductivity zones and improved the resolution of GPR profiles.

Ground penetrating radar; Attenuation compensation; Reverse time migration

Auken E, Christiansen A V. 2004. Layered and laterally constrained 2D inversion of resistivity data.Geophysics, 69(3):752-761.Bano M. 1996. Modeling and inverseQimaging of ground penetrating radar waves in 1 and 2D.GeophysicalResearchLetters, 23(22):3123-3126.Bitri A, Grandjean G. 1998. Frequency-wavenumber modelling and migration of 2D GPR data in moderately heterogeneous dispersive media.GeophysicalProspecting, 46(3):287-301.

10.6038/cjg20161034.

國家自然科學(xué)基金項目(41574104,41274075)資助.

朱尉強，男，1989年生，2013年畢業(yè)于北京大學(xué)，現(xiàn)北京大學(xué)地球物理專業(yè)研究生在讀.

黃清華，男，教授，1990年畢業(yè)于中國科學(xué)技術(shù)大學(xué)，1999年獲日本大阪大學(xué)博士學(xué)位.主要從事地球電磁學(xué)、地震物理學(xué)方面的教學(xué)和科研工作，E-mail：huangq@pku.edu.cn

10.6038/cjg20161034

P631

2015-12-21，2016-03-19收修定稿

朱尉強， 黃清華. 2016. 探地雷達衰減補償逆時偏移成像方法. 地球物理學(xué)報，59(10):3909-3916,

Zhu W Q, Huang Q H. 2016. Attenuation compensated reverse time migration method of ground penetrating radar signals.ChineseJ.Geophys. (in Chinese),59(10):3909-3916,doi:10.6038/cjg20161034.