雷林林, 劉四新, 傅磊*, 吳俊軍

1 吉林大學地球探測科學與技術(shù)學院, 長春 130026 2 中國石油集團東方地球物理勘探有限責任公司 新興物探開發(fā)處， 河北涿州 072751

?

基于全波形反演的探地雷達數(shù)據(jù)逆時偏移成像

雷林林1, 劉四新1, 傅磊1*, 吳俊軍2

1 吉林大學地球探測科學與技術(shù)學院, 長春 130026 2 中國石油集團東方地球物理勘探有限責任公司 新興物探開發(fā)處， 河北涿州 072751

逆時偏移成像(RTM)常用來處理復雜速度模型，包括陡傾角及橫向速度變化劇烈的模型.與常規(guī)偏移成像方法(如Kirchhoff偏移)相比，逆時偏移成像能提供更好的偏移成像結(jié)果，近些年逆時偏移成像越來越廣泛地應(yīng)用到勘探地震中，它逐漸成為石油地震勘探中的一種行業(yè)標準.電磁波和彈性波在動力學和運動學上存在相似性，故本文開發(fā)了基于麥克斯韋方程組的電磁波逆時偏移成像算法，并將其應(yīng)用到探地雷達數(shù)據(jù)處理中.時間域有限差分(FDTD)用于模擬電磁波正向和逆向傳播過程，互相關(guān)成像條件用于獲得最終偏移結(jié)果.逆時偏移成像算法中，偏移成像結(jié)果受初始模型影響較大，而其中決定電磁波傳播速度的介電常數(shù)的影響尤為重要.本文基于時間域全波形反演(FWI)算法反演獲得了更為精確的地下介電常數(shù)模型，并將其反演結(jié)果作為逆時偏移成像的初始介電常數(shù)模型.為了驗證此算法的有效性，首先構(gòu)建了一個復雜地質(zhì)結(jié)構(gòu)模型，合成了共偏移距及共炮點探地雷達數(shù)據(jù)，分別應(yīng)用常規(guī)Kirchhoff偏移算法及逆時偏移成像算法進行偏移處理，成像結(jié)果顯示由逆時偏移成像算法得到的偏移結(jié)果與實際模型具有較高的一致性；此外本文在室內(nèi)沙槽中進行了相關(guān)的物理模擬實驗，采集了共偏移距及共炮點探地雷達數(shù)據(jù)，分別應(yīng)用Kirchhoff和疊前逆時偏移成像算法進行處理，結(jié)果表明疊前逆時偏移成像在實際應(yīng)用中能獲得更好的成像效果.

逆時偏移； 探地雷達； 全波形反演； 物理模擬實驗

1 引言

偏移是一種地震成像技術(shù)，它在空間域或時間域?qū)⒂涗浀降耐蜉S歸位到真實的反射界面.傳統(tǒng)的地震偏移成像方法基于單程波方程，并將多次散射波當做噪聲處理，但對于復雜地質(zhì)結(jié)構(gòu)及大傾度界面如裂縫，合理的利用多次散射信息至關(guān)重要.傳統(tǒng)波動方程偏移方法如Kirchhoff積分偏移(Gray, 1986; Gray and May, 1994; Docherty)難以適應(yīng)橫向速度變化；有限差分波動方程偏移(Vidale, 1988; Mufti et al., 1996)近似求解波動方程，存在傾角限制；頻率波數(shù)域偏移方法(Loewenthal and Mufti, 1983)是基于常速度模型，當?shù)叵陆橘|(zhì)速度存在變化時會出現(xiàn)明顯的精度問題.有別于傳統(tǒng)的基于單程波方程的偏移方法，逆時偏移(Baysal et al., 1983; Levin, 1984; Chang and McMechan, 1986，1994; Sun and McMechanz, 2001; Yoon et al., 2003; 劉欣欣等, 2013; 王保利等, 2012)基于全波方程，能夠精確的處理沿任意方向傳播的波場，包括耦合波及橫縱波之間的轉(zhuǎn)換波，從而可更精準地將波場歸位到其真實的地下空間位置.雖然逆時偏移技術(shù)的成功應(yīng)用實例并不多見，但它在勘探地震中受到越來越多的關(guān)注，逆時偏移包括疊前逆時偏移和疊后逆時偏移.

探地雷達電磁波與地震波在運動學和動力學存在一定的相似性，此外，采集到的探地雷達數(shù)據(jù)和地震數(shù)據(jù)都是空間位置的時域電壓信號，故而在一定程度上勘探地震中的數(shù)據(jù)處理技術(shù)可以運用到探地雷達數(shù)據(jù)處理中.Fisher(1992)提出用基于張量波動方程的逆時偏移算法對單通道的共偏移距的探地雷達數(shù)據(jù)進行了偏移處理，處理后的結(jié)果在一定程度顯示了該算法的優(yōu)越性.Feng和Sato(2004)利用疊前逆時偏移算法，對地雷探測中采集到的多偏移距數(shù)據(jù)進行處理，一方面強散射波得到壓制，另一方面獲得了地下更為精細的雷達剖面圖像.Zhou等(2005)應(yīng)用一種基于遺傳算法的偏移速度分析方法，對合成及野外實測的共炮點雷達數(shù)據(jù)進行分析獲得了地下介電常數(shù)分布，并以其作為疊前逆時偏移的初始速度模型，取得了較好的結(jié)果.

