機(jī)載探地雷達(dá)數(shù)值模擬及逆時(shí)偏移成像

傅磊，劉四新＊，劉瀾波，吳俊軍

1吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

2康涅狄格大學(xué)城市與環(huán)境工程系，Storrs 06269

3中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司新興物探開(kāi)發(fā)處，河北涿州 072751

1 引言

探地雷達(dá)（Ground Penetrating Radar）簡(jiǎn)稱(chēng)GPR，是一種利用電磁波進(jìn)行無(wú)損探測(cè)的淺層地球物理技術(shù)，被廣泛應(yīng)用于考古學(xué)，冰川學(xué)，礦產(chǎn)勘查，以及環(huán)境與工程等領(lǐng)域.Norhidayahti（2012）等利用探地雷達(dá)在馬來(lái)西亞的仁嶺進(jìn)行考古研究；Arcone（1991）研究了南極洲冰雪覆蓋的3.2km深度內(nèi)的地質(zhì)結(jié)構(gòu)；Moore等（1999）用50MHz和200MHz雷達(dá)天線(xiàn)對(duì)挪威斯瓦爾巴特群島南部的多熱型冰川進(jìn)行了高精度勘察；Singh（2007）利用探地雷達(dá)對(duì)大型礦砂礦床資源量進(jìn)行評(píng)估.對(duì)于植被嚴(yán)重覆蓋的區(qū)域，或人類(lèi)無(wú)法到達(dá)的危險(xiǎn)地區(qū)如戰(zhàn)場(chǎng)、雷場(chǎng)等，常規(guī)的探地雷達(dá)則顯得無(wú)能為力，然而機(jī)載探地雷達(dá)（Damm et al.，2006；劉四新等，2012）（Airborne GPR）卻是一種有效的探測(cè)手段.

目前，機(jī)載探地雷達(dá)的發(fā)展和其在工程上的應(yīng)用仍然受很多因素制約.美軍寒地研究及工程實(shí)驗(yàn)室（CRREL）通過(guò)在直升機(jī)上懸吊傳統(tǒng)探地雷達(dá)系統(tǒng)（Arcone，1991），在阿拉斯加及南極洲對(duì)冰層厚度及冰層中裂隙進(jìn)行了探測(cè)研究.德國(guó)地球科學(xué)和自然資源研究所（BGR）開(kāi)發(fā)了一套特殊的機(jī)載探地雷達(dá)系統(tǒng)（Krellmann et al.，2008），它被安裝在AS350直升機(jī)上，角反射器天線(xiàn)作為發(fā)射和接收天線(xiàn)，操作頻率為150MHz，該雷達(dá)系統(tǒng)成功應(yīng)用于地質(zhì)調(diào)查.美國(guó)德克薩斯大學(xué)奧斯汀分校地球物理研究所（UTIG）開(kāi)發(fā)了固定翼機(jī)載探地雷達(dá)系統(tǒng)（Holt et al.，2006），該平臺(tái)包括了重力儀、磁力儀、激光高度計(jì)、導(dǎo)航系統(tǒng)、數(shù)碼相機(jī)以及冰雷達(dá)，該系統(tǒng)為窄帶系統(tǒng)，特點(diǎn)是速度快、功率大.

由于機(jī)載探地雷達(dá)的天線(xiàn)距離地表較高，其電磁波傳播要比常規(guī)探地雷達(dá)復(fù)雜.常規(guī)探地雷達(dá)天線(xiàn)通過(guò)與地面耦合向地下發(fā)射電磁波，而機(jī)載探地雷達(dá)天線(xiàn)首先與空氣耦合，電磁波能量在空氣中由于幾何擴(kuò)散而衰減，到達(dá)地表時(shí)，一部分能量反射（散射）回空氣中，剩下的能量繼續(xù)向地下傳播，遇到地下介質(zhì)不連續(xù)面便產(chǎn)生反射波，該反射波經(jīng)過(guò)地層衰減以及幾何擴(kuò)散，最終到達(dá)接收天線(xiàn).可見(jiàn)，相比于常規(guī)探地雷達(dá)，機(jī)載探地雷達(dá)的探測(cè)機(jī)理變得更為復(fù)雜.

