高陡起伏地形GPR偏移的波場(chǎng)延拓及精度分析

石明，戴前偉，馮德山，張彬，尹小波(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙，410083；2.貴州省有色金屬和核工業(yè)地質(zhì)勘查局物化探總隊(duì)，貴州 都勻，558004；3.有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410083)

高陡起伏地形GPR偏移的波場(chǎng)延拓及精度分析

石明1,2，戴前偉1,3，馮德山1,3，張彬1,3，尹小波1,3
(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙，410083；
2.貴州省有色金屬和核工業(yè)地質(zhì)勘查局物化探總隊(duì)，貴州 都勻，558004；
3.有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410083)

基于高陡起伏地形對(duì)探地雷達(dá)波場(chǎng)的成像精度產(chǎn)生影響，將傳統(tǒng)的單微分算子改進(jìn)為交叉微分算子(CDO);基于交叉微分算子的差分逆時(shí)偏移改變傳統(tǒng)差分偏移算法的波場(chǎng)延拓方向，給出無(wú)近似的波場(chǎng)延拓差分公式、上行波邊界條件和實(shí)現(xiàn)過(guò)程，得到能適應(yīng)于高陡起伏地形的探地雷達(dá)偏移算法，然后以橫向速度差異較大的高陡起伏反射模型為例，從波場(chǎng)延拓精度與解析結(jié)果、振幅保真加權(quán)值、高陡界面能量歸位能力、成像剖面4方面分析交叉微分算子偏移法與傳統(tǒng)方法的成像精度。研究結(jié)果表明：基于交叉微分算子的探地雷達(dá)偏移算法能夠?qū)崿F(xiàn)高陡界面的正確歸位。

探地雷達(dá)；偏移成像；高陡起伏地形；交叉微分算子；波場(chǎng)延拓

探地雷達(dá)(ground penetrating radar,GPR)數(shù)據(jù)偏移成像是淺層地球物理勘探精細(xì)化反演的重要組成部分，反演結(jié)果為偏移成像提供原始的初始速度模型，偏移成像為反演結(jié)果提供有效的校正依據(jù)[1-2]。偏移成像的目的是使界面上的反射波歸位，繞射點(diǎn)收斂[3]。探地雷達(dá)數(shù)據(jù)偏移從自激自收剖面的最后1個(gè)采樣點(diǎn)開(kāi)始計(jì)算，沿著負(fù)時(shí)間方向進(jìn)行延拓，根據(jù)CLAERBOUT[4]的爆炸反射界面成像原理，取零時(shí)刻的深度剖面為偏移成像剖面。復(fù)雜地質(zhì)構(gòu)造的波動(dòng)方程偏移成像方法主要有頻率波數(shù)域F-K偏移、Kirchhoff積分法偏移、差分偏移和有限元偏移。MUFTI等[5-7]采用基于單程標(biāo)量波動(dòng)方程的差分偏移法，該方法可靈活地處理橫向速度劇烈變化的傾斜界面，對(duì)小傾角起伏的復(fù)雜地層，該算法對(duì)原波動(dòng)方程進(jìn)行了相應(yīng)的階層近似，由于略去了波場(chǎng)沿深度方向上的高階偏導(dǎo)，這種近似處理對(duì)偏移效果影響嚴(yán)重，使得該方法受反射界面的傾角限制，同時(shí)解的精度和穩(wěn)定性也容易受空間采樣率、近似的階乘等級(jí)和實(shí)現(xiàn)方法的嚴(yán)重影響。GRAY等[8-11]采用基于射線理論的Kirchhoff積分偏移法，將分散在各道且來(lái)自同一個(gè)繞射點(diǎn)的能量進(jìn)行收攏，該方法計(jì)算效率高，但需要較精確的速度模型，也不易確定對(duì)應(yīng)的偏移孔徑和格林函數(shù)，難以處理橫向速度劇烈變化的高陡界面，對(duì)復(fù)雜地質(zhì)構(gòu)造的空間歸位能力較差。STOLT等[12-15]采用改進(jìn)的頻率波數(shù)域相移偏移法經(jīng)過(guò)2次傅里葉變換，將像空間e( x,z=0,t)變換為E( kx,z=0,w)，再對(duì)E( kx,z=0,w)進(jìn)行波場(chǎng)延拓，外推至深度為z處的波場(chǎng)E( kx,z,w)，最后再進(jìn)行2次一維傅里葉逆變換e( x,z,t=0)。該方法假定波速為恒定值，能夠得到精確解，具有精度高、不受地層起伏限制等優(yōu)點(diǎn)，但該方法均要求地下介質(zhì)的橫向速度不變，難以適應(yīng)橫向速度劇烈變化的高陡起伏地形。由此可見(jiàn)，將各類(lèi)偏移方法應(yīng)用于探地雷達(dá)領(lǐng)域都有其適用范圍和局限性。傳統(tǒng)的差分偏移法求解近似波動(dòng)方程，對(duì)雷達(dá)波場(chǎng)向下進(jìn)行延拓成像，能夠適應(yīng)劇烈的橫向速度變化，但近似的差分偏移方程受地形起伏的影響嚴(yán)重，對(duì)高陡地形反射界面的偏移效果較差，尤其對(duì)高陡且橫向速度劇烈變化的反射界面，常規(guī)時(shí)間偏移方法難以使其中的復(fù)雜構(gòu)造正確歸位。為此，本文作者從傳統(tǒng)差分偏移法中上行波方程產(chǎn)生的誤差根源出發(fā)，將傳統(tǒng)的單微分算子改進(jìn)為交叉微分算子(crossed differential operator,CDO)，對(duì)波動(dòng)方程進(jìn)行交叉坐標(biāo)變換，并采用全新的交叉微分算子差分格式和上行波邊界條件，改變雷達(dá)波場(chǎng)的延拓方向，以地表雷達(dá)波場(chǎng)數(shù)據(jù)作為初始數(shù)據(jù)，根據(jù)逆時(shí)成像原理[16-18]，從地表記錄的波場(chǎng)開(kāi)始，將某個(gè)時(shí)間的時(shí)間切片逆時(shí)傳播進(jìn)入地下，并由波場(chǎng)從切片時(shí)間點(diǎn)的反向進(jìn)行傳播，在時(shí)間軸上實(shí)現(xiàn)外推，恢復(fù)零時(shí)刻各深度上的雷達(dá)波場(chǎng)值，并確保整個(gè)波場(chǎng)延拓的偏移過(guò)程采用非近似的精確公式。最后以廣泛存在于淺地表的高陡起伏地形反射界面為算例，建立高陡起伏的反射模型，重點(diǎn)討論橫向速度劇烈變化情況下高陡界面的雷達(dá)波場(chǎng)延拓，并從波場(chǎng)延拓精度與解析解、振幅保真率、高陡界面能量歸位能力、成像剖面等方面著重分析不同方法的成像精度。