本文開發(fā)了電磁波逆時偏移算法，電磁波場的正向及逆向傳播是基于時間域有限差分方法，單軸各向異性完全匹配層(UPML)(Fang and Ying, 2009; Feng and Dai, 2011)吸收邊界條件用于吸收邊界的反射波，互相關(guān)成像條件用于獲得空間圖像.應(yīng)用電磁波全波形反演算法，對合成及實測的共偏移距、多偏移距探地雷達數(shù)據(jù)全波形反演，獲得了與實際情況更為接近的地下介質(zhì)分布情況.通過全波形反演獲得的介電常數(shù)和電導率模型作為逆時偏移的初始模型，逆時偏移的結(jié)果與基于共偏移距雷達數(shù)據(jù)Kirchhoff偏移成像結(jié)果對比發(fā)現(xiàn)，逆時偏移成像效果明顯好于傳統(tǒng)偏移成像結(jié)果.

2 逆時偏移基本原理及電磁波逆時偏移算法

傳統(tǒng)的偏移算法基本是基于單程波方程深度外推，所以在對多分量地震資料處理過程中首先要進行縱、橫波分離，然后用現(xiàn)有的偏移算法分別處理縱波場和橫波場.它的缺陷在于：1)將多分量地震資料當作幾個標量波的簡單疊加，忽略了地震波的矢量特性，從而影響地震多分量資料的處理與解釋精度；2)假如地震波場橫波與縱波沒有完全分離(目前這種現(xiàn)象廣泛存在)，而在實際資料的偏移處理過程中又假設(shè)為標量波場，這必然會使偏移結(jié)果產(chǎn)生許多人為假象，影響處理精度.逆時偏移與傳統(tǒng)的偏移方法有很大的不同，逆時偏移可以看成是波場在時間軸上逆向傳播過程，當波場倒退到零時刻，則所有反射波與繞射波的能量都回到最初被反射和繞射的空間位置，應(yīng)用成像條件即可得到最終的偏移成像剖面.逆時偏移的波動理論方法基于全波方程，波場可以沿各個方向傳播，因此該方法不受介質(zhì)傾角限制，適用于速度任意變化的模型.

逆時偏移算法的實現(xiàn)概括起來主要包括三個步驟：第一是震源波場正向延拓；第二是接收波場的逆時延拓；第三是應(yīng)用成像條件成像.前兩個步驟的實現(xiàn)可以基于相同的數(shù)值方法，本文應(yīng)用時間域有限差分法(劉四新和曾昭發(fā), 2007; 劉四新等, 2006; Taflove and Hagness, 2005; Liu and Sato, 2005; Liu et al., 2007; Liu et al., 2011)實現(xiàn)電磁波場在空間的正向和逆向傳播.影響疊前逆時偏移成像效果的一個非常重要因素是成像條件(Chattopadhyay and McMechan, 2008)的應(yīng)用，Claerbout(1971)提出的互相關(guān)成像條件是最基礎(chǔ)的成像條件，它采用零延遲互相關(guān)成像條件，該成像條件的不足之處在于對首波及背景散射波做互相關(guān)會產(chǎn)生低頻噪聲，其次應(yīng)用該成像條件得到的成像結(jié)果與地下介質(zhì)的反射系數(shù)沒有對應(yīng)關(guān)系.鑒于零延遲互相關(guān)成像條件的弊端，后續(xù)逐步發(fā)展了反褶積成像條件(Guitton et al., 2007)、震源補償成像條件(Guitton et al., 2006)、擴展成像條件(Sava and Fomel, 2006)以及角度域成像條件(Yoon and Marfurt, 2006)等.標準的逆時偏移采用零延遲互相關(guān)成像條件，其計算公式為：

(1)

其中IMAG(x)為地下空間x的成像結(jié)果，S(x,t,si)為源波場，R(x,T-t,si)為接收波場，t和si分別代表時間和源序號.當接收波場傳播到時間t時，讀取震源波場t時刻的場值和對應(yīng)時刻接收波場的場值并將兩波場值相乘，最后將計算區(qū)域?qū)?yīng)網(wǎng)格點相加即可得到最終成像結(jié)果.

本文開發(fā)的電磁波逆時偏移成像算法如圖1所示，對于合成數(shù)據(jù)，首先給定真實模型，通過FDTD模擬電磁波在模型中的傳播過程，同時記錄接收波場.逆時偏移成像的第一步是通過FDTD進行激勵源波場的正向傳播模擬，使用雷克子波作為激勵源函數(shù)，電磁波最大傳播時間序列為nmax，為了節(jié)約存儲空間，將每個時刻模型空間的外圍波場存儲；第二步是電磁波波場的逆向傳播模擬，該步驟包括兩個方面的反向傳播過程：(一)利用包圍模型空間外圍全時序波場恢復電磁波場的正向傳播過程，這樣做的目的其實就是以時間換空間；(二)利用接收波場進行電磁波場的逆向傳播模擬；第三步是成像條件的應(yīng)用，因為互相關(guān)條件容易實現(xiàn)且無穩(wěn)定性問題，本文選取互相關(guān)作為成像條件.波場逆向傳播的過程中，利用成像條件獲得每個時刻的偏移圖像，當波場傳播到時間零點，將各個時刻的圖像疊加即可獲得單源逆時偏移成像結(jié)果.此外，通過去噪技術(shù)去除低頻影響，再將所有單源偏移結(jié)果疊加得到最終的逆時偏移成像結(jié)果.

