高陡起伏地形GPR偏移的波場延拓及精度分析

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

高陡起伏地形GPR偏移的波場延拓及精度分析

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

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

探地雷達；偏移成像；高陡起伏地形；交叉微分算子；波場延拓

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

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

逆時偏移是波向外傳播的逆過程，即倒退至零時刻的成像條件，將所有反射與繞射的能量外推至零時刻所處位置的波場[19-20]。對于探地雷達，逆時偏移的定解問題可描述成

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

代入方程(1)中，得

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

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

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

相應地，存在交叉逆變換：

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

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

2　交叉微分算子偏移的波場延拓

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

該式的通解形式為

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

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

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

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

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

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

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

3　數(shù)值算例

3.1高陡起伏地形反射模型

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

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

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

3.2波場延拓精度與解析結果

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

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

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

3.3振幅保真加權值

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

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

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

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

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

3.4高陡界面能量歸位能力

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

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

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

3.5成像剖面

波場延拓精度、振幅保真加權值、高陡界面能量歸位能力最終影響偏移成像質量，4種偏移方法的成像結果見圖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)可見：第1個高陡界面均基本恢復初始形態(tài)，僅頂點處產生較強繞射波并未收斂；第2個高陡界面則成像十分模糊，難以分辨。Kirchhoff法的成像剖面提高了第2個高陡界面的成像質量，基本形態(tài)能夠識辨，僅在繞射界面及繞射頂點周圍存在一部分散亂能量，如圖9(c)所示。圖9(d)所示為交叉微分算子差分法的成像剖面，2個高陡界面形態(tài)清晰，且2個頂點處能量完全收斂，還原了初始模型中高陡起伏的形態(tài)，成像質量較高。

4　結論

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

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

[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]雷林林,劉四新,傅磊,等.基于全波形反演的探地雷達數(shù)據(jù)逆時偏移成像[J].地球物理學報,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]馮德山,戴前偉,何繼善.探地雷達的正演模擬及有限差分波動方程偏移處理[J].中南大學學報(自然科學版),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偏移法在探地雷達正演圖像處理中的應用[J].地球物理學進展,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波動方程偏移的頻率域插值方法[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]傅磊,劉四新,劉瀾波,等.機載探地雷達數(shù)值模擬及逆時偏移成像[J].地球物理學報,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

尹小波，教授，從事淺層地球物理試驗檢測方法及理論研究；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

國家自然科學基金資助項目(41374118，41574116)；中國鐵路總公司科技研究開發(fā)計劃項目(2014G005-B)；湖南省交通科技計劃項目(201423)；湖南省住房和城鄉(xiāng)建設廳科技計劃項目(BZ201408，BZ201411)；中南大學創(chuà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)