1　傳統(tǒng)差分與交叉微分算子差分法

逆時(shí)偏移是波向外傳播的逆過(guò)程，即倒退至零時(shí)刻的成像條件，將所有反射與繞射的能量外推至零時(shí)刻所處位置的波場(chǎng)[19-20]。對(duì)于探地雷達(dá)，逆時(shí)偏移的定解問(wèn)題可描述成

其中： v( x, z)為介質(zhì)速度，由記錄的雷達(dá)波記錄P( x, z=0,t)逆時(shí)向下延拓，傳統(tǒng)差分法對(duì)方程(1)進(jìn)行浮動(dòng)變換，針對(duì)不同的起伏界面進(jìn)行相應(yīng)的階乘近似：

代入方程(1)中，得

方程(5)未經(jīng)任何近似，系數(shù)α由變換模式和速度參數(shù)決定。假定存在類(lèi)似的交叉變換：

其中：系數(shù)c1，c2，d1和d2為常數(shù)。聯(lián)合方程(1)和(6)得：

對(duì)比方程(5)和方程(7)，得

相應(yīng)地，存在交叉逆變換：

即得交叉微分算子的二階方程：

方程(11)是波動(dòng)方程(1)的精確表達(dá)式，也是交叉微分算子偏移方程的定義式，未進(jìn)行任何階乘近似，可適合高陡的起伏界面偏移成像。

2　交叉微分算子偏移的波場(chǎng)延拓

在一維情況下，方程(1)和(11)均有如下形式：

該式的通解形式為

對(duì)于交叉微分算子偏移方法，在處理雷達(dá)上行波和下行波時(shí)，方程(13)有雙程解，分別為和，且式(13)同時(shí)具有雷達(dá)上行波解和下行波解，在進(jìn)行波場(chǎng)延拓偏移時(shí)，給出邊界條件：