圖1 電磁波逆時偏移成像算法流程Fig.1 RTM algorithm flow for electromagnetic wave

3 數(shù)值模擬實驗

圖2 相對介電常數(shù)模型Fig.2 Relative dielectric constant model

文中構(gòu)建了一個復雜地質(zhì)結(jié)構(gòu)模型如圖2所示，X方向長度為1.2 m，Z方向長度為0.7 m.圖中藍色區(qū)域代表空氣，暗紅色區(qū)域代表土壤層，其相對節(jié)電常數(shù)為6，電導率為5 mS·m-1，平均厚度約為0.12 m，該土壤層的底界面起伏變化；黃色區(qū)域代表干沙層，其相對介電常數(shù)為4，電導率為1 mS·m-1，平均厚度約為0.48 m.由于土壤層與干沙層之間存在大傾角起伏界面，故而該模型中存在變化劇烈的橫向速度.干沙中埋有2個異常體，其中藍色方形區(qū)域代表一個高速異常體，其相對介電常數(shù)為1；圓心位于(0.8, 0.4)m處圓形暗紅色區(qū)域代表一個低速異常體，其相對介電常數(shù)為6，電導率為4 mS·m-1.

3.1 雷達合成數(shù)據(jù)

典型的雷達反射測量中，電場與測線往往是垂直的，因此我們應(yīng)用二維TM模FDTD求解Maxwell方程組，同時采用UPML吸收邊界條件來減少邊界處反射.為確保穩(wěn)定性及減少數(shù)值頻散，設(shè)置網(wǎng)格大小為0.004 m，發(fā)射源為主頻1600 MHz的雷克子波，時間取樣間隔0.008 ns，采樣點為1024.首先合成了共偏移距剖面探地雷達數(shù)據(jù)，第一道信號發(fā)射天線的水平距離位于0.05 m，接收天線與發(fā)射天線間距為0.05 m，發(fā)射天線的移動步長為0.012 m，沿地面總該采集了88道雷達信號，水平方向上覆蓋了1.05 m的水平距離.圖3a是合成的共偏移距雷達剖面，橫坐標代表道號，縱坐標代表雷達波的雙程走時，從圖可知在2.5～4.5 ns之間的雙曲線同向軸是來自土壤層和沙層界面的反射波，值得注意的是出現(xiàn)了來自不規(guī)則界面拐點處的若干散射波，這些散射雜波會對地下目標的識別造成困擾.水平位置位于第20道，頂點約在6 ns出現(xiàn)的雙曲線同向軸來自方形高速異常體的反射波，而水平位置位于第60道，頂點約為4.8 ns處出現(xiàn)的雙曲線同向軸則是來自圓形低速異常體的反射波.合成的共炮點雷達數(shù)據(jù)如圖3b所示，水平軸代表道號，縱軸是雷達波的雙程走時.共炮點雷達數(shù)據(jù)集與共偏移距雷達剖面在某種程度上存在一定的相似性，來自界面及異常的反射波出現(xiàn)的雙程走時相當，此外反射波的曲線形態(tài)也存在一定的相似性.

3.2 基于全波形反演的初始模型建立

逆時偏移成像算法中，影響最終成像效果的一個最重要因素是初始模型.為了獲得一個較高精確度的初始模型，本文采用時間域全波形反演技術(shù)，反演地下介質(zhì)的相對介電常數(shù)分布，從而為逆時偏移成像提供更為精準的初始模型.全波形反演技術(shù)試圖利用包括走時、幅度、相位等全波形信息來反演地下參數(shù)分布，它比常規(guī)的基于射線的成像方法能提供更加可靠的參數(shù)信息，本文應(yīng)用的時間域全波形反演技術(shù)細節(jié)可參考(Meles et al., 2010, 2011)，其主要目的是建立包含觀測數(shù)據(jù)和模擬數(shù)據(jù)的目標函數(shù)(2)式：

(2)

時間域全波形反演算法中，不管是單參數(shù)分別反演還是介電常數(shù)和電導率聯(lián)合反演，初始介電常數(shù)和電導率模型的選取至關(guān)重要.考慮到如圖2所示的真實相電常數(shù)模型中存在起伏較大的界面，地下介電常數(shù)模型不易給定，故而這里選取的初始介電常數(shù)模型如圖4a所示，藍色區(qū)域代表空氣；黃色區(qū)域代表地下介質(zhì)，其相對介電常數(shù)為4.5.圖4b是時間域全波形反演25次迭代反演后得到的相對介電常數(shù)結(jié)果，首先土壤與沙層之間的起伏界面得到了清晰的刻畫，起伏界面的形狀及位置與真實模型匹配較好，另外在反演的界面上方相對介電常數(shù)要高于界面下方的相對介電常數(shù)；反演圖像中埋于沙層中的兩個異常體得到了較好的刻畫，兩個異常體的水平位置與實際模型相對應(yīng).右側(cè)低速異常體的位置與真實模型一致，而且其內(nèi)部相對介電常數(shù)值高于周圍背景值；左側(cè)高速異常體的位置與真實模型一致，其內(nèi)部相對介電常數(shù)值低于背景值.

