基于歸一化互相關(guān)成像條件的GPR逆時(shí)偏移成像

陳德鵬，戴前偉，馮德山，王洪華，張彬

?

基于歸一化互相關(guān)成像條件的GPR逆時(shí)偏移成像

陳德鵬1, 2, 3，戴前偉1, 2，馮德山1, 2，王洪華3，張彬1, 2

(1. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙，410083;2. 有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410083；3. 桂林理工大學(xué) 地球科學(xué)學(xué)院，廣西 桂林，541004)

針對(duì)基于零時(shí)刻成像條件的探地雷達(dá)(GPR)逆時(shí)偏移精度低、難以對(duì)復(fù)雜結(jié)構(gòu)進(jìn)行高精度成像的缺點(diǎn)，將歸一化互相關(guān)成像條件應(yīng)用于GPR疊前逆時(shí)偏移成像。從二維GPR電磁波方程出發(fā)，采用基于完全匹配層(PML)邊界條件的時(shí)域有限單元法(FETD)模擬電磁波正向和逆向傳播，采用歸一化互相關(guān)成像條件獲取疊前逆時(shí)偏移的偏移結(jié)果，將空間高通濾波用于壓制互相關(guān)過程中產(chǎn)生的低頻噪聲，然后編制相應(yīng)的GPR疊前逆時(shí)偏移程序。在此基礎(chǔ)上，建立2個(gè)復(fù)雜的GPR模型，利用基于歸一化互相關(guān)成像條件的GPR逆時(shí)偏移程序進(jìn)行計(jì)算，并與基于零時(shí)刻成像條件的GPR逆時(shí)偏移剖面進(jìn)行對(duì)比。研究結(jié)果表明：與基于零時(shí)刻成像條件的GPR逆時(shí)偏移剖面相比，基于歸一化互相關(guān)成像條件的GPR疊前逆時(shí)偏移剖面能更清晰地反映異常體空間形態(tài)和內(nèi)部結(jié)構(gòu)信息，其分辨率和成像精度更高。

探地雷達(dá)；逆時(shí)偏移；歸一化互相關(guān)成像條件；零時(shí)刻成像條件