即在上行波邊界條件(14)的約束下，交叉微分算子偏移方程僅有上行波解。

式(15)中右邊第1項(xiàng)為上行波解，第2項(xiàng)為與時(shí)間無(wú)關(guān)的深度函數(shù)()zΨ，通過(guò)略去高階微分算子，通常采用穩(wěn)定的對(duì)稱(chēng)隱式Crank-Nicolson差分格式進(jìn)行延拓。顯然，盡管傳統(tǒng)的差分偏移方法消除了下行波的影響，但是額外增加了1個(gè)與時(shí)間無(wú)關(guān)的深度函數(shù)Ψ(z)。

傳統(tǒng)差分偏移方法與交叉微分算子法的差分網(wǎng)格及波場(chǎng)的延拓方向分別如圖1和圖2所示。

圖1　傳統(tǒng)差分偏移法的波場(chǎng)延拓方向Fig.1　Extrapolation direction of wave field in traditional difference migration algorithms

圖2　交叉微分算子差分偏移法的波場(chǎng)延拓方向Fig.2　Extrapolation direction of wave field in crossed differential operator migration algorithms

得到交叉微分算子偏移方程的隱式 Crank-Nicolson差分格式：

3　數(shù)值算例

3.1高陡起伏地形反射模型

為分析交叉微分算子差分偏移算法對(duì)起伏界面的延拓精度，設(shè)置1個(gè)陡緩程度差異較大的起伏地形反射模型，如圖3所示。其中計(jì)算區(qū)域是1個(gè)長(zhǎng)×寬為3.0 m×1.5 m的矩形，起伏地形的陡緩程度分別為50° 和80°，高陡兩側(cè)的橫向速度差異較大，介質(zhì)的主要介電參數(shù)分別為：相對(duì)介電常數(shù)1ε=5.0，2ε=15.0，電導(dǎo)率1σ=0.4 mS/m，2σ=0.2 mS/m，磁導(dǎo)率1μ=2.0 H/m，2μ=1.0 H/m。正演計(jì)算區(qū)域的邊界采用PML吸收邊界條件[20](4個(gè)方向上的邊界厚度均占據(jù)5個(gè)網(wǎng)格)。

圖3　高陡起伏地形的反射模型Fig.3　Reflection model with steep-dip rough terrains

正演計(jì)算結(jié)果如圖4所示，反射模型中起伏界面的形態(tài)、陡緩程度、水平界面形態(tài)均存在不同程度的畸變，尤其是高陡界面兩側(cè)及高陡界面的頂點(diǎn)繞射現(xiàn)象嚴(yán)重，以致第2個(gè)高陡界面兩側(cè)的水平界面更難以分辨，且反射同相軸的振幅也存在嚴(yán)重縮減現(xiàn)象。

3.2波場(chǎng)延拓精度與解析結(jié)果

波場(chǎng)延拓的外推精度直接影響成像質(zhì)量，相鄰2個(gè)高陡界面的橫向距離較小，使得高陡界面的波場(chǎng)延拓值存在較大誤差。根據(jù)初始模型中各介質(zhì)的主要介電參數(shù)特征，建立相匹配的速度模型，針對(duì)第2個(gè)高陡起伏地形的左側(cè)界面，給出了該界面附近區(qū)域的延拓結(jié)果，即為偏移剖面的第80~90道和0~15 ns區(qū)域內(nèi)不同深度的波場(chǎng)延拓值，結(jié)果如圖5所示。

圖4　高陡起伏地形反射模型的正演結(jié)果Fig.4　Simulation results of reflection model with steep-dip rough terrains

圖5(a)所示為速度模型得到的理論延拓結(jié)果，將延拓值矩陣在橫縱向上的長(zhǎng)度進(jìn)行歸一化，計(jì)算反射同相軸的相對(duì)傾斜度，解析結(jié)果中顯示的同相軸傾斜度與模型中設(shè)置的高陡界面傾角相符。圖5(b)所示為傳統(tǒng)差分偏移法的延拓結(jié)果，同相軸的相對(duì)軸傾斜度約為40°，且在高陡界面的底部，延拓結(jié)果往淺部上移，使得道集的能量并未收斂至高陡界面上，在傾斜度和道集的能量收斂方面均存在較大誤差，該現(xiàn)象在F-K方法所得結(jié)果中更為嚴(yán)重，道集的能量分布更散，如圖5(c)所示，可見(jiàn)第80道波的9~11 ns深度范圍內(nèi)仍存在較大的延拓幅值。圖5(d)所示為Kirchhoff法的延拓結(jié)果，計(jì)算得到的同相軸相對(duì)傾斜度約為44°，并且道集的能量收斂程度也有較大提高，但第83~87道波的9~10 ns深度上出現(xiàn)了“虛假”的反射同相軸。圖5(e)所示為交叉微分算子差分法的延拓結(jié)果，在同相軸的相對(duì)傾斜度和道集的能量收斂方面取得了較好的平衡，且在高陡界面的底部并未出現(xiàn)延拓值明顯較小的趨勢(shì)，延拓精度大幅度提高。