圖5給出了真實模型、初始模型及反演結(jié)果之間的詳細數(shù)值比較，黑色實線代表真實模型的相對介電常數(shù)值，藍色虛線代表初始模型相對介電常數(shù)值，紅色實線是時間域全波形反演迭代25次后的反演結(jié)果.圖5a是經(jīng)過左側(cè)圓形異常體A1-A2線(見圖2)的相對介電常數(shù)比較，橫軸為相對介電常數(shù)值，縱軸對應(yīng)的空間位置.小于0.1m時，反應(yīng)在真實模型中是相對介電常數(shù)為1空氣，可以發(fā)現(xiàn)三條曲線完全重合.0.1～0.3m之間，給出的初始模型與真實模型之間存在1.5的差異，而全波形反演結(jié)果基本介于初始模型與真實模型之間.在0.3m處是界面突變點，反演結(jié)果中在同樣位置出現(xiàn)了拐點.0.5m附近是異常體所在的位置，反演結(jié)果的數(shù)值大小及拐點位置與真實模型基本一致.圖5b是經(jīng)過右側(cè)圓形異常體B1-B2線(見圖2)的相對介電常數(shù)比較，從中可以得到與左圖相似的結(jié)論：全波反演的結(jié)果較為準確的刻畫了地下相對介電常數(shù)分布情況，可以為逆時偏移提供更為精準的初始模型.

圖4 (a)初始介電常數(shù)模型； (b)全波形反演25次迭代后的介電常數(shù)反演結(jié)果Fig.4 (a) The initial permittivity model for the FWI inversion； (b) the inversion permittivity results after 25 iterations

圖5 真實模型，初始模型及時間域全波反演結(jié)果在不同位置的介電常數(shù)比較(a) A1-A2線(見圖2)介電常數(shù)比較; (b) B1-B2線(見圖2)介電常數(shù)比較.Fig.5 The permittivity comparison of true model, initial model and FWI inversion results at different position(a) Permittivity comparison of A1-A2 line (the left black dash line indicated in Fig.2); (b) Permittivity comparison of B1-B2 line (the right black dash line indicated).

3.3 合成數(shù)據(jù)逆時偏移成像分析

為了評價逆時偏移算法應(yīng)用到探地雷達數(shù)據(jù)處理中的效果，選取Kirchhoff偏移算法與之比較.對于逆時偏移成像算法，我們可以將其處理共偏移距數(shù)據(jù)，亦可以處理共炮點數(shù)據(jù).

對于如圖3a所示的共偏移距雷達數(shù)據(jù)，基于真速度模型，用Kirchhoff偏移成像算法獲得的偏移成像結(jié)果圖6a所示.從圖可知，來自不規(guī)則界面的反射波得到一定程度的抑制，可以大致識別出土壤與沙層的界面，然而土壤與沙層界面的散射波沒有完全收斂，仍然存在較嚴重的雜波干擾.偏移剖面中，來自兩個異常體的反射波能量基本上得到了壓制.左側(cè)高速異常體水平位置與真實模型一致，但是垂直位置與實際模型存在一定的偏差；右側(cè)低速異常體的位置與實際模型基本一致.此外需要注意的是成像結(jié)果中出現(xiàn)了一些干擾雜波，這些雜波的存在將使得雷達數(shù)據(jù)的解釋工作更加復雜.因此，對于存在橫向速度變化較大的模型，即便使用真實速度模型進行Kirchhoff偏移，也得不到良好的偏移結(jié)果.

圖6 不同偏移方法應(yīng)用不同初始模型得到的成像結(jié)果(a) 基于共偏移距數(shù)據(jù)的Kirchhoff偏移結(jié)果(真介電常數(shù)模型作為初始模型); (b)基于共偏移距數(shù)據(jù)的RTM結(jié)果(真介電常數(shù)模型作為初始模型); (c) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(地下相對介電常數(shù)取4作為初始模型); (d) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(地下相對介電常數(shù)取6作為初始模型); (e) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(FWI反演結(jié)果作為初始模型); (f) 基于共炮點雷達數(shù)據(jù)集的RTM結(jié)果(真介電常數(shù)模型作為初始模型). 所有圖像都進行了振幅歸一化.Fig.6 Imaging results from different migration methods with different initial model(a) The Kirchhoff migration result based on common offset radar data (true dielectric constant model as the initial model); (b) RTM result based on common offset radar data (true dielectric constant model as the initial model); (c) RTM result based on common offset radar data (dielectric constant of 4 in the underground as the initial model); (d) RTM result based on common offset radar data (dielectric constant of 6 in the underground as the initial model); (e) RTM result based on common offset radar data (FWI result as the initial model); (f) RTM result based on common shot radar data (true dielectric constant model as the initial model). All figures are normalized by their maximum amplitude.