探地雷達(dá)(ground penetrating radar, GPR)作為一種采用高頻脈沖電磁波探測(cè)地下地質(zhì)體結(jié)構(gòu)及物性參數(shù)的重要淺部地球物理方法[1]，具有無損性、分辨率和效率均較高等優(yōu)點(diǎn)，在工程與環(huán)境等領(lǐng)域得到廣泛應(yīng)用[2?4]。偏移是實(shí)測(cè)GPR數(shù)據(jù)向真實(shí)地下結(jié)構(gòu)轉(zhuǎn)換的最重要環(huán)節(jié)，偏移質(zhì)量將直接影響后續(xù)解釋的準(zhǔn)確性。目前，常用的GPR偏移方法主要有Kirchhoff偏移 法[5]、頻率?波數(shù)域偏移法[6]和逆時(shí)偏移法[7]等。在眾多偏移方法中，逆時(shí)偏移方法具有精度高、相位準(zhǔn)確、不受橫向變速和高陡傾角影響、能夠適應(yīng)沿任意方向傳播波場(chǎng)等優(yōu)點(diǎn)，可更精準(zhǔn)地將電磁波場(chǎng)歸位到其真實(shí)位置，因此，得到越來越多研究者的廣泛關(guān)注[8]。與彈性波逆時(shí)偏移類似，GPR逆時(shí)偏移根據(jù)不同的數(shù)據(jù)采集方式可采用不同的成像條件。目前，常用的GPR逆時(shí)偏移成像條件主要有零時(shí)刻成像條件、零延時(shí)互相關(guān)成像條件和歸一化互相關(guān)成像條件。零時(shí)刻成像條件的基本原理是根據(jù)爆炸反射面模型，將地表接收的GPR數(shù)據(jù)進(jìn)行逆時(shí)外推，接收波場(chǎng)外推至零時(shí)刻的結(jié)果即為逆時(shí)偏移成像剖面。目前，該成像條件主要用于等偏移距GPR數(shù)據(jù)逆時(shí)偏移處理。FISHER等[9]借鑒彈性波逆時(shí)偏移算法，將零時(shí)刻成像條件應(yīng)用于GPR數(shù)據(jù)的逆時(shí)偏移成像。隨后，LEUSCHEN等[10?14]對(duì)零時(shí)刻成像條件的GPR逆時(shí)偏移算法進(jìn)行完善，并應(yīng)用于實(shí)測(cè)數(shù)據(jù)處理。由于現(xiàn)有的商業(yè)GPR系統(tǒng)大都采用剖面法測(cè)量，因此，基于零時(shí)刻成像條件的逆時(shí)偏移算法被廣泛應(yīng)用于實(shí)測(cè)GPR數(shù)據(jù)處理。然而，采用剖面法測(cè)量時(shí)，偏移距較小，難以對(duì)高陡地層構(gòu)造產(chǎn)生的反射波進(jìn)行有效采集，因此，基于零時(shí)刻成像條件的GPR逆時(shí)偏移難以對(duì)復(fù)雜結(jié)構(gòu)進(jìn)行高精度成像。近年來，隨著陣列天線GPR系統(tǒng)的成功研制和應(yīng)用，多偏移距數(shù)據(jù)的快速采集成為可能，GPR疊前逆時(shí)偏移方法得到人們的重視。它以采集的多偏移距數(shù)據(jù)為基礎(chǔ)，根據(jù)成像點(diǎn)的入射波到時(shí)和反射波到時(shí)相同的原理，計(jì)算源點(diǎn)正傳電磁波場(chǎng)和接收點(diǎn)反傳電磁波場(chǎng)的互相關(guān)，從而對(duì)地下復(fù)雜結(jié)構(gòu)進(jìn)行成像。朱尉強(qiáng)等[8, 15?16]將基于互相關(guān)成像條件的逆時(shí)偏移算法引入GPR數(shù)據(jù)處理，針對(duì)GPR電磁波的傳播特點(diǎn)和介質(zhì)特性對(duì)逆時(shí)偏移算法進(jìn)行完善，能適應(yīng)多天線或陣列天線觀測(cè)模式。然而，現(xiàn)有的互相關(guān)成像條件應(yīng)用于GPR疊前逆時(shí)偏移處理時(shí)，在發(fā)射源點(diǎn)附近會(huì)產(chǎn)生很強(qiáng)的低頻噪聲，并且深部結(jié)構(gòu)的反射能量異常微弱。為此，CLAERBOUT等[17]根據(jù)彈性波逆時(shí)偏移算法提出了歸一化互相關(guān)成像條件，其基本原理是利用源點(diǎn)電磁波場(chǎng)或接收點(diǎn)電磁波場(chǎng)對(duì)兩者的互相關(guān)成像結(jié)果進(jìn)行歸一化。該成像條件可有效減弱發(fā)射源點(diǎn)附近的低頻噪聲，且大大加強(qiáng)深部結(jié)構(gòu)的反射能量，目前被廣泛應(yīng)用于聲波和彈性波逆時(shí)偏移成像。本文在借鑒前人逆時(shí)偏移理論基礎(chǔ)上，將歸一化互相關(guān)成像條件應(yīng)用于GPR疊前逆時(shí)偏移成像。將時(shí)間域有限單元法(FETD)用于模擬電磁波正向和逆向傳播過程，歸一化互相關(guān)成像條件用于獲取逆時(shí)偏移結(jié)果。設(shè)計(jì)2個(gè)復(fù)雜GPR模型，利用歸一化互相關(guān)成像條件的逆時(shí)偏移算法進(jìn)行計(jì)算，并與零時(shí)刻成像條件的逆時(shí)偏移結(jié)果進(jìn)行對(duì)比。

1 GPR逆時(shí)偏移成像方法