3.3振幅保真加權(quán)值

為了更細(xì)致地比較各種方法的振幅保真能力，計(jì)算模型中第2個(gè)高陡界面頂點(diǎn)(能量繞射聚焦點(diǎn))位置的單道信號(hào)，并重點(diǎn)考察深度上第650~1 050個(gè)樣本的振幅。

圖5　高陡起伏地形的波場(chǎng)延拓與解析結(jié)果的對(duì)比Fig.5　Comparison of wavefield extrapolation in steep-dip rough terrains with analytical results using different methods

圖6　不同方法在第2個(gè)高陡界面繞射點(diǎn)的單道波信號(hào)Fig.6　Trace signal of the second diffracting point in steep-dip rough interface with different methods

該振幅均已進(jìn)行均一化處理，計(jì)算結(jié)果見(jiàn)圖6。從圖6可見(jiàn)：在第650~1 050個(gè)樣本深度范圍內(nèi)，傳統(tǒng)FD法和F-K法在信號(hào)振幅、信號(hào)位置方面取得的效果相當(dāng)，Kirchhoff法則在第550~750個(gè)樣本深度上的信號(hào)位置出現(xiàn)了偏差，信號(hào)振幅保真卻比前2種略有提高，而CDO差分方法的最大振幅則約為其他3種方法的2倍，表明振幅得到了較好保真。

將延拓精度對(duì)應(yīng)的解析夾角、偏移夾角和相對(duì)振幅值建立振幅保真的加權(quán)函數(shù)[21]：

3.4高陡界面能量歸位能力

為了更直觀地比較各類(lèi)偏移方法對(duì)高陡界面的能量收斂能力，將提取的各單道信號(hào)能量進(jìn)行疊加堆積，得到不同偏移方法所得的剖面能量圖譜，如圖8所示。能量圖譜顯示了不同偏移方法對(duì)計(jì)算區(qū)域內(nèi)繞射波的收斂能力，尤其體現(xiàn)了對(duì)高陡界面反射能量的歸位水平。

圖8(a)所示為傳統(tǒng)差分偏移法的能量分布結(jié)果。從8(a)可見(jiàn)：受高陡地形的影響，第1個(gè)界面能量并未完全收斂，反射界面能量軸的附近存在能量損失，界面附近出現(xiàn)能量重疊的影子。該“重疊現(xiàn)象”同樣出現(xiàn)在F-K偏移法所得結(jié)果中，并且在模型頂部的兩端存在較明顯的剩余能量，如圖8(b)所示。圖8(c)所示為Kirchhoff偏移法的能量分布結(jié)果，2個(gè)高陡界面的能量基本收斂，但形態(tài)模糊，界面上的能量分布較分散，表明在較精確的速度模型約束下，Kirchhoff偏移法對(duì)橫向變速介質(zhì)中的高陡界面具有較強(qiáng)的歸位能力。圖8(d)所示為交叉微分算子差分偏移法的能量分布結(jié)果，顯然，界面上的能量均已集中，繞射頂點(diǎn)的能量也基本收斂，且計(jì)算區(qū)域內(nèi)并未見(jiàn)剩余能量，高陡界面的能量歸位程度較高。

圖7　不同偏移方法在不同傾角反射界面的歸一化振幅曲線Fig.7　Normalized amplitude along different angle reflection interface with different migration method

3.5成像剖面

波場(chǎng)延拓精度、振幅保真加權(quán)值、高陡界面能量歸位能力最終影響偏移成像質(zhì)量，4種偏移方法的成像結(jié)果見(jiàn)圖9。

圖8　不同偏移方法的能量圖譜剖面Fig.8　Energy spectrum section with different migration methods

圖9　不同方法的偏移成像剖面Fig.9　Migration imaging section with different method