同樣，對于如圖3a所示的共偏移距雷達數(shù)據(jù)，將真介電常數(shù)模型作為初始模型參數(shù)，忽略電導率的影響，用逆時偏移成像獲得的成像結(jié)果如圖6b所示，可以發(fā)現(xiàn)來自起伏界面的散射波已被完全壓制住了，土壤與沙層間的界面可以被清晰的識別出來，成像后的界面形狀及位置與真實模型完全一致.來自兩個異常體的反射波聚焦較好，成像結(jié)果顯示兩個異常體的形態(tài)大小，空間位置與真實模型完全一致，此外偏移剖面中存在的雜波干擾較小.總體來說，逆時偏移成像效果明顯好于Kirchhoff偏移結(jié)果.

逆時偏移成像算法中，初始模型的選擇至關(guān)重要，初始模型的建立可以結(jié)合已知的先驗信息、CMP速度分析、或者反射層析成像結(jié)果.為了研究初始模型對成像精度的影響，文中分別選取地下相對介電常數(shù)為4的均勻半空間，地下相對介電常數(shù)為6的均勻半空間以及選取全波形反演的結(jié)果作為初始速度模型.

對于共偏移距探地雷達數(shù)據(jù)，當選取地下相對介電常數(shù)為4的均勻半空間作為初始介電常數(shù)模型時，采用逆時偏移技術(shù)獲得的成像結(jié)果如圖6c所示，從圖可知來自起伏界面的繞射波沒有完全歸位，雜波干擾嚴重，來自地下異常體的散射波能量聚焦效果較差.當選取地下相對介電常數(shù)為6的均勻半空間作為初始介電常數(shù)模型時，逆時偏移成像結(jié)果如圖6d所示，來自起伏界面的反射波能量聚焦較好，能清晰的識別起伏界面的位置及形態(tài)，但是來自地下異常體的散射波能量沒有完全聚焦.當選取全波形反演的結(jié)果作為初始介電常數(shù)模型時，逆時偏移成像的結(jié)果如圖6e所示，從圖可知，來自界面的反射波能量聚焦較好，能夠清晰的識別起伏界面的形態(tài)和位置；來自地下的兩個異常體的散射波能量聚焦較好，其空間位置與模型基本一致.由此可知，當初始模型越接近于真實模型時，逆時偏移成像結(jié)果越準確；全波形反演能夠為逆時偏移成像提供更為準確的初始模型，將全波形反演的結(jié)果作為逆時偏移成像的初始模型是行之有效的.

對于如圖3b所示的共炮點探地雷達數(shù)據(jù)，采用逆時偏移成像算法對其進行偏移成像，其結(jié)果如圖6d所示，圖中來自起伏界面的反射波能量得到聚焦，起伏界面的空間位置及形態(tài)與實際模型基本一致；來自地下異常體的散射波能量聚焦較好，其空間位置與模型一致.偏移剖面中，在地面處出現(xiàn)了若干個和源位置一致的強能量點；其反射波子波略寬于如圖6b中的反射波子波寬度.

4 物理實驗及數(shù)據(jù)分析

4.1 物理實驗?zāi)Ｐ图皵?shù)據(jù)采集

本文在實驗室沙槽進行了物理模擬實驗，如圖7a所示，沙槽長8 m，寬1.6 m，深0.8 m，沙槽四壁為混泥土，一個PVC盒以45°傾斜角埋于沙中，其埋深約0.24 m.圖7c展示的是PVC盒埋于干沙中的示意圖，PVC盒厚度為0.04 m，其中間為空氣，頂端位于(0.75, 0.24)m處，干沙相對介電常數(shù)是4，PVC的相對介電常數(shù)是3.1.圖7b展示了本實驗中采用的 MALA ProEx雷達控制單元及兩套中心頻率為1600 MHz的高頻天線.當采集共炮點雷達數(shù)據(jù)時，雷達主機外接兩個高頻模塊，同時將兩套高頻天線分別接于高頻模塊；當采集共偏移距雷達數(shù)據(jù)時，雷達主機外接一個高頻模塊，同時只需要一個高頻天線.

首先用1600 MHz高頻天線進行了共偏移距測量，測量過程中天線緊貼干沙表面，圖8a所示是實測的共偏移距探地雷達剖面，橫軸為道號，縱軸為雷達波雙程走時.圖中0～1 ns之間是天線耦合直達波，頂點位于3.2 ns的雙曲線是來自PVC盒子頂端的散射波.同時可以發(fā)現(xiàn)剖面中出現(xiàn)了一些微弱的散射雜波，雜波是由于干沙局部不均勻性造成的.

其次，采集了共炮點雷達數(shù)據(jù)集，第一個發(fā)射天線位于起始位置，第一個接收位置距離發(fā)射天線的間距為0.1 m，每一個發(fā)射位置對應(yīng)21個接收位置，即單炮雷達數(shù)據(jù)包含了21道信號，接收器間的間隔是0.01 m，發(fā)射天線的移動間距為0.05 m.圖8b是實測的共炮點探地雷達數(shù)據(jù)集，它包含36個共炮點探地雷達剖面數(shù)據(jù).這些數(shù)據(jù)將用于時間域全波形反演，以此來反演地下介質(zhì)分布，從而為疊前逆時偏移成像提供初始速度模型.