GPR逆時(shí)偏移的基本原理是將地表采集的接收電磁波場(chǎng)在時(shí)間軸上進(jìn)行逆向傳播，當(dāng)電磁波場(chǎng)逆推至零時(shí)刻時(shí)，所有反射波與繞射波的能量都回到最初被反射和繞射的空間位置，然后應(yīng)用成像條件可獲取最終的成像剖面。因此，GPR逆時(shí)偏移實(shí)現(xiàn)過程可分為以下3步：1) 在時(shí)間域進(jìn)行源點(diǎn)電磁波場(chǎng)的正向傳播；2) 在時(shí)間域進(jìn)行接收點(diǎn)電磁波場(chǎng)逆時(shí)外推；3) 利用成像條件獲得成像結(jié)果。其中，前2步可利用相同的數(shù)值模擬方法進(jìn)行計(jì)算，本文采用基于三角形網(wǎng)格剖分的FETD來模擬電磁波的正向及逆時(shí)外推。

1.1 電磁波場(chǎng)正向傳播及逆向外推

考慮二維地電模型的走向?yàn)檩S，根據(jù)電磁波理論，TM模式下電場(chǎng)分量滿足的波動(dòng)方程為[18]

式中：E為方向電場(chǎng)強(qiáng)度(V/m)；，和0分別為介質(zhì)的介電常數(shù)(F/m)、電導(dǎo)率(S/m)和真空磁導(dǎo)率(H/m)；為激勵(lì)源函數(shù)。利用伽遼金有限元法可推導(dǎo)式(1)的有限元方程[18]為

其中：，′和分別為質(zhì)量矩陣、阻尼矩陣和剛度矩陣，它們只與介質(zhì)的物性參數(shù)和幾何分布有關(guān)；為源向量。采用三角形網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行剖分(如圖1所示，其中N和N分別為和方向的網(wǎng)格數(shù))，每個(gè)矩形網(wǎng)格被2條對(duì)角線再剖分成4個(gè)三角形網(wǎng)格，以提高復(fù)雜GPR模型的正向與逆向電磁波場(chǎng)的計(jì)算精度。同時(shí)，在計(jì)算區(qū)域的外邊界采用完全匹配層吸收邊界條件[19]處理截?cái)噙吔缣幍某瑥?qiáng)反射波。

圖1 有限元網(wǎng)格剖分及節(jié)點(diǎn)編號(hào)示意圖

1.2 空間高通濾波

在實(shí)現(xiàn)逆時(shí)偏移過程中，源點(diǎn)電磁波場(chǎng)和接收點(diǎn)電磁波場(chǎng)的互相關(guān)計(jì)算不僅會(huì)在反射界面處產(chǎn)生強(qiáng)振幅，而且會(huì)在電磁波傳播整個(gè)路徑的非反射點(diǎn)處產(chǎn)生強(qiáng)振幅，若將這些振幅沿著時(shí)間求和，則會(huì)在剖面頂部出現(xiàn)弧形強(qiáng)干擾現(xiàn)象，特別是電性差異的界面會(huì)形成呈低波數(shù)特征的串?dāng)_噪聲，這些低頻噪聲的存在嚴(yán)重影響成像質(zhì)量。為有效消除這些低頻噪聲，目前主要采用的去燥方法有：在電磁波場(chǎng)傳播過程中消除；應(yīng)用去噪成像條件，成像后采用濾波法去噪等[20]。其中，采用空間高通濾波法壓制低頻噪聲，具有理論簡(jiǎn)單、易操作的優(yōu)點(diǎn)，被廣泛應(yīng)用于GPR逆時(shí)偏移成像去燥。為此，本文利用空間高通濾波法[21]對(duì)成像結(jié)果進(jìn)行去燥處理，其高通濾波模板計(jì)算公式為

2 歸一化互相關(guān)成像條件

成像條件的選取是開展高效高精度GPR逆時(shí)偏移的關(guān)鍵之一，其選取質(zhì)量直接影響成像質(zhì)量。目前，GPR逆時(shí)偏移中常用的成像條件主要3種：零時(shí)刻成像條件、零延時(shí)互相關(guān)成像條件和歸一化互相關(guān)成像條件。

2.1 零時(shí)刻成像條件

零時(shí)刻成像條件主要應(yīng)用于等偏移距GPR數(shù)據(jù)的逆時(shí)偏移成像。根據(jù)爆炸反射面原理，將等偏移距GPR數(shù)據(jù)作為邊值條件，并以介質(zhì)半速度進(jìn)行逆時(shí)外推。傳播到零時(shí)刻的電磁波場(chǎng)就是爆炸反射面的 像[22]，其成像表達(dá)式為