影響機(jī)載探地雷達(dá)探測(cè)效果的因素主要有天線(xiàn)極化方向、天線(xiàn)高度以及地表粗糙度等，為了研究這些因素與機(jī)載探地雷達(dá)探測(cè)效果之間的關(guān)系，采用三維時(shí)間域有限差分（3D-FDTD）方法對(duì)沙漠地區(qū)地下空洞掩體機(jī)載探地雷達(dá)探測(cè)進(jìn)行數(shù)值模擬，并分析各因素對(duì)探測(cè)效果的影響.機(jī)載探地雷達(dá)探測(cè)過(guò)程中，起伏的地表會(huì)導(dǎo)致地表散射雜波的存在，這些散射雜波的能量一般較強(qiáng)，在地下目標(biāo)體埋深較淺的情況下，目標(biāo)體信號(hào)會(huì)被地表雜波掩蓋.為了消除地表雜波的影響，張蓓等（2005）提出子空間地表雜波抑制方法，并取得了一定的效果，本文使用逆時(shí)偏移成像（Chang and McMechan，1994；Wu et al.，1996；Sun et al.，2006；龍桂華等，2011；王童奎等，2012；徐興榮等，2012；孫小東等，2012；王保利等，2012；胡明順等，2013；張巖和吳國(guó)忱，2013）的方法來(lái)消除地表雜波的影響.

2 高斯隨機(jī)粗糙表面

傳統(tǒng)探地雷達(dá)探測(cè)中，與天線(xiàn)接觸的地表一般比較平坦，而機(jī)載探地雷達(dá)探測(cè)過(guò)程中遇到的地形情況比較復(fù)雜.為了模擬粗糙不平的起伏地表，引入高斯隨機(jī)粗糙表面（Kobayashi et al.，2002；法文哲和金亞秋，2010），它的功率譜表示如下：

其中σ為粗糙表面的均方根高度，它決定了粗糙表面在垂直方向的尺度；kx和ky分別為粗糙面在x和y方向的波數(shù)；lx和ly分別為沿x和y方向的相關(guān)長(zhǎng)度，它與粗糙表面在水平方向的尺度有關(guān).對(duì)于一定的水平尺度lx和ly，均方根高度越大意味著表面越粗糙.我們用自由空間中的波數(shù)k與均方根高度σ的乘積大小來(lái)度量粗糙面的粗糙度.圖1所示分別是三種不同粗糙度的高斯隨機(jī)粗糙表面，本文模擬中，選取水平相關(guān)長(zhǎng)度lx和ly均為1.25m.圖1a是kσ1為0.5的高斯隨機(jī)粗糙表面，其最大高程為1.14m，最小高程為-0.91m；圖1b是kσ2為1的高斯隨機(jī)粗糙表面，其最大高程為2.28m，最小高程為-1.82m；圖1c是kσ3為2的高斯隨機(jī)粗糙表面，其最大高程為4.56m，最小高程為-3.64m.從圖可知，當(dāng)度量值kσ越大，隨機(jī)粗糙表面的粗糙度越大.

3 數(shù)值模擬及分析

通過(guò)對(duì)典型地質(zhì)模型進(jìn)行雷達(dá)波的正演模擬，可以了解電磁波在介質(zhì)中的傳播過(guò)程，有助于對(duì)雷達(dá)數(shù)據(jù)的解釋.從麥克斯韋微分方程組出發(fā)，Yee（1966）提出了時(shí)間域有限差分求解法，它將麥克斯韋旋度方程組進(jìn)行差分離散，在時(shí)間上迭代推進(jìn)求解空間中每一個(gè)時(shí)刻的電磁場(chǎng)值.

本文利用三維時(shí)間域有限差分法（劉四新等，2005；2007；2011；2012；Taflove et al.，2005）對(duì)起伏的沙漠地區(qū)地下空洞掩體機(jī)載探地雷達(dá)探測(cè)進(jìn)行了數(shù)值模擬.圖2a是建立的三維地質(zhì)模型，其空間尺度在x方向?yàn)?16m，在y方向?yàn)?0m，在z方向?yàn)?0m.用度量值kσ2=1高斯隨機(jī)粗糙表面代表起伏的沙漠表面，沙漠表面所在的平均高程為15m，即干沙的厚度約為15m，其相對(duì)節(jié)電常數(shù)為4，電導(dǎo)率為0.1mS·m-1.在干沙地層內(nèi)部有兩個(gè)沿y走向空洞掩體，分別用充滿(mǎn)空氣的兩個(gè)圓柱體表示.圖2b是y=15m處的二維切片模型，從圖可知，兩個(gè)空洞的xz截面中心坐標(biāo)分別為（48，10.5）和（68，7.5）m，空洞的截面半徑均為1.5m.圖中虛線(xiàn)代表天線(xiàn)的飛行軌跡，即合成數(shù)據(jù)的記錄位置，其距離平均地表的平均高度約12m，水平方向30～90m，每隔0.5m記錄一道數(shù)據(jù).

圖2 （a）沙漠地區(qū)空洞三維模型；（b）相應(yīng)的二維模型，色標(biāo)為相對(duì)介電常數(shù)，無(wú)量綱Fig.2 （a）Three-dimensional model of cavity embedded in the sand；（b）Corresponding two-dimensional model

3.1 波場(chǎng)快照分析

