基于時間窗互相關成像條件的GPR逆時偏移成像

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

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

0 引 言

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

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

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

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

1 GPR逆時偏移成像原理

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

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

(1)

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

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

(2)

(3)

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

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

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

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

(4)

其中:d(x,z)為成像點(x,z)到發射天線的距離;vk為電磁波在地下第k種介質中的傳播速度(k=1, 2, …,K);tmax為時窗長度。 假設整個探測空間被第k種均勻介質填充, 則成像點(x,z)處產生的主反射時間可估計為d(x,z)/vk。 由于電磁波在實際介質中的傳播速度在[min(vk) max(vk)]區內波動,故成像點(x,z)處的主反射產生時間介于[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發射信號的脈沖寬度,k=1, 2, …,K。 在時間窗[tm1tm2]內, 該成像點的正傳和反傳電磁波場歸一化互相關計算公式可表示為

(8)

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

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

(9)

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

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

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

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

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

3 數值算例

3.1 模型1

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

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

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

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

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

圖5 模型1中第3個發射源的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—歸一化互相關成像條件; b—歸一化互相關成像條件+濾波; c—時間窗互相關成像條件; d—時間窗互相關成像條件+濾波

圖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—歸一化互相關成像條件;b—歸一化互相關成像條件+濾波;c—時間窗互相關成像條件;d—時間窗互相關成像條件+濾波

3.2 模型2

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

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

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

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

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

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

4 實測數據的應用

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

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

5 結 論

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