其中：(，)為逆時(shí)偏移成像剖面；E(，，)為地表接收到GPR數(shù)據(jù)。采用剖面法進(jìn)行GPR數(shù)據(jù)采集時(shí)，偏移距一般都較小，難以對(duì)高陡地層構(gòu)造產(chǎn)生的反射波進(jìn)行有效采集，因此，基于零時(shí)刻成像條件的GPR逆時(shí)偏移算法難以對(duì)復(fù)雜結(jié)構(gòu)進(jìn)行高精度成像。

2.2 歸一化互相關(guān)成像條件

互相關(guān)成像條件最早在彈性波逆時(shí)偏移領(lǐng)域被提出，其基本原理[17]是以時(shí)間一致性原理為基礎(chǔ)，通過計(jì)算源點(diǎn)正傳電磁波場(chǎng)和接收點(diǎn)反傳電磁波場(chǎng)的互相關(guān)，從而獲得成像剖面。標(biāo)準(zhǔn)的零延遲互相關(guān)成像條件公式[7]為

由式(6)和(7)可知：歸一化互相關(guān)成像條件是在零延時(shí)互相關(guān)成像條件的基礎(chǔ)上，用源點(diǎn)或接收點(diǎn)電磁波場(chǎng)進(jìn)行歸一化，可削弱源強(qiáng)振幅的影響，提升深部反射界面的能量。CHATTOPADHYAY等[23]在彈性波逆時(shí)偏移中對(duì)上述成像條件進(jìn)行對(duì)比分析，認(rèn)為源歸一化互相關(guān)成像條件能夠更精確反映地下復(fù)雜結(jié)構(gòu)，可提供較準(zhǔn)確的成像結(jié)果。因此，本文采用源歸一化互相關(guān)成像條件進(jìn)行GPR疊前逆時(shí)偏移成像。

3 數(shù)值算例

3.1 模型一

圖2所示為1.80 m×1.25 m(長(zhǎng)×寬)的三圓GPR模型示意圖，在埋深約0.30 m處存在1條起伏界面。起伏界面上方介質(zhì)的相對(duì)介電常數(shù)為5，電導(dǎo)率為0.001 S/m；下方介質(zhì)的相對(duì)介電常數(shù)為10，電導(dǎo)率為0.000 5 S/m。起伏界面的下方埋有3個(gè)半徑為0.05 m的圓形空洞異常體，其圓心位置分別為(0.40, 0.65)，(0.90, 0.90)和(1.40, 1.15) m。利用FETD進(jìn)行等偏移距和多偏移距測(cè)量方式進(jìn)行正演計(jì)算時(shí)的發(fā)射源信號(hào)采用中心頻率為500 MHz的雷克子波，采樣時(shí)間間隔為0.01 ns，時(shí)窗長(zhǎng)度為35 ns。采用等偏移距測(cè)量方式正演時(shí)，收、發(fā)天線的距離為0.01 m。采用多偏移距測(cè)量方式正演時(shí)，在地表布設(shè)15個(gè)發(fā)射天線，發(fā)射天線的間距為0.10 m，第1個(gè)發(fā)射天線放置在(0.1,0) m位置處；每個(gè)發(fā)射天線右側(cè)布置50道接收天線，接收點(diǎn)的間距為0.01 m，水平方向上覆蓋長(zhǎng)度為0.50 m。

圖2 GPR模型一示意圖

圖3(a)和(b)所示分別為模型一的多偏移距和等偏移距正演剖面。由圖3可見：起伏界面和3個(gè)圓狀異常體的反射波清晰，但這3個(gè)圓狀異常體產(chǎn)生的反射波不能準(zhǔn)確定位異常體的真實(shí)空間位置。由此可見，開展GPR逆時(shí)偏移成像準(zhǔn)確定位異常體的真實(shí)空間位置非常必要。

圖3 模型一GPR正演剖面