整個(gè)模型空間離散為0.2m的立方體單元格，選取中心頻率為50MHz的雷克子波作為激勵(lì)源，時(shí)間采樣間隔為0.2ns，時(shí)窗為250ns.當(dāng)激勵(lì)源位于合成記錄位置的中點(diǎn)，采用x方向的電場(chǎng)極化，截取不同時(shí)刻的波場(chǎng)快照，見(jiàn)圖3，白色星號(hào)表示激勵(lì)源所在位置，白色曲線(xiàn)代表起伏的沙漠地表，兩個(gè)圓形分別代表地下目標(biāo)體.觀察y=15m截面不同時(shí)刻x方向的電場(chǎng)，圖3a是60ns時(shí)刻的波場(chǎng)快照，可以看到此時(shí)波前剛剛傳到起伏的沙漠表面，在高程大于平均高程（15m）的60m處，已經(jīng)有部分能量散射到空氣中；圖3b是100ns時(shí)刻的波場(chǎng)快照，從圖可知，由于起伏地表的影響，一部分能量散射回空氣中，并且散射波雜亂交錯(cuò)，剩下的能量穿透沙漠地表繼續(xù)往前傳播，其波前面已經(jīng)觸及到淺部異常體，并且波前面變的不規(guī)整；圖3c是140ns時(shí)刻的波場(chǎng)快照，此時(shí)波前面已經(jīng)傳到深部異常體，發(fā)現(xiàn)兩個(gè)異常體的反射波能量弱于空氣中的散射雜波能量；圖3d是180ns時(shí)刻的波場(chǎng)快照，其波前面可以看做是在起伏地表上的許多二次源的波前面疊加效果，可見(jiàn)起伏地表會(huì)造成地表散射雜波，而且其能量往往強(qiáng)于地下目標(biāo)體反射回波信號(hào).

3.2 不同天線(xiàn)極化方向?qū)μ綔y(cè)效果影響分析

機(jī)載雷達(dá)探測(cè)實(shí)際過(guò)程中，由于直升機(jī)的飛行速度、當(dāng)?shù)仫L(fēng)向及風(fēng)速的變化等影響，天線(xiàn)的姿態(tài)會(huì)發(fā)生變化，如左右旋轉(zhuǎn)，上下?lián)u擺.根據(jù)機(jī)載雷達(dá)的特殊性，本文研究了四種典型的天線(xiàn)配置方式，如圖4所示，天線(xiàn)的飛行軌跡為沿x方向的帶右端箭頭的虛線(xiàn)，其中FBy代表天線(xiàn)平行并分別配置于直升機(jī)的頭部和尾部，天線(xiàn)的極化方向?yàn)閥方向；FBx代表天線(xiàn)同軸共線(xiàn)配置于直升機(jī)的頭部和尾部，天線(xiàn)的極化方向?yàn)閤方向；LRy代表天線(xiàn)共線(xiàn)配置于直升機(jī)的左右兩翼，天線(xiàn)的極化方向?yàn)閥方向；LRx代表天線(xiàn)平行配置于直升機(jī)的左右兩翼，天線(xiàn)的極化方向?yàn)閤方向.

圖3 不同時(shí)刻的波場(chǎng)快照Fig.3 Wave field snapshots at different time

圖4 機(jī)載探地雷達(dá)四種不同的天線(xiàn)配置方式Fig.4 Four different antenna configurations for airborne GPR survey

根據(jù)四種典型的天線(xiàn)配置方式，研究了不同極化方式對(duì)機(jī)載探地雷達(dá)探測(cè)地下異常體的影響.選取粗糙度為kσ2=1的高斯隨機(jī)表面代表起伏的沙漠表面，固定天線(xiàn)距離沙漠地表高度為12m，模型參數(shù)與圖2所示三維模型保持一致.圖5是四種典型的天線(xiàn)配置情況下獲得的機(jī)載探地雷達(dá)探測(cè)地下空洞掩體的雷達(dá)剖面，圖中的直達(dá)波均已去除.圖5a是天線(xiàn)分別配置于直升機(jī)頭部和尾部，y方向極化時(shí)獲得的模擬雷達(dá)剖面，從圖可知，地表的散射雜波主要集中在110～130ns，其能量遠(yuǎn)遠(yuǎn)強(qiáng)于地下目標(biāo)體能量，淺部異常體反射信號(hào)完全淹沒(méi)于雜波中，但是來(lái)自深部異常體的雙曲形態(tài)反射信號(hào)仍然能被識(shí)別；圖5b是天線(xiàn)同軸共線(xiàn)分別位于直升機(jī)頭部和尾部，x方向極化時(shí)獲得的模擬雷達(dá)剖面，從中可知，地表的散射雜波能量仍然比較強(qiáng)，但是淺部的異常體目標(biāo)反射信號(hào)以及深部的目標(biāo)反射信號(hào)均可以被識(shí)別；圖5c是天線(xiàn)分別位于直升機(jī)左右兩側(cè)，x方向極化時(shí)獲得的模擬雷達(dá)剖面，強(qiáng)地表雜波淹沒(méi)了淺部目標(biāo)體信號(hào)；圖5d是天線(xiàn)分別位于直升機(jī)左右兩側(cè)，y方向極化時(shí)獲得的模擬雷達(dá)剖面，圖中來(lái)自地表的散射雜波能量主要集中在110～130ns范圍內(nèi)，淺部目標(biāo)體仍然可以識(shí)別.

