王敏玲,廖天元,王洪華,龔俊波,張 智,熊 彬

(桂林理工大學(xué) a.地球科學(xué)學(xué)院;b.廣西隱伏金屬礦產(chǎn)勘查重點實驗室,廣西 桂林 541006)

0 引 言

偏移作為實測數(shù)據(jù)向地下真實結(jié)構(gòu)轉(zhuǎn)換的重要技術(shù),是探地雷達(ground penetrating radar, GPR)高精度成像方法之一,偏移質(zhì)量的優(yōu)劣將直接影響后續(xù)解釋的準(zhǔn)確性[1]。目前,常用的GPR偏移方法主要有Kirchhoff偏移[2-3]、頻率域波數(shù)偏移[4-5]、逆時偏移[6-7]等。眾多偏移方法中,逆時偏移具有精度高、相位準(zhǔn)確、不受橫向變速和高陡傾角限制、能精確地將任意方向傳播的電磁波場歸位到其真實地下空間位置等優(yōu)點,在GPR成像領(lǐng)域得到廣泛關(guān)注[8]。Fisher等[9]首次將逆時偏移技術(shù)應(yīng)用于GPR數(shù)據(jù)處理中;傅磊等[10]通過對比逆時偏移與Kirchhoff偏移結(jié)果，驗證了逆時偏移算法的優(yōu)越性；王敏玲等[11]為避免基于互相

關(guān)成像條件的逆時偏移計算中源點正傳電磁波場的重建，將激發(fā)振幅成像條件應(yīng)用于GPR逆時偏移成像；Bradford等[12]針對起伏地表模型的逆時偏移問題，提出了一種考慮電磁波振幅補償?shù)哪鏁r偏移算法并應(yīng)用于合成及實測GPR數(shù)據(jù)處理中；Lu等[13]為克服地表起伏對逆時偏移成像的影響，將基于貼體網(wǎng)格剖分技術(shù)應(yīng)用于源點正傳和接收點反傳電磁波場的計算，提高了起伏地表GPR模型的成像精度；朱尉強等[14]基于時間反轉(zhuǎn)不變原則，通過修改逆時傳播方程，提出了一種考慮GPR衰減補償逆時偏移方法。

上述GPR逆時偏移大都采用互相關(guān)成像條件,需要計算地下介質(zhì)中每個成像點處所有時刻的源點正傳電磁波場和接收點反傳電磁波場的互相關(guān)。由于地下介質(zhì)參數(shù)存在差異, 成像點處的源點正傳電磁波場和接收點反傳電磁波場中不但包含目標(biāo)體的反射信號, 而且也包含介質(zhì)之間相互反射和背景散射引起的雜波。 因此, 易在非反射點處產(chǎn)生較強的低頻噪聲, 降低目標(biāo)體的成像質(zhì)量[15]。 目前, 去除逆時偏移剖面中的雜波干擾和低頻噪聲方法有3種[16]:第1種是濾波[17-19],即考慮低頻噪聲與有效信號在波數(shù)域內(nèi)的差異,利用導(dǎo)數(shù)濾波,Laplace濾波、空間濾波等濾波方法對噪聲與有效信號進行分離,這些濾波方法在一定程度上可以去除假象,但對背景散射產(chǎn)生的雜波干擾的去除顯得無能為力。第2種是控制方向成像[20-21],即從假象產(chǎn)生的原因出發(fā),在互相關(guān)成像條件基礎(chǔ)上,利用不同傳播方向的正傳和反傳電磁波場進行成像,主要有坡印廷矢量法和電磁波場分解法。 當(dāng)速度模型復(fù)雜時,空間上某一點的電磁波場是很多方向電磁波場的疊加,采用坡印廷矢量不能分別描述復(fù)雜結(jié)構(gòu)的每個方向的電磁波場;電磁波場分解法需要對每個時刻、空間每個點進行傅里葉變換,計算量巨大。 第3種是根據(jù)偏移模型的速度信息,估計各成像點的主反射發(fā)生時間界限,并根據(jù)時間界限設(shè)計出該成像點的時間窗, 對時間窗內(nèi)的正傳和反傳電磁波場作歸一化互相關(guān)計算, 并將計算結(jié)果作為該成像點的零延遲互相關(guān)逆時偏移成像的權(quán)值系數(shù)以壓制時間窗外的雜波干擾[15]。 Yang 等[22]將此方法應(yīng)用于鉆孔雷達的逆時偏移成像中, 數(shù)值驗證該方法具有計算速度快、 壓制低頻噪聲和雜波干擾效果好、 分辨率高等優(yōu)點。