圖4(a)和圖4(b)所示分別為采用基于歸一化互相關(guān)成像條件和零時(shí)刻成像條件的GPR逆時(shí)偏移算法對(duì)圖3(a)和圖3(b)中的數(shù)據(jù)進(jìn)行逆時(shí)偏移，獲得成像剖面。由圖4可見：這2種成像條件的逆時(shí)偏移獲得的剖面中起伏界面形態(tài)都非常清晰，繞射波收斂且與真實(shí)位置相符。由于模型兩端在基于歸一化互相關(guān)成像條件的疊前逆時(shí)偏移成像過程中疊加次數(shù)為1，而中間區(qū)域進(jìn)行了多次疊加，因此，界面兩端的成像不如中間界面成像清晰。偏移剖面中3個(gè)圓狀異常體產(chǎn)生的繞射波大部分收斂，并歸位到其真實(shí)位置。但是對(duì)比圖3(a)和圖4(b)可知：圖4(a)中圓狀異常體的疊前逆時(shí)偏移結(jié)果能清晰地反映空洞模型的空間形態(tài)和內(nèi)部結(jié)構(gòu)信息，而圖4(b)中基于零時(shí)刻成像條件的逆時(shí)偏移結(jié)果中3個(gè)圓狀異常體的繞射波能量大都收斂于圓形位置，不能清晰地反映空洞模型的空間分布特征。

為了更好地分析歸一化互相關(guān)成像條件的優(yōu)勢(shì)，分別提取歸一化成像條件和零時(shí)刻條件的逆時(shí)偏移結(jié)果在水平位置0.40 m處的單道波形，如圖5所示。由圖5可知：基于歸一化成像條件的逆時(shí)偏移結(jié)果的單道波形在空洞異常體的頂部和底部都能很好地成像，能大致圈定異常體的垂直分布范圍；而基于零時(shí)刻成像條件逆時(shí)偏移結(jié)果的單道波形只能對(duì)空洞異常體的底部很好地成像，并且地表產(chǎn)生的直達(dá)波非常強(qiáng)烈，壓制了后續(xù)異常體成像的振幅。對(duì)比結(jié)果表明：相比于零時(shí)刻成像條件的逆時(shí)偏移結(jié)果，基于互相關(guān)成像條件的逆時(shí)偏移成像結(jié)果分辨率更高，成像質(zhì)量更好。

(a) 歸一化互相關(guān)成像條件；(b) 零時(shí)刻成像條件

1—?dú)w一化互相關(guān)成像條件；2—零時(shí)刻成像條件。

Fig. 5 Comparison of single waveform of RTM profiles of model 1 at horizontal position of 0.40 m

3.2 模型二

模型二為1.80 m×1.25 m(長(zhǎng)×寬)的矩形區(qū)域，埋深約0.30 m處存在1條起伏界面，如圖6所示。起伏界面上方和下方介質(zhì)的相對(duì)介電常數(shù)和電導(dǎo)率與模型一的相同。起伏界面下方分別埋有1個(gè)半徑為0.05 m的圓狀空洞異常體和1個(gè)邊長(zhǎng)為0.10 m的正方形空洞異常體，圓狀空洞異常體的圓心位置為(0.40, 0.65) m，正方形空洞異常體的中心位置為(1.40, 0.65) m。利用FETD進(jìn)行等偏移距和多偏移距測(cè)量方式進(jìn)行正演計(jì)算時(shí)的發(fā)射源信號(hào)采用中心頻率為500 MHz的雷克子波，采樣時(shí)間間隔為0.01 ns，時(shí)窗長(zhǎng)度為35 ns。采用等偏移距測(cè)量方式正演時(shí)，收發(fā)天線的距離為 0.01 m。采用多偏移距測(cè)量方式正演時(shí)，在地表布設(shè)15個(gè)發(fā)射天線，發(fā)射天線的間距為0.10 m，第1個(gè)發(fā)射天線放置在(0.10,0) m位置；每個(gè)發(fā)射天線右側(cè)布置50道接收天線，接收點(diǎn)的間距為0.01 m，水平方向上覆蓋長(zhǎng)度為0.50 m。

圖6 GPR模型二示意圖