4.2 雷達實測數(shù)據(jù)的全波形反演

實測數(shù)據(jù)必須經(jīng)過必要的數(shù)據(jù)預處理，包括對共炮點探地雷達數(shù)據(jù)進行直流偏移和帶通濾波，同時考慮到本文所用的時間域全波形反演算法是二維的，而實測雷達數(shù)據(jù)是三維的，故采用(Bleistein, 1986)提出的三維轉(zhuǎn)二維技術(shù)，將三維實測數(shù)據(jù)轉(zhuǎn)為二維情況下的數(shù)據(jù)，從而消除它們在幅度和相位上的差異.由于在沙中僅埋有一個PVC盒子，故而在時間域全波形反演中，采用相對介電常數(shù)為4均勻模型作為的初始相對介電常數(shù)模型(圖9a).以實測雷達數(shù)據(jù)和計算雷達數(shù)據(jù)的均方根作為反演迭代次數(shù)的目標函數(shù)，如圖10所示，經(jīng)過15次迭代后，目標函數(shù)的誤差改變小于1%.圖9b是迭代15次后的相對介電常數(shù)反演結(jié)果，圖中可清晰辨別出傾斜掩埋的PVC盒子，其形狀和位置與實際模型基本一致，傾角也非常接近45°，此外反演結(jié)果中小散射點清晰可見.

圖7 物理模擬實驗(a) PVC盒子傾斜掩埋在沙中; (b) 實驗用到的儀器：Mala ProEx控制單元及兩套1600 MHz的高頻天線;(c) PVC盒子埋于干沙中的示意圖.Fig.7 Laboratory experiment setup(a) The PVC box tilted buried in sand; (b) Instruments used to carry out radar measurements consisting of Mala ProEx control unit and two sets of 1.6 GHz antenna; (c) A sketch of the PVC box embedded in sand.

圖8 物理實驗實測數(shù)據(jù)(a) 實測的共偏移距剖面; (b) 實測的包含36條共炮點雷達剖面的雷達數(shù)據(jù)集，所有圖像均進行了歸一化.Fig.8 Physical experiments radar data(a) The measured common-offset profile; (b) The measured radar gather consisting of 36 common-shot radar profiles, all figures are normalized by their maximum amplitude.

圖9 FWI的初始介電常數(shù)模型(a) 相對介電常數(shù)定為4; (b) 15次反演迭代后的介電常數(shù)結(jié)果.Fig.9 The initial permittivity model for the FWI(a) The dielectric constant is set to 4; (b) The inversion permittivity results after 15 iterations.

圖10 以實測雷達數(shù)據(jù)和計算雷達數(shù)據(jù)間均方根值作為迭代次數(shù)的函數(shù)(均方根值經(jīng)過以做歸一化，而且15次迭代后均方根改變值小于1%)Fig.10 RMS values as a function of the iteration number for the observed radar data and the calculated data (The RMS is normalized to its maximum and after 15 iterations; the RMS misfit changes less than 1%)

圖11 (a)實測共偏移距剖面Kirchhoff偏移結(jié)果; (b)實測共炮點雷達數(shù)據(jù)得到的逆時偏移成像結(jié)果，兩圖均已歸一化處理Fig.11 (a) Kirchhoff migration result from measured common-offset profile, (b) RTM result from measured common-shot radar gather; all figures are normalized by their maximum amplitude

4.3 探地雷達實測數(shù)據(jù)的逆時偏移結(jié)果

通過時間域全波形反演獲得了地下介電常數(shù)分布情況，基于反演結(jié)果給定的初始介電常數(shù)模型，分別進行了Kirchhoff偏移成像和逆時偏移成像，其結(jié)果見圖11.圖11a是對實測的共偏移距雷達數(shù)據(jù)進行Kirchhoff偏移成像得到的結(jié)果，顯然來自PVC盒頂?shù)碾p曲線同向軸能量收斂到一個點，該聚焦點幾乎與PVC盒頂端位置吻合，然而PVC盒身在偏移圖像中并沒有得到成像，此外，成像結(jié)果中出現(xiàn)了許多散射雜點，這些散射雜點不利于雷達數(shù)據(jù)的解釋.圖11b是對實測共炮點探地雷達數(shù)據(jù)進行逆時偏移成像處理獲得的結(jié)果，從圖可知，成像結(jié)果中傾斜PVC盒子整體(白色橢圓框所示)形態(tài)及傾角與實際模型基本一致，其偏移結(jié)果明顯好于Kirchhoff偏移成像結(jié)果.

5 結(jié)論與討論

本文開發(fā)了電磁波二維逆時偏移成像算法，其中電磁波場的正向傳播和反向傳播均采用時間域有限差分方法，互相關(guān)成像條件用于獲得成像結(jié)果.考慮到全空間波場存儲內(nèi)存開銷較大，文中將成像空間外邊界做為存儲空間，記錄其各個時刻的場值，記錄的場值完全能夠恢復電磁波場的正向傳播過程.