綜上可知，當(dāng)天線(xiàn)的極化方向與異常體的走向平行時(shí)，地下淺部目標(biāo)體的信號(hào)被來(lái)自地表的強(qiáng)雜波覆蓋，地下淺部異常體目標(biāo)很難被識(shí)別，而當(dāng)天線(xiàn)的極化方向與異常體垂直時(shí)，地下異常體目標(biāo)反射波能量增強(qiáng)，雖然地表雜波能量較強(qiáng)，但是淺部異常體目標(biāo)仍然可以被識(shí)別.因此，在機(jī)載探地雷達(dá)實(shí)際探測(cè)過(guò)程中，在保證測(cè)線(xiàn)與異常體走向垂直的同時(shí)，天線(xiàn)飛行過(guò)程中應(yīng)盡量保證天線(xiàn)的極化方向與異常體走向垂直.

3.3 天線(xiàn)高度對(duì)探測(cè)效果影響分析

圖5 四種不同天線(xiàn)配置方式下模擬的雷達(dá)剖面Fig.5 Simulated profile results for four antenna configurations of Airborne GPR

機(jī)載探地雷達(dá)探測(cè)過(guò)程中，為了保證天線(xiàn)的飛行姿態(tài)穩(wěn)定，天線(xiàn)的飛行軌跡應(yīng)避開(kāi)地表植被以及地面建筑物，即天線(xiàn)距離地面高度不能過(guò)低.為了研究天線(xiàn)距離地表高度對(duì)探測(cè)結(jié)果的影響，選取粗糙度為kσ2=1的高斯隨機(jī)粗糙表面代表起伏的沙漠表面，天線(xiàn)位于直升機(jī)左右兩側(cè)，激勵(lì)源的極化方向?yàn)閤方向極化，模型中的其他參數(shù)保持不變，改變天線(xiàn)的高度，分別從6m到8m、10m、12m（如圖6所示）模擬了不同高度情況下沙漠地區(qū)地下空洞掩體機(jī)載雷達(dá)探測(cè)雷達(dá)剖面，圖中直達(dá)波均已去除.圖6a是當(dāng)天線(xiàn)距離地表平均高度為6m時(shí)得到的合成雷達(dá)剖面，從圖可知，位于淺部的異常體及深部的異常體反射波雙曲形態(tài)同向軸清晰可見(jiàn)，來(lái)自地表的散射雜波能量較強(qiáng)，圖6b是天線(xiàn)距離地表平均高度為8m時(shí)的合成雷達(dá)剖面，圖6c是天線(xiàn)距離地表平均高度為10m時(shí)的合成雷達(dá)剖面，圖6d是當(dāng)天線(xiàn)距離地表平均高度為12m時(shí)獲得的雷達(dá)剖面，從圖可知，隨著天線(xiàn)距離地表的距離增加，來(lái)自地下異常體的反射波的能量越來(lái)越弱，其雙曲線(xiàn)型同向軸的振幅變小，同向軸越來(lái)越模糊，而來(lái)自地表的反射雜波的能量仍然很強(qiáng).對(duì)于機(jī)載探地雷達(dá)探測(cè)地下異常體，發(fā)射天線(xiàn)發(fā)射的電磁波能量在空氣中由于幾何擴(kuò)散而衰減，其能量與距離的二次方成反比，隨著天線(xiàn)高度增加，到達(dá)地表的能量減少，從而降低了信噪比，不利于地下異常體的探測(cè).機(jī)載雷達(dá)野外實(shí)際探測(cè)中，為了得到更高分辨率的雷達(dá)圖像，在保證天線(xiàn)飛行姿態(tài)穩(wěn)定的情況下，天線(xiàn)高度盡可能降低.

3.4 粗糙表面對(duì)探測(cè)效果影響分析