圖7(a)和圖7(b)所示分別為模型二的多偏移距和等偏移距正演剖面。由圖7可見：起伏界面、圓狀異常體和正方形空洞異常體的反射波清晰可見，但是圓狀異常體和正方形空洞產(chǎn)生的反射波不能準(zhǔn)確定位其真實(shí)空間位置。為此，對(duì)圖7中的多偏移距和等偏移距數(shù)據(jù)分別利用基于歸一化互相關(guān)成像條件和基于零時(shí)刻成像條件的GPR逆時(shí)偏移算法進(jìn)行逆時(shí)偏移成像。圖8(a)和圖8(b)所示為分別利用圖7(a)和圖7(b)中的數(shù)據(jù)進(jìn)行GPR逆時(shí)偏移成像的結(jié)果。

由圖8可見：偏移剖面中3個(gè)圓狀異常體產(chǎn)生的繞射波大部分收斂，并歸位到其真實(shí)位置。但對(duì)比圖8(a)和圖8(b)可知：圖8(a)中圓狀異常體的疊前逆時(shí)偏移結(jié)果能清晰地反映空洞模型的空間形態(tài)和內(nèi)部結(jié)構(gòu)信息，而圖8(b)中疊后逆時(shí)偏移結(jié)果中3個(gè)圓狀異常體的繞射波能量大都收斂于圓形位置，不能清晰地反映空洞模型的空間分布特征。圖9所示為2種成像條件的逆時(shí)偏移結(jié)果在水平位置0.40 m處的單道波形。由圖9可知：基于歸一化成像條件的逆時(shí)偏移結(jié)果的單道波形在空洞異常體的頂部和底部都能很好地成像，能大致圈定異常體的垂直分布范圍；而基于零時(shí)刻成像條件逆時(shí)偏移結(jié)果的單道波形只能對(duì)空洞異常體的底部很好地成像，并且地表產(chǎn)生的直達(dá)波非常強(qiáng)烈，壓制了后續(xù)異常體成像的振幅。對(duì)比結(jié)果表明：相比于零時(shí)刻成像條件的逆時(shí)偏移結(jié)果，基于互相關(guān)成像條件的逆時(shí)偏移成像結(jié)果的分辨率更高，成像質(zhì)量更好。

(a) 多偏移距GPR正演剖面；(b) 等偏移距GPR正演剖面

(a) 歸一化互相關(guān)成像條件；(b) 零時(shí)刻成像條件

1—?dú)w一化互相關(guān)成像條件；2—零時(shí)刻成像條件。

4 結(jié)論

1) 為提高GPR數(shù)據(jù)的逆時(shí)偏移計(jì)算精度，將歸一化互相關(guān)成像條件引入到GPR的逆時(shí)偏移成像中。

2) 采用基于三角形剖分的FETD模擬電磁波場(chǎng)的正向及逆向傳播，可獲得基于歸一化互相關(guān)成像條件的GPR逆時(shí)偏移成像剖面。

3) 基于歸一互相關(guān)成像條件的GPR逆時(shí)偏移方法能更精細(xì)地刻畫出模型中異常體的空間形態(tài)和內(nèi)部結(jié)構(gòu)特征，提高了GPR對(duì)異常體的探測(cè)分辨率。

[1] 劉瀾波, 錢榮毅. 探地雷達(dá): 淺表地球物理科學(xué)技術(shù)中的重要工具[J]. 地球物理學(xué)報(bào), 2015, 58(8): 2606?2617.LIU Lanbo, QIAN Rongyi. Ground penetrating radar: a critical tool in near-surface geophysics[J]. Chinese Journal of Geophysics, 2015, 58(8): 2606?2617.

[2] ALANI A M, ABOUTALEBI M, KILIC G. Applications of ground penetrating radar(GPR) in bridge deck monitoring and assessment[J]. Journal of Applied Geophysics, 2013, 97(13): 45?54.

[3] CHANG C W, LIN C H, LIEN H S. Measurement radius of reinforcing steel bar in concrete using digital image GPR[J]. Construction and Building Materials, 2009, 23(2): 1057?1063.