數(shù)值模擬實驗中，構(gòu)建了橫向速度變化較大的模型，合成了共偏移距和共炮點雷達數(shù)據(jù).對于合成的共偏移距雷達數(shù)據(jù)，同時采用Kirchhoff偏移成像和逆時偏移成像算法，對比成像結(jié)果發(fā)現(xiàn)，逆時偏移成像剖面中，來自起伏界面及地下異常體的空間位置和幾何形態(tài)與實際模型一致，而Kirchhoff偏移剖面中來自起伏界面和異常體的反射波(散射波)能量聚焦不完全，剖面中仍然存在明顯的雜波，可見逆時偏移成像算法明顯好于Kirchhoff偏移算法.然而需要注意的是，逆時偏移成像算法是基于正演計算，每完成一次偏移需要進行兩次以上正演計算，而正演計算非常耗時，故而其計算效率遠低于Kirchhoff偏移算法.采用逆時偏移成像處理共炮點雷達合成數(shù)據(jù)，得到的成像剖面與該算法處理共偏移距雷達數(shù)據(jù)獲得的成像剖面基本一致，不同之處在于前者成像剖面中的子波寬度略寬于后者.

影響逆時偏移成像效果的關(guān)鍵因素之一是包括介電常數(shù)和電導率在內(nèi)的初始模型的選擇，數(shù)值實驗表明更精準的初始模型能獲得更高精度的偏移成像結(jié)果.為了獲得更精確的初始模型，本文通過時間域全波形反演算法對探地雷達數(shù)據(jù)進行全波形反演，反演后的地下相對介電常數(shù)模型作為逆時偏移成像的初始速度模型.數(shù)值模擬實驗及實際物理模擬實驗表明，通過全波形反演獲取的信息作為逆時偏移成像的初始模型是行之有效的，同時逆時偏移成像的效果在一定程度上檢驗了全波反演結(jié)果的準確性.

Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.

Bleistein N. 1986. Two-and-one-half dimensional in-plane wave propagation.GeophysicalProspecting, 34(5): 686-703.

Chang W F, McMechan G A. 1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition.Geophysics, 51(1): 67-84.

Chang W F, McMechan G A. 1994. 3-D elastic prestack, reverse-time depth migration.Geophysics, 59(4): 597-609.

Chattopadhyay S, McMechan G A. 2008. Imaging conditions for prestack reverse-time migration.Geophysics, 73(3): S81-S89.

Claerbout J F. 1971. Toward a unified theory of reflector mapping.Geophysics, 36(3): 467-481.

Docherty P C. 1991. A brief comparison of some Kirchhoff integral formulas for migration and inversion.Geophysics, 56(8): 1164-1169.

Fang N S, Ying L A. 2009. Stability analysis of FDTD to UPML for time dependent Maxwell equations.ScienceinChinaSeriesA, 52(4): 794-816.

Feng D S, Dai Q W. 2011. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.NDT&EInternational, 44(6): 495-504.Feng X, Sato M. 2004. Pre-stack migration applied to GPR for landmine detection.InverseProblems, 20(6): S99-S115.

Fisher E, McMechan G A, Annan A P, et al. 1992. Examples of reverse-time migration of signal-channel ground-penetrating radar profiles.Geophysics, 57(4): 577-586.

Gray S H. 1986. Efficient traveltime calculations for Kirchhoff migration.Geophysics, 51(8): 1685-1688.

Gray S H, May W P. 1994. Kirchhoff migration using eikonal equation traveltimes.Geophysics, 59(5): 810-817.

Guitton A, Kaelin B, Biondi B. 2006. Least-squares attenuation of reverse-time-migration artifacts.Geophysics, 72(1): S19-S23.

Guitton A, Valenciano A, Bevc D, et al. 2007. Smoothing imaging condition for shot-profile migration.Geophysics, 72(3): S149-S154. Levin S A. 1984. Principle of reverse-time migration.Geophysics, 49(5): 581-583.

Liu S X, Sato M. 2005. Transient radiation from an unloaded, finite dipole antenna in a borehole: experimental and numerical results.Geophysics, 70(6): K43-K51.

Liu S X, Zeng Z F, Xu B. 2006. FDTD simulation of ground penetrating radar signal in 3-Dimensional dispersive medium.JournalofJilinUniversity(EarthScienceEdition) (in Chinese), 36(1): 123-127.

Liu S X, Zeng Z F, Deng L. 2007. FDTD simulations for ground

penetrating radar in urban applications.JournalofGeophysicsandEngineering, 4(3): 262-267.

Liu S X, Zeng Z F.2007.Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.ChineseJ.Geophys. (in Chinese), 50(1): 320-326.

Liu S X, Wu J J, Zhou J F, et al. 2011. Numerical simulations of borehole radar detection for metal ore.IEEEGeosci.RemoteSens.Lett., 8(2): 308-312.

Liu X X, Yin X Y, Liu B H. 2013. Prestack reverse-time migration of elastic waves in porous media.ProgressinGeophys. (in Chinese), 28(6): 3165-3173.

Loewenthal D, Mufti I R. 1983. Reversed time migration in spatial frequency domain.Geophysics, 48(5): 627-635.

Meles G, Greenhalgh S, van der Kruk J, et al. 2011. Taming the non-linearity problem in GPR full-waveform inversion for high contrast media.JournalofAppliedGeophysics, 73(2): 174-186. 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.

Mufti I R, Pita J A, Huntley R W. 1996. Finite-difference depth migration of exploration-scale 3-D seismic data.Geophysics, 61(3): 776-794.