傳統(tǒng)探地雷達(dá)探測(cè)過(guò)程中，天線(xiàn)與地面耦合，大部分能量直接向地下傳播；而機(jī)載探地雷達(dá)探測(cè)過(guò)程中，天線(xiàn)位于距離地面較高的空中，發(fā)射天線(xiàn)發(fā)出的電磁波在空氣中由于幾何擴(kuò)撒而損耗能量，當(dāng)電磁波到達(dá)地表時(shí)，粗糙不平的地表會(huì)導(dǎo)致電磁波的空間各個(gè)不同方向的散射，為了研究不同粗糙度表面對(duì)探測(cè)結(jié)果的影響，把天線(xiàn)距離地面的高度固定為12m，天線(xiàn)位于直升機(jī)左右兩側(cè)，激勵(lì)源為x方向極化，其他模型參數(shù)保持不變，改變沙漠地表的粗糙度.圖7是沙漠地表不同粗糙度情況下獲得的機(jī)載探地雷達(dá)探測(cè)地下空洞的模擬結(jié)果.圖中水平方向是沿x方向的距離，垂直方向是電磁波雙程傳播時(shí)間，圖中所有剖面的直達(dá)波都已去除.

圖6 不同天線(xiàn)高度下模擬的雷達(dá)剖面Fig.6 Simulated radar profile of different antenna height

圖7a是沙漠地表水平情況下得到的雷達(dá)剖面，從圖可知，水平同向軸是來(lái)自沙漠地表的反射波，而兩條雙曲線(xiàn)分別代表來(lái)自地下空洞的反射波，兩條雙曲線(xiàn)的頂點(diǎn)位置反映了異常體所在的水平位置，從剖面可知，第一個(gè)異常體水平位置為48m，第二個(gè)異常體的水平位置為68m，這與模型是一致的.圖7b是沙漠表面粗糙度為kσ1=0.5時(shí)的雷達(dá)剖面，相比于圖7a，可以發(fā)現(xiàn)，雖然來(lái)自沙漠表面以及地下空洞的反射波均能清晰地識(shí)別，但是雷達(dá)剖面在110～130ns之間出現(xiàn)了很多地表散射雜波，這些雜波的存在勢(shì)必給剖面的解釋帶來(lái)困惑，另外來(lái)自地下空洞的反射波的振幅稍有變小.圖7c是沙漠表面粗糙度為kσ2=1時(shí)的雷達(dá)剖面，仍然可以清晰地識(shí)別來(lái)自沙漠表面的反射波，然而由于沙漠表面粗糙度的增加，導(dǎo)致電磁波在地面的散射更嚴(yán)重，穿透到地下的電磁波能量更少，所以來(lái)自地下空洞的反射同向軸變得模糊，位于110～130ns的散射波能量變的更強(qiáng).圖7d是沙漠表面粗糙度為kσ3=2的雷達(dá)剖面，來(lái)自地表的散射雜波能量進(jìn)一步增強(qiáng)，而來(lái)自地下異常體的反射波的雙曲線(xiàn)形態(tài)已經(jīng)很難被識(shí)別.從該剖面中，我們很難定位地下異常體的位置.

由此可見(jiàn)，粗糙的地表會(huì)導(dǎo)致電磁波的散射，地表散射而損耗的電磁波能量越多，散射雜波的能量越強(qiáng)，同時(shí)，穿透到地下的電磁波能量相應(yīng)的就越少，來(lái)自地下目標(biāo)體的反射波能量就越弱.反映在雷達(dá)剖面中圖像分辨率降低，不利于機(jī)載探地雷達(dá)探測(cè)的開(kāi)展.

4 逆時(shí)偏移成像

雷達(dá)電磁波與地震波在運(yùn)動(dòng)學(xué)上存在相似性，反射地震上的偏移技術(shù)可以借鑒到探地雷達(dá)數(shù)據(jù)處理中.偏移的目的是為了使雷達(dá)剖面中傾斜界面共深度映像聚焦，使繞射波歸位.為了消除機(jī)載探地雷達(dá)探測(cè)中由于起伏地形產(chǎn)生的地表散射雜波影響，借鑒反射地震中的逆時(shí)偏移成像（RTM）技術(shù).

圖7 機(jī)載探地雷達(dá)于沙漠地區(qū)探測(cè)地下空洞的模擬結(jié)果Fig.7 Radar profile of airborne GPR detecting cavities at desert area