圖9(a)和圖9(b)所示分別為傳統(tǒng)差分法和F-K法的偏移成像剖面。從圖9(a)和圖9(b)可見(jiàn)：第1個(gè)高陡界面均基本恢復(fù)初始形態(tài)，僅頂點(diǎn)處產(chǎn)生較強(qiáng)繞射波并未收斂；第2個(gè)高陡界面則成像十分模糊，難以分辨。Kirchhoff法的成像剖面提高了第2個(gè)高陡界面的成像質(zhì)量，基本形態(tài)能夠識(shí)辨，僅在繞射界面及繞射頂點(diǎn)周?chē)嬖谝徊糠稚y能量，如圖9(c)所示。圖9(d)所示為交叉微分算子差分法的成像剖面，2個(gè)高陡界面形態(tài)清晰，且2個(gè)頂點(diǎn)處能量完全收斂，還原了初始模型中高陡起伏的形態(tài)，成像質(zhì)量較高。

4　結(jié)論

1)為提高探地雷達(dá)對(duì)高陡起伏界面的偏移歸位能力，將傳統(tǒng)的單微分算子改進(jìn)為交叉微分算子，提出基于交叉微分算子的差分偏移方法。該方法采用非近似的公式進(jìn)行波場(chǎng)延拓，改變雷達(dá)波場(chǎng)的延拓方向，并給出交叉微分算子偏移法的波場(chǎng)延拓過(guò)程和相應(yīng)的差分格式，避免了傳統(tǒng)差分偏移法中略去高階偏導(dǎo)。

2)重點(diǎn)討論了橫向變速介質(zhì)中高陡起伏地形的雷達(dá)波波場(chǎng)延拓算例，并從波場(chǎng)延拓精度與解析解、成像的振幅保真加權(quán)值、高陡界面能量歸位能力和成像剖面4個(gè)方面，比較了傳統(tǒng)差分法、F-K頻率波數(shù)域法、Kirchhoff積分法和交叉微分算子差分法的偏移成像效果，分析了各類(lèi)方法的成像精度。結(jié)果表明交叉微分算子差分偏移法優(yōu)于傳統(tǒng)的差分偏移法，提高了雷達(dá)波場(chǎng)的延拓精度，且振幅得到較好地保真，更有效地指導(dǎo)探地雷達(dá)高陡起伏地形的正確歸位。