[4] SOLLA M, LORENZO H, RIAL F I, et al. GPR evaluation of the Roman masonry arch bridge of Lugo (Spain)[J]. Nondestructive Testing and Evaluation International, 2011, 44(1): 8?12.

[5] 戴前偉, 馮德山, 何繼善. Kirchhoff偏移法在探地雷達(dá)正演圖像處理中的應(yīng)用[J].地球物理學(xué)進(jìn)展, 2005, 20(3): 849?853. DAI Qianwei, FENG Deshan, HE Jishan. The application of Kirchhoff's migration method in the image processing of the ground penetrating radar forward simulation[J]. Progress in Geophysics, 2005, 20(3): 849?853.

[6] BITRI A, GRANDJEAN G. Frequency-wavenumber modelling and migration of 2D GPR data in moderately heterogeneous dispersive media[J]. Geophysical Prospecting, 1998, 46(3): 287?301.

[7] 雷林林, 劉四新, 傅磊, 等.基于全波形反演的探地雷達(dá)數(shù)據(jù)逆時(shí)偏移成像[J]. 地球物理學(xué)報(bào), 2015, 58(9): 3346?3355. LEI Linlin, LIU Sixin, FU Lei, et al. Reverse time migration applied to GPR data based on full wave inversion[J]. Chinese Journal of Geophysics, 2015, 58(9): 3346?3355.

[8] 朱尉強(qiáng), 黃清華.探地雷達(dá)衰減補(bǔ)償逆時(shí)偏移成像方法[J]. 地球物理學(xué)報(bào), 2016, 59(10): 3909?3916. ZHU Weiqiang, HUANG Qinghua. Attenuation compensated reverse time migration method of ground penetrating radar signals[J]. Chinese Journal of Geophysics, 2016, 59(10): 3909?3916.

[9] FISHER E, MCMECHAN G A, ANNAN A P, et al. Examples of reverse-time migration of single-channel, ground-penetrating radar profiles[J]. Geophysics, 1992, 57(4): 577?586.

[10] LEUSCHEN C, PLUMB R. A matched-filter approach to wave migration[J]. Journal of Applied Geophysics, 2000, 43(2): 271?280.

[11] ZHOU H, SATO M, LIU H. Migration velocity analysis and prestack migration of common-transmitter GPR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(1): 86?91.