為此,本文在上述理論基礎(chǔ)上,將時間窗互相關(guān)成像條件應(yīng)用于地面GPR隨機介質(zhì)逆時偏移成像中。正傳和反傳電磁波場采用時域有限差分法(finite difference time domain method, FDTD)進行計算, 單軸各向異性完全匹配層(uniaxial perfectly matched layer, UPML)邊界條件用于吸收模型截斷邊界處的超強反射波。時間窗內(nèi)的歸一化互相關(guān)成像結(jié)果作為所有時刻的歸一化互相關(guān)成像結(jié)果的權(quán)值,獲取逆時偏移剖面。應(yīng)用基于時間窗互相關(guān)成像條件的GPR逆時偏移算法對典型GPR模型的多偏移距數(shù)據(jù)和實測數(shù)據(jù)進行逆時偏移成像,并與歸一化互相關(guān)成像條件的逆時偏移結(jié)果進行對比,分析了基于時間窗互相關(guān)成像條件的GPR逆時偏移算法在壓制低頻噪聲、去除雜波干擾和提高成像精度方面的優(yōu)越性。

1 GPR逆時偏移成像原理

GPR逆時偏移的基本原理是將地表記錄到的接收點電磁波場在時間軸上進行逆向傳播,當(dāng)電磁波場逆推至零時刻,則所有反射波與繞射波的能量都回到最初被反射和繞射的空間位置,然后應(yīng)用成像條件可獲取最終的偏移成像剖面[11]。因此,逆時偏移的實現(xiàn)過程主要分為3步:1)源點正傳電磁波場的計算;2)接收點反傳電磁波場的計算;3)成像條件的應(yīng)用。其中,前兩步可利用相同的數(shù)值模擬方法進行計算,本文采用基于UPML邊界條件的FDTD來實現(xiàn)源點正傳電磁波場和接收點反傳電磁波場的計算。

成像條件的應(yīng)用是開展高效高精度GPR逆時偏移的關(guān)鍵,它的好壞將直接影響成像質(zhì)量。目前GPR逆時偏移中應(yīng)用最廣泛的成像條件是Clearbout等[23]提出的零延時互相關(guān)成像條件,其基本原理是根據(jù)成像點的入射波與反射波到時相同的原理,通過計算源點正傳電磁波場和接收點反傳電磁波場的零延時互相關(guān),從而獲得成像剖面。標(biāo)準(zhǔn)的零延遲互相關(guān)成像條件公式為

(1)

零延遲互相關(guān)成像條件具有理論簡單、實現(xiàn)簡易、計算量小等優(yōu)點,但其最主要不足在于直達波和背景散射波的互相關(guān)會產(chǎn)生很強的低頻噪聲,尤其在偏移速度模型中存在相對介電常數(shù)比較明顯的分界面時,這種噪聲會嚴(yán)重干擾淺地表結(jié)構(gòu)的描述,使得成像結(jié)果模糊不可辨認(rèn),且深部結(jié)構(gòu)的成像能量非常微弱。

鑒于互相關(guān)成像條件的不足,Schleicher等[24]在互相關(guān)成像條件的基礎(chǔ)上,提出了歸一化互相關(guān)成像條件,其表達式為

(2)

(3)

由式(2)、 (3)可知:歸一化互相關(guān)成像條件是在互相關(guān)成像條件的基礎(chǔ)上,用源點電磁波場或接收點電磁波場進行歸一化,可對成像結(jié)果整體上起到能量調(diào)解的作用。該成像條件具有削弱源點振幅影響、提升深部結(jié)構(gòu)能量等優(yōu)點。

2 基于時間窗互相關(guān)成像條件的GPR逆時偏移算法