[1]FENG X,SATO M,LIU C,et al.Topographic correction of elevated GPR[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,7(3):799-804.

[2]雷林林,劉四新,傅磊,等.基于全波形反演的探地雷達(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.

[3]CHARLES P O,MICHAEL H P,DAVID L W,et al.Improving GPRimageresolutioninlossygroundusingdispersive migration[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45(8):2492-2500.

[4]CLAERBOUT J F.Imaging the earth’s interior[M].Boston: Blackwell Scientific Publications,1985:78-92.

[5]MUFTI I R,PITA J A,HUNTLEY R W.Finite-difference depth migration of exploration-scale 3-D seismic data[J].Geophysics, 1996,61(3):776-794.

[6]馮德山,戴前偉,何繼善.探地雷達(dá)的正演模擬及有限差分波動(dòng)方程偏移處理[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2006, 37(2):361-365. FENG Deshan,DAI Qianwei,HE Jishan.Forward simulation of ground penetrating radar and its finite difference method waveequation migration processing[J].Journal of Central South University(Science and Technology),2006,37(2):361-365.

[7]COSTA C,CAMPOS I,COSTA J C,et al.Iterative methods for 3D implicit finite-difference migration using the complex Padé approximation[J].Journal of Geophysics and Engineering,2013, 10(4):45011-45027.

[8]GRAY S H,MAY W P.Kirchhoff migration using eikonal equation travel time[J].Geophysics,1994,59(5):810-817.

[9]戴前偉,馮德山,何繼善.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 penetratingradar forward simulate[J].Progressin Geophysics,2005,20(3):849-853.

[10]CHENG C.Three-dimensional pre-stack depth migration of receiver functions with the fast marching method:a Kirchhoff approach[J].GeophysicalJournal International,2016,205: 819-829.

[11]LIN L C,SHI B P,AN P.Multiwavelet prestack Kirchhoff migration[J].Geophysics,2016,81(3):s79-s85.

[12]STOLT R H.Migration by Fourier transforms[J].Geophysics, 1978,43(1):23-48.

[13]LOEWENTHAL D,MUFTI I R.Reversed time migration in spatial frequency domain[J].Geophysics,1983,48(5):627-635.

[14]GAZDAG J,SGUAZZERO P.Migration of seismic data by phase shift plus interpolation[J].Geophysics,1994,49(1): 124-131.

[15]賀振華,GARDNER G H.F-K波動(dòng)方程偏移的頻率域插值方法[J].石油物探,1985,24(1):1-22. HE Zhenhua,GARDNER G H.Interpolation in the Fourier transformdomainwithapplicationstoF-Kmigration[J]. Geophysical Prospecting for Petroleum,1985,24(1):1-22.

[16]SUN R,MCMECHAN G A.Scalar reverse-time depth migration of prestack elastic seismic data[J].Geophysics,2001,66(5): 1519-1527.

[17]DAVYDENKOM,VERSCHUURDJ.Full-wavefield migration:using surface and internal multiples in imaging[J]. Geophysical Prospecting,2016,64(2):505-513.

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

[19]BAYSAL E,KOSLOFF D D,SHERWOOD J W C.Reversetime migration[J].Geophysics,1983,48(11):1514-1524.

[20]BERENGER J P.A perfectly matched layer for the absorption of electromagnetic-waves[J].Journal of Computational Physics, 1994,114:185-200.

[21]EKREN B O,URSIN B.True-amplitude frequency-wavenumber constant offset migration[J].Geophysics,1999,64:915-924.

(編輯陳燦華)

Precision analysis of wave-field extrapolation in migration with steep-dip rough terrains for GPR

SHI Ming1,2,DAI Qianwei1,3,FENG Deshan1,3,ZHANG Bin1,3,YIN Xiaobo1,3
(1.School of Geosciences and Info-Physics,Central South University,Changsha 410083,China;
2.Geophysical and Geochemical Prospecting Team,Nonferrous Metals and Nuclear Industry Geological Exploration Bureau of Guizhou,Duyun 558004,China;
3.Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring, Ministry of Education,Changsha 410083,China)

Considering that ground penetrating radar(GPR)is obviously influenced by the inclination angle of reflecting interface,in order to improve the imaging accuracy of steep-dip rough terrains in GPR migration,the single differential operator was replaced by a cross differential operator(CDO),and a reverse-time migration method based on CDO was proposed,which completely changed the direction of wave-field extrapolation,and the implementation of CDO FD migration method was eventually achieved with the precise difference equations and corresponding boundary conditions of up-going GPR wave.Then a laterally variable velocity reflection model with steep-dip rough terrains was presented, particular attention was paid to the precision of wave-field extrapolation in laterally variable velocity,and from wave-field extrapolation precision,fidelity weight of wave amplitude,ability of energy convergence in steep-dip case and migration imaging section to analyze the imaging precision obtained by CDO FD migration method,FD method,F-K method and Kirchhoff method.The results show that the diffraction and scattering can be mostly converged with COD FD migration method,which shows that the COD FD migration method can deal with steep-dip rough terrains of GPR data.

ground penetrating radar(GPR);migration imaging;steep-dip rough terrains;crossed differential operator (CDO);wave field extrapolation

尹小波，教授，從事淺層地球物理試驗(yàn)檢測(cè)方法及理論研究；E-mail:yinxiaobo@csu.edu.cn

P631

A

1672-7207(2016)07-2448-08

10.11817/j.issn.1672-7207.2016.07.036

2016-03-14；

2016-05-15

國(guó)家自然科學(xué)基金資助項(xiàng)目(41374118，41574116)；中國(guó)鐵路總公司科技研究開(kāi)發(fā)計(jì)劃項(xiàng)目(2014G005-B)；湖南省交通科技計(jì)劃項(xiàng)目(201423)；湖南省住房和城鄉(xiāng)建設(shè)廳科技計(jì)劃項(xiàng)目(BZ201408，BZ201411)；中南大學(xué)創(chuàng)新驅(qū)動(dòng)項(xiàng)目(2015CX008) (Projects(41374118,41574116)supported by the National Natural Science Foundation of China;Project(2014G005-B)supported by the Key Program of Science and Technology of China Railway Corporation;Project(201423)supported by the Traffic Department of Science and Technology Program of Hunan Province;Projects(BZ201408,BZ201411)supported by the Science and Technology Plan of Housing and Construction Department of Hunan Province;Project(2015CX008)supported by Innovation-Driven Plan of Central South University)