地震上逆時(shí)偏移可分為疊前逆時(shí)偏移和疊后逆時(shí)偏移.疊前偏移是把共炮點(diǎn)道集記錄或共偏移距道集記錄中的反射波歸位到產(chǎn)生它們的反射界面上并使繞射波收斂到產(chǎn)生它的繞射點(diǎn)上，最后得到能夠反映界面反射系數(shù)特點(diǎn)的并正確歸位了的地震波形剖面.疊后偏移是在水平疊加剖面的基礎(chǔ)上進(jìn)行的，針對(duì)水平疊加剖面上存在的傾斜反射層不能正確地歸位和繞射波不能完全收斂的問(wèn)題，采用了爆炸反射面的概念來(lái)實(shí)現(xiàn)傾斜反射層的正確歸位和繞射波的完全收斂.相對(duì)其他方法而言，逆時(shí)偏移用全程波方程對(duì)波場(chǎng)延拓，避免對(duì)波動(dòng)方程的近似，因此沒(méi)有傾角限制，原理上可以利用轉(zhuǎn)換波、棱鏡波、或者多次反射波成像，并獲得更精確的振幅等動(dòng)力學(xué)信息，實(shí)現(xiàn)保幅成像，還可以更好地對(duì)復(fù)雜速度成像進(jìn)行更細(xì)化更精確的估計(jì)，成像方法不受介質(zhì)速度變化的影響，能夠?qū)?fù)雜區(qū)域進(jìn)行較準(zhǔn)確的成像.

本文借鑒反射地震疊前偏移技術(shù)，采用和勘探地震中類(lèi)似的共炮集觀測(cè)方式合成雷達(dá)數(shù)據(jù).為了合成共炮集數(shù)據(jù)，采用1個(gè)發(fā)射天線(xiàn)及8個(gè)接收天線(xiàn)，接收天線(xiàn)間距為1m，發(fā)射天線(xiàn)與第一個(gè)接收天線(xiàn)的偏移距為3m，發(fā)射天線(xiàn)每次移動(dòng)距離為1m，這樣可以獲得一個(gè)共炮集數(shù)據(jù)集.合成數(shù)據(jù)時(shí)，天線(xiàn)距離地面的高度固定為12m，激勵(lì)源為x方向極化，其他模型參數(shù)保持不變，改變沙漠地表的粗糙度，獲得不同粗糙度情況下共炮集雷達(dá)數(shù)據(jù)集.

圖8 不同地表粗糙度情況下雷達(dá)剖面偏移成像結(jié)果Fig.8 GPR migration results for different roughness surface

采用二維RTM算法，對(duì)共炮集雷達(dá)數(shù)據(jù)集進(jìn)行逆時(shí)偏移成像，同時(shí)為了對(duì)比說(shuō)明RTM的成像效果，本文給出了基爾霍夫偏移成像（Duquet et al.，2000；Eaton et al.，2002；Zhu and Lines，1998）結(jié)果，如圖8所示.圖中左邊一欄四個(gè)圖是來(lái)自基爾霍夫偏移成像的結(jié)果，右邊一欄是來(lái)自二維逆時(shí)偏移成像的結(jié)果.圖8a是水平地表情況下基爾霍夫偏移成像結(jié)果，其地下目標(biāo)體的反射波聚焦較好，聚焦位置基本反映了地下目標(biāo)體的空間位置，圖8e是地表情況下逆時(shí)偏移成像結(jié)果，從中可知，地下目標(biāo)體的空間位置、幾何形態(tài)及尺寸與實(shí)際模型基本一致；圖8b是在kσ1=0.5的粗糙表面情況下獲得的基爾霍夫偏移成像結(jié)果，圖中散射雜波仍然雜亂無(wú)章，散射雜波并沒(méi)有得到聚焦，雜波的存在仍然會(huì)影響地下異常體的判斷，圖8f是對(duì)應(yīng)的逆時(shí)偏移結(jié)果，來(lái)自地表的散射雜波歸位，地下目標(biāo)體清晰可見(jiàn)，地下目標(biāo)體的空間位置、幾何形態(tài)、尺寸以及起伏地表與實(shí)際模型基本一致；同樣的，圖8g和圖8f中的偏移成像結(jié)果明顯好于圖8c和圖8d的.

基爾霍夫偏移技術(shù)只能讓地下反射目標(biāo)體信號(hào)聚焦，聚焦位置大致反映目標(biāo)體空間位置，無(wú)法聚焦的散射雜波仍然會(huì)對(duì)數(shù)據(jù)解釋造成困擾，隨著地表粗糙度增加，地下目標(biāo)體聚焦能量減弱，而逆時(shí)偏移成像技術(shù)不但可以使地下目標(biāo)體反射波聚焦，刻畫(huà)地下目標(biāo)體空間位置、幾何形態(tài)及尺寸大小，而且可以讓地表散射雜波歸位，同時(shí)起伏的地形亦得到刻畫(huà).

5 結(jié)論

（1）數(shù)值模擬試驗(yàn)表明天線(xiàn)的極化方向與異常體走向垂直時(shí)更利于探測(cè)地下目標(biāo)體，因此機(jī)載雷達(dá)野外探測(cè)時(shí)，為了得到更高分辨率的雷達(dá)圖像，測(cè)線(xiàn)與異常體走向垂直的同時(shí)，盡可能地保證天線(xiàn)的極化方向與異常體走向垂直.