由于實際地下介質(zhì)的速度存在差異,各層介質(zhì)之間會產(chǎn)生相互反射和背景雜波干擾,利用互相關(guān)成像條件進行成像時,各成像點的正傳和反傳電磁波場既包含目標(biāo)反射信號,也包含介質(zhì)之間相互反射引起的雜波干擾[15]。這些干擾信號的互相關(guān)會降低目標(biāo)體的成像質(zhì)量。因此,可根據(jù)偏移速度模型估計每個成像點的主反射信號產(chǎn)生時間,設(shè)計合理的時間窗以壓制時間窗外的雜波干擾[22]。對于偏移速度模型中的某一個成像點(x,z), GPR信號主反射發(fā)生時間tk(x,z)可表示為

tk(x,z)=tmax-d(x,z)/vk,

(4)

其中:d(x,z)為成像點(x,z)到發(fā)射天線的距離;vk為電磁波在地下第k種介質(zhì)中的傳播速度(k=1, 2, …,K);tmax為時窗長度。 假設(shè)整個探測空間被第k種均勻介質(zhì)填充, 則成像點(x,z)處產(chǎn)生的主反射時間可估計為d(x,z)/vk。 由于電磁波在實際介質(zhì)中的傳播速度在[min(vk) max(vk)]區(qū)內(nèi)波動,故成像點(x,z)處的主反射產(chǎn)生時間介于[min(d(x,z)/vk) max(d(x,z)/vk)]。

min(d(x,z)/vk)≤tk(x,z)≤max(d(x,z)/vk)。

(5)

為壓制該成像點附近處的雜波干擾,該點的時間窗上、下界限可定義為

tm1(x,z)=max(min(tk(x,z)-tw),0),

(6)

tm2(x,z)=min(max(tk(x,z)),tmax)。

(7)

其中：tw為GPR發(fā)射信號的脈沖寬度,k=1, 2, …,K。 在時間窗[tm1tm2]內(nèi), 該成像點的正傳和反傳電磁波場歸一化互相關(guān)計算公式可表示為

(8)

其中, TG(x,z)將作為成像點(x,z)處反射強度的加權(quán)系數(shù)以抑制時間窗外的雜波干擾。因此,成像點(x,z)處的最終逆時偏移結(jié)果可表示為

NRTTG(x,z)=NRT(x,z)·TG(x,z),

(9)

其中，NRT(x,z)表示該成像點處的歸一化互相關(guān)成像結(jié)果。 由式(9)可知,基于時間窗函數(shù)的GPR逆時偏移算法主要包括逆時偏移的實現(xiàn)和權(quán)值系數(shù)TG(x,z)的計算。 對于逆時偏移的實現(xiàn)過程: 首先, 利用鉆孔取樣、 層析成像、 全波形反演等技術(shù)獲取實際地下介質(zhì)的電磁波速度參數(shù), 建立逆時偏移所需的速度模型； 然后,計算正傳電磁波場S(x,z,t)和反傳電磁波場R(x,z,t)； 最后，應(yīng)用歸一化互相關(guān)成像條件獲取逆時偏移結(jié)果。

對于權(quán)值系數(shù)TG(x,z)的計算過程：根據(jù)式(6)、 (7)計算偏移速度模型各成像點的主反射時間及時間窗上、 下界限。利用式(8)求解每個成像點處的權(quán)值系數(shù)。 圖1展示了某成像點處正傳和反傳電磁波形, 正傳電磁波形到達目標(biāo)的時刻為目標(biāo)的主反射信號發(fā)射時刻,tmax為雷達的最大采樣時間。 在獲得歸一化互相關(guān)成像條件的GPR逆時偏移結(jié)果和權(quán)值系數(shù)之后, 利用式(9),計算基于時間窗互相關(guān)成像條件的GPR逆時偏移剖面。

基于時間窗互相關(guān)成像條件的GPR逆時偏移成像流程如圖2所示,其中內(nèi)虛線框中為逆時偏移所需步驟,外虛線框中為基于時間窗函數(shù)的GPR逆時偏移所需步驟[22]。

