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

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

?

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

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

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

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

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

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

1 GPR逆時偏移成像方法

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

1.1 電磁波場正向傳播及逆向外推

考慮二維地電模型的走向為軸，根據(jù)電磁波理論，TM模式下電場分量滿足的波動方程為[18]

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

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

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

1.2 空間高通濾波

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

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

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

2.1 零時刻成像條件

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

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

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

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

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

3 數(shù)值算例

3.1 模型一

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

圖2 GPR模型一示意圖

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

圖3 模型一GPR正演剖面

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

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

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

1—歸一化互相關(guān)成像條件；2—零時刻成像條件。

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

圖6 GPR模型二示意圖

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

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

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

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

1—歸一化互相關(guān)成像條件；2—零時刻成像條件。

4 結(jié)論

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

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

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

[1] 劉瀾波, 錢榮毅. 探地雷達(dá): 淺表地球物理科學(xué)技術(shù)中的重要工具[J]. 地球物理學(xué)報, 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é)進展, 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ù)逆時偏移成像[J]. 地球物理學(xué)報, 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] 朱尉強, 黃清華.探地雷達(dá)衰減補償逆時偏移成像方法[J]. 地球物理學(xué)報, 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é)報(自然科學(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é)報(自然科學(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] 傅磊, 劉四新, 劉瀾波, 等.機載探地雷達(dá)數(shù)值模擬及逆時偏移成像[J]. 地球物理學(xué)報, 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é)進展, 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] 劉紅偉, 劉洪, 鄒振, 等.地震疊前逆時偏移中的去噪與存儲[J]. 地球物理學(xué)報, 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] 張智, 劉有山, 徐濤, 等. 彈性波逆時偏移中的穩(wěn)定激發(fā)振幅成像條件[J]. 地球物理學(xué)報, 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

國家自然科學(xué)基金資助項目(41374118, 41574116, 41604102, 41704128)；廣西自然科學(xué)基金資助項目(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