（2）隨著天線(xiàn)高度的增加，從發(fā)射天線(xiàn)發(fā)射出的電磁波能量由于幾何擴(kuò)散而衰減更多，不利于探測(cè)，實(shí)際野外探測(cè)中，保證天線(xiàn)姿態(tài)穩(wěn)定的情況下，飛行高度盡可能地降低.

（3）地表粗糙度的增加會(huì)導(dǎo)致地表散射雜波增強(qiáng)，造成傳播到地下的電磁波能量減少，削弱地下目標(biāo)體信號(hào)，雷達(dá)剖面的分辨率降低，不利于地下異常體目標(biāo)識(shí)別.

（4）成功借鑒反射地震上逆時(shí)偏移成像技術(shù)，使地下目標(biāo)體反射波聚焦，聚焦位置很好地刻畫(huà)了異常體的空間位置、幾何形態(tài)及尺寸大小，同時(shí)，可以使地表散射雜波歸位，起伏的地形在偏移后的剖面中得到刻畫(huà)；隨著地表粗糙度的增加，基爾霍夫偏移成像結(jié)果更差，而逆時(shí)偏移成像結(jié)果仍然能夠很好地刻畫(huà)起伏地形及異常體.

Arcone S A.1991.Dielectric constant and layer-thickness interpretation of helicopter-borne short-pulse radar waveforms reflected from wet and dry river-ice sheets.IEEE Transactions on Geoscience and Remote Sensing，29（5）：768-777.

Arcone S A.1996.High resolution of glacial ice stratigraphy；A ground-penetrating radar study of Pegasus Runway，McMurdo Station，Antarctica.Geophysics，61（6）：1653-1663.

Chang W F，McMechan G A.1994.3-D elastic prestack，reversetime depth migration.Geophysics，59（4）：597-609.

Damm V，Casassa G，Blindow N，et al.2006.Airborne GPR measurements in the Atacama Desert first results and constraints using a 150MHz pulse radar for groundwater exploration.11th International Conference on Ground Penetrating Radar，19-22.

Docherty P.1991.A brief comparison of some Kirchhoff Integral formulas for migration and inversion.Geophysics，56（8）：1164-1169.

Duquet B，Marfurt K J，Dellinger J A.2000.Kirchhoff modeling，inversion for reflectivity，and subsurface illumination.Geophysics，65（4）：1195-1209.

Eaton D W，Hope J，Bohlen T.2002.An elastic Kirchhoff method for computing synthetic seismograms：Modeling the receiverfunction response of an irregular Moho surface.Seismol.Res.Lett.，73（2）：232.

Fa W Z，Jin Y Q.2010.Simulation of radar sounder echo from lunar surface structure.Sci.China Earth Sci.，40（4）：473-485.

Holt J W，Peters M E，Morse D L，et al.2006.Identifying and characterizing subsurface echoes in airborne radar sounding data from a high-clutter environment in Taylor Valley，Antarctica.∥Proc.the 11th Int.Conf.on Ground Penetrating Radar，AIR.3-5.

Hu M S，Pan D M，Dong S H，et al.2013.Optimization for random boundary building in reverse-time migration based on scattering theory.Progress in Geophys.（in Chinese），28（4）：2069-2077.

Kobayashi T，Oya H，Ono T.2002.A-scope analysis subsurface radar sounding of lunar mare region.Earth，Planets and Space，54（10）：973-982.

Krellmann Y，Lentz H，Triltzsch G.2008.Stepped-frequency radar system in gating mode：An Experiment as a new helicopterborne GPR system for geological applications.IEEE Geosci.and Remote Sens.Symp.1：153-156.

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.Journal of Jilin University （Earth Science Edition），36（1）：123-127.

Liu S X，F(xiàn)eng Y Q，F(xiàn)u L，et al.2012.Advances and numerical simulation of airborne ground penetrating radar.Progress in Geophys.（in Chinese），27（2），0727-0735.

Long G H，Zhao Y B，Li X F，et al.2011.Accelerating 3D staggered-grid finite-difference seismic wave modeling on GPU cluster.Progress in Geophys.（in Chinese），26（6）：1938-1949.

Liu S X，Zeng Z F，Deng L.2007.FDTD simulations for ground penetrating radar in urban applications.Journal of Geophysics and Engineering，4（3）：262-267.

Liu S X，Wu J J，Zhou J F，et al.2011.Numerical simulations of borehole radar detection for metal ore.IEEE Geosci.and Remote Sens.Lett.，8（2）：308-312.

Liu S X，Wu J J，Dong H，et al.2012.The experimental results and analysis of a borehole radar prototype.Journal of Geophysics and Engineering，9（2）：201.