圖1 逆時偏移算法中某成像點處入射信號和反向信號傳播示意圖Fig.1 Timing of transmitted and reflected signals in RTM of GPR

圖2 基于時間窗互相關(guān)成像條件的GPR逆時偏移成像流程圖Fig.2 GPR reverse time migration imaging flow chart based on time window cross-correlation imaging condition

3 數(shù)值算例

3.1 模型1

模型1是大小為1.8 m×1.0 m的圓形空洞模型, 如圖3所示。 背景介質(zhì)設(shè)置為干土壤,其相對介電常數(shù)為5, 電導(dǎo)率為0.001 S/m; 土壤中埋有一個半徑為0.1 m的圓形空洞, 其圓心位置為(0.9 m, 0.5 m)。 采用二維FDTD正演算法計算正傳與反傳電磁波場, 計算區(qū)域被剖分成360×200的正方形網(wǎng)格單元, 網(wǎng)格間距為0.005 m×0.005 m; 模型外施加完全匹配層, 厚度為0.1 m。 采用共源觀測方式, 發(fā)射源信號采用中心頻率為500 MHz的雷克子波, 采樣時間間隔0.01 ns, 時窗長度16 ns。 地表布設(shè)5個發(fā)射源(三角形表示), 發(fā)射源的間距為0.2 m; 第1個發(fā)射源放置在水平距離0.5 m處; 每個發(fā)射源兩側(cè)各布置50個接收源,接收點的間距為0.01 m,水平方向上覆蓋了0.5 m。

圖4為(0.55 m, 0.4 m)處的源點正傳電磁波和接收點反傳電磁波波形對比, 兩條豎直虛線分別表示時間窗上、 下界限。 時間窗[4.12 ns 8.34 ns]內(nèi)的正傳電磁波波形與反傳電磁波波形具有很強的相關(guān)性,在此時間窗內(nèi)進行互相關(guān)計算,可以壓制其他時刻的正傳和反傳電磁波場互相關(guān)運算產(chǎn)生的噪聲干擾。

分別利用基于時間窗互相關(guān)成像條件和歸一化互相關(guān)成像條件在第3個發(fā)射源點的偏移孔徑內(nèi)進行偏移成像獲得的成像剖面如圖5所示。 逆時偏移剖面中的圓形空洞均清晰可見, 驗證了本文編制的基于歸一化互相關(guān)成像條件和時間窗互相關(guān)成像條件的GPR逆時偏移算法的正確性和有效性。 對比圖5a、 c 可知, 拉普拉斯濾波前的歸一化互相關(guān)成像結(jié)果在空洞附近出現(xiàn)較強的背景噪聲, 且地表附近受發(fā)射源強振幅影響范圍更大; 而時間窗互相關(guān)成像條件的逆時偏移剖面中, 由于采用時間窗內(nèi)的歸一化互相關(guān)運算結(jié)果作為權(quán)值系數(shù), 壓制了目標(biāo)體附近的噪聲和發(fā)射源強振幅對空洞目標(biāo)體的影響。 圖5b、 d 為利用拉普拉斯濾波方法對圖5a、 c進行濾波后獲得的成像剖面。 對比濾波前后的逆時偏移剖面可知, 濾波可壓制低頻噪聲和發(fā)射源強振幅的影響, 提高目標(biāo)體的成像精度。 對比圖5b、d 可知, 濾波后的時間窗互相關(guān)成像條件的GPR逆時偏移剖面的分辨率更高。

圖3 圓形空洞模型示意圖Fig.3 Schematic map of the circle void model

圖4 x=0.55 m, y=0.4 m處的源點正傳電磁波(a)與接收點反傳電磁波(b)波形Fig.4 Forward and reverses electromagnetic wave forms at position of x=0.55 m, y=0.4 m

為更好地說明時間窗互相關(guān)成像條件的優(yōu)越性, 分別提取圖5中地表正中心位置處的單道波形能量進行對比, 圖6a、 b分別為拉普拉斯濾波前、 濾波后的歸一化互相關(guān)成像條件和時間窗互相關(guān)成像條件的GPR逆時偏移波形對比。 雖然時間窗互相關(guān)成像條件的逆時偏移結(jié)果中波形能量更弱些, 但波形的持續(xù)長度更小, 分辨率和精度更高。