Sava P, Fomel S. 2006. Time-shift imaging condition in seismic migration.Geophysics, 71(6): S209-S217.

Sun R, McMechanz G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data.Geophysics, 66(5): 1519-1527. Taflove A, Hagness S. 2005. Computational Electrodynamics - The Finite-Difference Time-Domain Method, 3rd Ed. London: Artech House.Vidale J E. 1988. Finite-difference calculation of travel times.Bull.Seismol.Soc.Am., 78(6): 2062-2076.

Wang B L, Gao J H, Chen W C, et al. 2012. Efficient boundary storage strategies for seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 55(7): 2412-2421, doi: 10.6038/j.issn.0001-5733.2012.07.025.

Yoon K, Marfurt K J. 2006. Reverse-time migration using the Poynting vector.ExplorationGeophysics, 37(1): 102-107.

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.

Zhou H, Sato M, Liu H J. 2005. Migration velocity analysis and prestack migration of common-transmitter GPR data.IEEETransactionsonGeoscienceandRemoteSensing, 43(1): 86-91.

附中文參考文獻

劉四新, 曾昭發(fā), 徐波. 2006. 三維頻散介質(zhì)中地質(zhì)雷達信號的FDTD數(shù)值模擬. 吉林大學學報(地球科學版), 36(1): 123-127. 劉四新, 曾昭發(fā). 2007. 頻散介質(zhì)中地質(zhì)雷達波傳播的數(shù)值模擬. 地球物理學報, 50(1): 320-326.

劉欣欣, 印興耀, 劉百紅. 2013. 孔隙介質(zhì)彈性波疊前逆時偏移方法. 地球物理學進展, 28(6): 3165-3173.

王保利, 高靜懷, 陳文超等. 2012. 地震疊前逆時偏移的有效邊界存儲策略. 地球物理學報, 55(7): 2412-2421, doi: 10.6038/j.issn.0001-5733.2012.07.025.

(本文編輯 汪海英)

Reverse time migration applied to GPR data based on full wave inversion

LEI Lin-Lin1, LIU Si-Xin1, FU Lei1*, WU Jun-Jun2

1CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China2BGPInc.,ChinaNationalPetroleumCorporation,ZhuozhouHebei072751,China

Reverse-time migration (RTM) is used for subsurface imaging to handle complex velocity models including steeply dipping interfaces and dramatic lateral variations, and promises better imaging results compared to traditional migration methods such as Kirchhoff migration algorithm. RTM has been increasingly applied to seismic surveys for hydrocarbon resource explorations. Based on the similarity of kinematics and dynamics between electromagnetic waves and elastic waves, we develop a pre-stack RTM method and apply it to processing ground penetrating radar (GPR) data.The finite-difference time domain (FDTD) numerical method is used to simulate the electromagnetic wave propagation including forward and backward extrapolation. The cross-correlation imaging condition is used to obtain the final image. In order to provide a velocity model with relatively higher accuracy as the initial velocity model for RTM, we apply a full waveform inversion (FWI) in the time domain to estimate the subsurface velocity structure based on reflection radar data. For testing the effectiveness of the algorithm, we have constructed a complex geological model; and synthesized common-offset radar data and common-shot profile (CSP) radar reflection data. All data are migrated with the traditional Kirchhoff migration method and pre-stack RTM method separately. The migration results from pre-stack RTM show better coincidence with the true model. Furthermore, we have performed a physical experiment in a sandbox where a polyvinyl chloride (PVC) box is buried in the sand. The obtained common-offset radar data and common-shot radar data are migrated by using the Kirchhoff migration method and per-stack RTM algorithm separately. The per-stack RTM result shows that the RTM algorithm can obtain better imaging results.In the numerical experiments, velocity variation exists in the geological model due to the irregular interface between soil and sand. The imaging result from pre-stack RTM is better than that of the traditional Kirchhoff migration algorithm because pre-stack RTM can handle the model with horizontal velocity variation. In the RTM imaging, the shape and position of the interface and their anomalies match well with the true model, while the imaging of Kirchhoff migration cannot do so. As for the physical experiment, the Kirchhoff migration algorithm and pre-stack RTM have been applied to process the measured radar data. The imaging results from RTM show its advantage compared with the Kirchhoff migration method.

Reverse time migration; GPR; Full waveform inversion; Physical model experiment

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

10.6038/cjg20150927.

Lei L L, Liu S X, Fu L, et al. 2015. Reverse time migration applied to GPR data based on full wave inversion.ChineseJ.Geophys. (in Chinese),58(9):3346-3355,doi:10.6038/cjg20150927.

10.6038/cjg20150927

P631

2014-03-14，2015-08-02收修定稿

國家自然科學基金項目(40874073,41074076),國家高技術(shù)研究發(fā)展計劃項目(2013AA064603)及中石油東方地球物理公司中青年科技創(chuàng)新基金項目(11-06-2013)共同資助.

雷林林，女，1987年生，博士生，主要從事地球物理數(shù)據(jù)處理與正反演研究.E-mail:leilljlu@126.com

*通訊作者 傅磊，男，1985年生，講師.E-mail:fuleijlu@163.com