Moore J C，P?lli A，Ludwig F，et al.1999.High-resolution hydrothermal structure of Hansbreen，Spitsbergen，mapped by ground-penetrating radar.Journal of Glaciology，45（151）：524-532.

Muztaza N M，Saidin M M，Azwin I N，et al.2012.Archaeological structure detection using 3DGPR survey in Jeniang，Kedah，Malaysia.∥International Proceedings of Chemical，Biological and Environmental Sciences，36：49-53.

Singh N K.2007.Ground penetrating radar（GPR）mineral base profiling and orebody optimization.∥The 6th International Heavy Minerals Conference.Hluhluwe，South Africa，185-194.

Sun R，McMechan G A，Lee C S，et al.2006.Prestack scalar reverse-time depth migration of 3Delastic seismic data.Geophysics，71（5）：S199-S207.

Sun X D，Li Z C，Wang X L.2012.Pre-stack reverse-time migration using a finite difference method based on triangular grids.Progress in Geophys.（in Chinese），27（5）：2077-2083.

Taflove A，Hagness S.2005.Computational Electrodynamics：The Finite-Difference Time-Domain Method （3rd ed）.London：Artech House.

Wang B L，Gao J H，Chen W C，et al.2012.Efficient boundary storage strategies for seismic reverse time migration.Chinese J.Geophys.（in Chinese），55（7）：2412-2421.

Wang T K，Li Y，Guo A H，et al.2012.The study of reverse-time migration technology in Nanpu 1st structure.Progress in Geophys.（in Chinese），27（6）：2541-2547.

Xu X R，Wang X W，Wang Y C，et al.2012.Study and application of imaging condition for reverse-time migration based on wavefields separation.Progress in Geophys.（in Chinese），27（5）：2084-2090.

Wu W J，Lines L R，Lu H X.1996.Analysis of higher-order，finite-difference schemes in 3-D reverse-time migration.Geophysics，61（3）：845-856.

Yee K.1966.Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEE Transactions on Antennas and Propagation，14（3）：302-307.

Zhang Y，Wu G C.2013.Review of prestack reverse-time migration in TTI media.Progress in Geophys.（in Chinese），28（1）：0409-0420.

Zhang B，Liu J X，Wu R B.2005.Research on subspace-based GPR ground bounce removal method.Signal Processing，21（4）：514-517.

Zhu J M，Lines L R.1998.Comparison of Kirchhoff and reversetime migration methods with applications to prestack depth imaging of complex structures.Geophysics，63（4）：1166-1176.

附中文參考文獻(xiàn)

法文哲，金亞秋.2010.雷達(dá)探測(cè)儀對(duì)月球次表層結(jié)構(gòu)的探測(cè)模擬方法.中國(guó)科學(xué)：地球科學(xué)，40（4）：473-485.

胡明順，潘冬明，董守華等.2013.基于散射理論的逆時(shí)偏移隨機(jī)邊界構(gòu)建策略?xún)?yōu)化分析.地球物理學(xué)進(jìn)展，28（4）：2069-2077.

劉四新，曾昭發(fā)，徐波.2006.三維頻散介質(zhì)中地質(zhì)雷達(dá)信號(hào)的FDTD數(shù)值模擬.吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），36（1）：123-127.

劉四新，馮彥謙，傅磊等.2012.機(jī)載探地雷達(dá)的進(jìn)展以及數(shù)值模擬.地球物理學(xué)進(jìn)展，27（2）：727-735.

龍桂華，趙宇波，李小凡等.2011.三維交錯(cuò)網(wǎng)格有限差分地震波模擬的GPU集群實(shí)現(xiàn).地球物理學(xué)進(jìn)展，26（6）：1938-1949.

孫小東，李振春，王小六.2012.三角網(wǎng)格有限差分法疊前逆時(shí)偏移方法研究.地球物理學(xué)進(jìn)展，27（5）：277-2083.

王保利，高靜懷，陳文超等.2012.地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略.地球物理學(xué)報(bào)，55（7）：2412-2421.

王童奎，李瑩，郭愛(ài)華.2012.逆時(shí)偏移技術(shù)在南堡1號(hào)構(gòu)造中的應(yīng)用研究.地球物理學(xué)進(jìn)展，27（6）：2541-2547.

徐興榮，王西文，王宇超等.2012.基于波場(chǎng)分離理論的逆時(shí)偏移成像條件研究及應(yīng)用.地球物理學(xué)進(jìn)展，27（5）：2084-2090.

張巖，吳國(guó)忱.2013.TTI介質(zhì)疊前逆時(shí)偏移成像研究綜述.地球物理學(xué)進(jìn)展，28（1）：409-420.

張蓓，劉家學(xué)，吳仁彪.2005.探地雷達(dá)子空間地雜波抑制方法研究.信號(hào)處理，21（Z1）：514-517.