圖7為5個發(fā)射源點的逆時偏移疊加成像剖面，圓狀空洞在逆時偏移剖面都清晰可見,且與空洞的真實位置相符。對比圖7a、c可知,基于時間窗互相關(guān)成像條件的GPR逆時偏移算法在一定程度上壓制了空洞附近的噪聲,偏移剖面的低頻噪聲更弱。從對比濾波前后的逆時偏移剖面(圖7b、d)可知,濾波后的時間窗互相關(guān)成像條件的GPR逆時偏移剖面中的空洞目標(biāo)體成像結(jié)果更接近真實情況,精度更好,分辨率更高。

圖5 模型1中第3個發(fā)射源的GPR逆時偏移剖面Fig.5 RTM profiles of Model 1 in the migration aperture of the 3th transmitting source with normalized cross correlation imaging condition and time gating cross correlation imaging conditiona—歸一化互相關(guān)成像條件; b—歸一化互相關(guān)成像條件+濾波; c—時間窗互相關(guān)成像條件; d—時間窗互相關(guān)成像條件+濾波

圖6 圖5中正中心位置(x=0.9 m, y=0.0 m)的單道波形能量對比Fig.6 Comparison of single waveform at center position of x=0.9 m, y=0.0 m in Fig.5

圖7 GPR模型1的逆時偏移成像剖面Fig.7 RTM profiles of GPR Model 1 with normalized cross correlation imaging condition and time gating cross correlation imaging conditiona—歸一化互相關(guān)成像條件;b—歸一化互相關(guān)成像條件+濾波;c—時間窗互相關(guān)成像條件;d—時間窗互相關(guān)成像條件+濾波

3.2 模型2

為更好地分析時間窗互相關(guān)成像條件在背景為隨機介質(zhì)中的目標(biāo)體的成像效果,分別建立三圓非隨機介質(zhì)和隨機介質(zhì)GPR模型,如圖8所示。模型大小為2.0 m×1.0 m,分上、下兩層介質(zhì),上層介質(zhì)與模型1相同,埋深為0.4 m;下層介質(zhì)為濕沙層,其相對介電常數(shù)為10,電導(dǎo)率為0.002 S/m,濕沙層中埋有3個半徑為0.06 m的圓形空洞,其圓心位置(圖8a)分別為(0.5m, 0.7m)、 (1.0m, 0.8m)、 (1.5m, 0.9m) 。在圖 8a基礎(chǔ)上,利用隨機過程中的譜分析和橢圓自相關(guān)函數(shù)理論,建立了相應(yīng)的隨機介質(zhì)模型, 上、 下層介質(zhì)的擾動分別是0.3和0.5, 如圖 8b所示。 利用二維FDTD正演算法計算電磁波在該模型中的正向及逆向傳播, 計算區(qū)域被剖分成400×200的正方形網(wǎng)格單元, 網(wǎng)格間距為0.005 m×0.005 m, 模型外完全匹配層厚度為0.1 m。 采用共源觀測方式,發(fā)射源信號采用中心頻率為500 MHz的雷克子波, 采樣時間間隔為0.01 ns, 時窗長度為20 ns。 在地表布設(shè)20個發(fā)射源,發(fā)射源的間距為0.1 m;第1個反射源放置在水平距離0 m處;每個發(fā)射源布置50個接收源,接收點間距為0.01 m,水平方向上覆蓋了0.5 m。

圖8 GPR模型2示意圖Fig.8 Diagrams of GPR Model 2a—非隨機介質(zhì);b—隨機介質(zhì)

圖9為非隨機介質(zhì)模型2的逆時偏移剖面圖,土壤與濕沙的分界面及3個空洞的成像清晰可見,且與真實空間位置相符。相比于濾波前的逆時偏移剖面,濾波后成像剖面的水平分界面和空洞更清晰,發(fā)射源強振幅的影響得到壓制,剖面分辨率更高。對比圖8b、d可知,基于時間窗互相關(guān)成像條件的GPR逆時偏移,由于采用時間窗內(nèi)的歸一化互相關(guān)結(jié)果作為權(quán)值系數(shù),壓制了分界面和3個空洞附近的雜波和低頻噪聲,偏移剖面中的低頻噪聲更微弱,分辨率更高。