[12] 馮德山, 戴前偉.基于小波多分辨探地雷達(dá)三維正演模擬及偏移處理[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 38(4): 758?763. FENG Deshan, DAI Qianwei. GPR three dimension forward simulation and migration processing with multi-resolution of wavelets[J]. Journal of Central South University(Science and Technology), 2007, 38(4): 758?763.

[13] BRADFORD J H. Reverse-time prestack depth migration of GPR data from topography for amplitude reconstruction in complex environments[J]. Journal of Earth Science, 2015, 26(6): 791?798.

[14] 戴前偉, 張彬, 馮德山, 等. 探地雷達(dá)各向異性介質(zhì)有限差分偏移線性變換[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012, 43(5): 1814?1820. DAI Qianwei, ZHANG Bin, FENG Deshan, et al. Linear transformation of finite difference migration of GPR in anisotropic medium[J]. Journal of Central South University(Science and Technology), 2012, 43(5): 1814?1820.

[15] 傅磊, 劉四新, 劉瀾波, 等.機(jī)載探地雷達(dá)數(shù)值模擬及逆時(shí)偏移成像[J]. 地球物理學(xué)報(bào), 2014, 57(5): 1636?1646. FU Lei, LIU Sixin, LIU Lanbo, et al. Airborne ground penetrating radar numerical simulation and reverse time migration[J]. Chinese Journal Geophysics, 2014, 57(5): 1636?1646.

[16] LIU Sixin, LEI Linlin, FU Lei, et al. Application of pre-stack reverse time migration based on FWI velocity estimation to ground penetrating radar data[J]. Journal of Applied Geophysics, 2014, 107(4): 1?7.

[17] CLAERBOUT J F, JOHNSON A G. Extrapolation of time-dependent waveforms along their Path of Propagation[J]. Geophysical Journal International, 1971, 26(1/2/3/4): 285?293.

[18] 戴前偉, 王洪華, 馮德山, 等.基于雙二次插值的探地雷達(dá)有限元數(shù)值模擬[J]. 地球物理學(xué)進(jìn)展, 2012, 27(2): 736?743. DAI Qianwei, WANG Honghua, FENG Deshan, et al. Finite element numerical simulation for GPR based on quadratic interpolation[J]. Progress in Geophysics, 2012, 27(2): 736?743.

[19] DROSSAERT F H, GIANNOPOULOS A. A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[J]. Geophysics, 2007, 72(2): T9?T17.

[20] 劉紅偉, 劉洪, 鄒振, 等.地震疊前逆時(shí)偏移中的去噪與存儲(chǔ)[J]. 地球物理學(xué)報(bào), 2010, 53(9): 2171?2180. LIU Hongwei, LIU Hong, ZHOU Zhen, et al. The problems of denoise and storage in seismic reverse time migration[J]. Chinese Journal Geophysics, 2010, 53(9): 2171?2180.

[21] 張智, 劉有山, 徐濤, 等. 彈性波逆時(shí)偏移中的穩(wěn)定激發(fā)振幅成像條件[J]. 地球物理學(xué)報(bào), 2013, 56(10): 3523?3533. ZHANG Zhi, LIU Youshan, XU Tao, et al. A stable excitation amplitude imaging condition for reverse time migration in elastic wave equation[J]. Chinese Journal Geophysics, 2013, 56(10): 3523?3533.

[22] SCHLEICHER J, COSTA J, NOVAIS A. A comparison of imaging conditions for wave-equation shot-profile migration[J]. Geophysics, 2008, 73(6): S219?S227.

[23] CHATTOPADHYAY S, MCMECHAN G A. Imaging conditions for prestack reverse-time migration[J]. Geophysics, 2008, 73(3): S81?S89.

(編輯 陳燦華)

Reverse time migration of ground penetrating radar based onnormalized cross correlation imaging condition

CHEN Depeng1, 2, 3, DAI Qianwei1, 2, FENG Deshan1, 2, WANG Honghua3, ZHANG Bin1, 2

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring,Ministry of Education, Changsha 410083, China;3. School of Earth Sciences, Guilin University of Technology, Guilin 541004, China)

On account of the low imaging precision of ground penetrating radar (GPR) reverse time migration (RTM) based on the zero time imaging condition and it is difficult to imaging complex structure with high resolution, normalized cross-correlation imaging conditions was applied in GPR pre-stack RTM imaging. Based on the two-dimensional GPR wave equation, the time domain finite element method based on perfectly matched layer (PML) boundary condition (FETD) was used to simulate the forward and reverse electromagnetic waves, normalized cross-correlation imaging condition was used to obtain the pre-stack RTM results, the spatial high-pass filtering was applied to suppress cross-correlation process of low frequency noise and the corresponding program of GPR pre-stack RTM based on normalized cross-correlation imaging condition was designed and implemented. After that, the two complex GPR models were established and calculated by GPR RTM program based on the normalized cross-correlation imaging condition and zero time imaging condition, respectively. The results show that, compared with imaging results with zero time imaging condition, the imaging results with those in the normalized cross-correlation imaging condition have higher image quality and can more clearly display the detailed spatial form and internal structure information of the abnormal body space.

ground penetrating radar; reverse time migration; normalized cross correlation imaging condition; zero time imaging condition

10.11817/j.issn.1672?7207.2018.05.025

P631

A

1672?7207(2018)05?1221?07

2017?10?10；

2017?12?08

國(guó)家自然科學(xué)基金資助項(xiàng)目(41374118, 41574116, 41604102, 41704128)；廣西自然科學(xué)基金資助項(xiàng)目(2016GXNSFBA380082, 2016GXNSFBA380215)(Projects(41374118, 41574116, 41604102, 41704128) supported by the National Natural Science Foundation of China; Projects(2016GXNSFBA380082, 2016GXNSFBA380215) supported by Guangxi Natural Science Foundation)

戴前偉，博士，教授，從事電磁方法及理論研究；E-mail: qwdai@csu.edu.cn