圖9 非隨機介質(zhì)模型2的GPR逆時偏移成像剖面Fig.9 RTM profiles for nonrandom medium of Model 2 with normalized cross correlation and time gating cross correlation imaging conditiona—歸一化互相關(guān)成像條件;b—歸一化互相關(guān)成像條件+濾波;c—時間窗互相關(guān)成像條件;d—時間窗互相關(guān)成像條件+濾波

圖10是分別利用歸一化互相關(guān)成像條件和時間窗互相關(guān)成像條件對隨機介質(zhì)模型2進行逆時偏移, 獲得的成像剖面。 可見, 土壤和濕沙的分界面和3個空洞的依然清晰可見, 但由于背景介質(zhì)分布的隨機性, 水平界面成像有微小擾動, 空洞異常體的形狀出現(xiàn)扭曲現(xiàn)象。基于歸一化互相關(guān)成像條件的逆時偏移結(jié)果中背景噪聲非常強烈, 對空洞成像造成了嚴(yán)重干擾(圖10a); 基于時間窗互相關(guān)成像條件的逆時偏移剖面中大部分背景噪聲得到壓制,成像精度和分辨率得到極大提高(圖10d)。

圖10 隨機介質(zhì)模型2的GPR逆時偏移成像剖面Fig.10 RTM profiles for random medium of Model 2 with normalized cross correlation and time gating cross correlation imaging conditiona—歸一化互相關(guān)成像條件; b—歸一化互相關(guān)成像條件+濾波; c—時間窗互相關(guān)成像條件; d—時間窗互相關(guān)成像條件+濾波

4 實測數(shù)據(jù)的應(yīng)用

為驗證基于時間窗互相關(guān)成像條件的逆時偏移算法應(yīng)用于實測GPR數(shù)據(jù)處理中的效果,開展地下管線GPR探測實驗, 獲得實測原始GPR剖面如圖11a所示: 原始剖面中管線產(chǎn)生的雙曲線反射波能量微弱, 難以識別及提取。 利用零時校正、 帶通濾波、 時變增益等常規(guī)處理方法進行處理后,獲得的剖面如圖11b所示: 經(jīng)過常規(guī)處理后的雷達剖面, 雙曲線反射波更清晰, 能量更強, 噪聲更低。為后續(xù)逆時偏移處理提供了高質(zhì)量的GPR數(shù)據(jù)。分別利用歸一化互相關(guān)成像條件和時間窗互相關(guān)成像條件對圖11b所示的雷達剖面進行處理,獲得的成像剖面見圖11c、d:經(jīng)過逆時偏移后的雙曲線反射波能量收斂,反射波歸位,能精確確定管線的埋深約為1.0 m,水平位置約為5.5 m。對比兩圖可知:相比于圖11c,圖11d中基于時間窗互相關(guān)成像條件的逆時偏移剖面的干擾波得到壓制,成像精度和分辨率得到極大提高。

圖11 實測數(shù)據(jù)的GPR剖面及其逆時偏移成像剖面Fig.11 GPR profiles and RTM profiles of observed dataa—原始GPR剖面; b—常規(guī)處理后的GPR剖面; c—歸一化互相關(guān)成像條件的逆時偏移剖面; d—時間窗互相關(guān)成像條件的逆時偏移剖面

5 結(jié) 論

從逆時偏移原理出發(fā),根據(jù)偏移模型的速度分布估算每個成像點的反射時間,設(shè)計了一個時間窗函數(shù)以壓制時間窗外的雜波干擾,提出了一種基于時間窗互相關(guān)成像條件的GPR逆時偏移算法。兩個典型GPR模型和實測數(shù)據(jù)的逆時偏移成像結(jié)果表明:相比于歸一化互相關(guān)成像條件的逆時偏移結(jié)果,基于時間窗互相關(guān)成像條件具有能明顯壓制背景低頻噪聲的優(yōu)點,成像精度更好,空間分辨率更高。