王敏玲， 龔俊波， 王洪華， 羅崇炎

(1. 桂林理工大學(xué) 地球科學(xué)學(xué)院，桂林 541004;2. 廣西隱伏金屬礦產(chǎn)勘查重點(diǎn)實(shí)驗(yàn)室，桂林 541004)

0 引言

探地雷達(dá)(Ground Penetrating Radar, GPR)是利用高頻脈沖電磁波來探測地下地質(zhì)體結(jié)構(gòu)及物性參數(shù)的一種重要淺部地球物理方法[1]，具有無損性、高分辨率、高效率等特點(diǎn)，在工程勘查、環(huán)境監(jiān)測、考古調(diào)查等領(lǐng)域中得到廣泛應(yīng)用[2-6]。逆時(shí)偏移成像[7]作為一種將實(shí)測雷達(dá)數(shù)據(jù)與真實(shí)地質(zhì)結(jié)構(gòu)轉(zhuǎn)換的關(guān)鍵技術(shù)，一直以來是GPR研究的重點(diǎn)，通過逆時(shí)偏移成像可獲得地下精細(xì)結(jié)構(gòu)分布，提高解譯精度[8]。

對GPR逆時(shí)偏移成像的研究，國內(nèi)、外學(xué)者已經(jīng)開展過相關(guān)工作。底青云[9]考慮地下介質(zhì)電導(dǎo)率對GPR信號的衰減，開展了衰減介質(zhì)GPR逆時(shí)偏移成像研究；魯興林和錢榮毅[10]通過數(shù)值算例分析了逆時(shí)偏移與克?；舴蚱频膬?yōu)缺點(diǎn)，得出了逆時(shí)偏移成像精度更高、效果更好的結(jié)論；Fu等[11]通過對傳統(tǒng)逆時(shí)偏移成像剖面進(jìn)行去模糊濾波處理，有效提高了逆時(shí)偏移的成像精度；Liu等[12]提出了一種基于分層介質(zhì)格林函數(shù)的頻率域GPR逆時(shí)偏移算法，提高了計(jì)算效率；王敏玲等[13]將激發(fā)振幅成像條件運(yùn)用于GPR逆時(shí)偏移成像中，有效提高了逆時(shí)偏移計(jì)算效率；Bradford等[14]為克服地形起伏對成像精度的影響，提出了一種適用于起伏地表的GPR逆時(shí)偏移算法，數(shù)值算例表明：該算法比傳統(tǒng)高程靜校正偏移方法的成像精度更好。

上述GPR逆時(shí)偏移成像研究大都將偏移速度模型看作均勻、分塊介質(zhì)模型，側(cè)重于理論數(shù)據(jù)的逆時(shí)偏移成像研究。然而實(shí)際地下介質(zhì)是隨機(jī)分布的，電磁波在其中傳播時(shí)會產(chǎn)生大量的不相干擾動，其實(shí)質(zhì)是源于介質(zhì)在小尺度上的非均勻性[15]。采用確定性的均勻介質(zhì)模型難以描述小尺度上的非均勻、隨機(jī)分布特性。因此，如何真實(shí)描述地下介質(zhì)的小尺度非均勻特性，成為當(dāng)前GPR研究的熱點(diǎn)。隨機(jī)介質(zhì)模型建立理論是將地下介質(zhì)的小尺度非均勻性性視為一種隨機(jī)過程，采用統(tǒng)計(jì)學(xué)方法加以描述介質(zhì)小尺度上的非均勻分布特征[16]。戴前偉等[17]采用隨機(jī)過程中的譜分解理論和橢圓自相關(guān)函數(shù)構(gòu)建隨機(jī)介質(zhì)GPR模型，并利用無單元法分析電磁信號在隨機(jī)介質(zhì)的傳播特征；宋二喬等[18]利用時(shí)域有限差分法模擬了雷達(dá)波在季節(jié)性凍土隨機(jī)介質(zhì)中的傳播特征；郭士禮等[19]采用多相隨機(jī)介質(zhì)模型建立方法構(gòu)建了混凝土隨機(jī)分布結(jié)構(gòu)，并分析了GPR信號的隨機(jī)擾動特征與多相離散隨機(jī)介質(zhì)模型參數(shù)之間的關(guān)系。筆者在上述理論研究基礎(chǔ)上，采用隨機(jī)過程中的譜分解理論和橢圓自相關(guān)函數(shù)構(gòu)建了2個(gè)典型的GPR隨機(jī)介質(zhì)模型，并采用時(shí)域有限差分法(Finite Difference Time Domain, FDTD)正演模擬和逆時(shí)偏移成像，分析了隨機(jī)介質(zhì)中的電磁波傳播特征和異常體的逆時(shí)偏移成像效果，為可否利用GPR方法實(shí)現(xiàn)地質(zhì)結(jié)構(gòu)的高精度探測提供理論參考。

1 隨機(jī)介質(zhì)GPR模型建立方法

隨機(jī)介質(zhì)模型主要由非均勻性大、小兩種尺度組成。大尺度非均勻性描述介質(zhì)的背景特性，小尺度非均勻性是加載在大尺度上的一種隨機(jī)擾動，通常用一個(gè)均值為零的二階平穩(wěn)隨機(jī)過程來表示[16]。GPR探測過程中，電磁波在地下介質(zhì)中傳播時(shí)主要受介質(zhì)相對介電常數(shù)的影響。二維隨機(jī)介質(zhì)中空間坐標(biāo)(x,z)處的相對介電常數(shù)ε(x,z)可表示為：

ε(x,z)=ε0(x,z)+δε(x,z)

(1)

其中：ε0(x,z)為背景介質(zhì)(大尺度)的相對介電常數(shù)；δε(x,z)為背景介質(zhì)上的空間相對隨機(jī)擾動量(小尺度)可表示為：

φ(x,z)=δε(x,z)/ε0

(2)

假設(shè)φ是具有零均值、一定方差及自相關(guān)函數(shù)的空間平穩(wěn)隨機(jī)過程[16]，則φ滿足[17]：

〈φ(x,z)〉=0

(3)

〈φ2(x,z)〉=α2

(4)

構(gòu)造二階平穩(wěn)過程ε(x,z)的步驟如下[16, 20]：

1) 選擇自相關(guān)函數(shù)。目前，描述隨機(jī)介質(zhì)的自相關(guān)函數(shù)有多種。其中，混合型自相關(guān)函數(shù)能較好地描述具有多尺度平滑的隨機(jī)介質(zhì)而在彈性波、電磁波正演模擬領(lǐng)域得到廣泛應(yīng)用，其表達(dá)式為：

(5)

式中：a,b分別是介質(zhì)在x，z方向上的自相關(guān)長度，是表示相似程度大小的參數(shù)；p為粗糙度因子。

2) 計(jì)算混合型自相關(guān)函數(shù)φ(x,z)的二維傅里葉變換φ(kx,kz)。

3) 用隨機(jī)數(shù)發(fā)生器生成0～2π之間服從獨(dú)立均值分布的二維隨機(jī)場φ(kx,kz)。

4) 計(jì)算隨機(jī)功率譜函數(shù)ψ(kx,kz)，其表達(dá)式為：

exp[iθ(kx,kz)]

(6)

式中，ω(kx,kz)為窗函數(shù)，通過窗函數(shù)的應(yīng)用可降低空間離散誤差，從而減少自相關(guān)函數(shù)的低頻能量[17]。

5) 計(jì)算ψ(kx,kz)的二維傅里葉逆變換，從而得到隨機(jī)擾動v(x,z)。

6) 計(jì)算隨機(jī)擾動v(x,z)的均值和方差：

μ=〈v(x,z)〉

(7)

d2=〈[v(x,z)-μ]2〉

(8)

7) 通過規(guī)范化生成均值為零、方差為α2的隨機(jī)擾動：

(9)

圖1為按照上述理論，構(gòu)建的不同自相關(guān)長度隨機(jī)相對介電常數(shù)模型。其中，背景介質(zhì)的相對介電常數(shù)為3，其隨機(jī)擾動量的方差為0.1。由圖1可知：隨機(jī)介質(zhì)中的相對介電常數(shù)呈隨機(jī)分布，局部區(qū)域存在微小異常；自相關(guān)長度a,b分別表示隨機(jī)介質(zhì)在x,z方向上擾動的平均尺度，自相關(guān)長度越大，則該方向上相對介電常數(shù)的連續(xù)性越好；當(dāng)a=∞，表示在x方向均勻分布，此時(shí)二維隨機(jī)介質(zhì)變成一維層狀隨機(jī)介質(zhì)，如圖1(d)所示。由此可見，筆者采用的隨機(jī)介質(zhì)建立方法能靈活、有效地描述地下實(shí)際介質(zhì)的非均勻隨機(jī)分布特性。

2 GPR逆時(shí)偏移成像方法

GPR逆時(shí)偏移[21]主要包括以下三個(gè)過程：源點(diǎn)電磁波場的正向傳播，接收點(diǎn)電磁波場的逆時(shí)外推和成像條件的應(yīng)用。前兩個(gè)過程的實(shí)現(xiàn)可利用相同的數(shù)值模擬方法，如FDTD[22]、時(shí)域有限單元法[23]、時(shí)域偽譜法[24]等，其中FDTD具有原理簡單、計(jì)算效率高、易實(shí)現(xiàn)等優(yōu)點(diǎn)而在GPR數(shù)值模擬中得到廣泛應(yīng)用。因此，我們采用FDTD計(jì)算源點(diǎn)正傳電磁波場和接收點(diǎn)反傳電磁波場。

圖1 不同相對常數(shù)的橢圓自相關(guān)函數(shù)產(chǎn)生的隨機(jī)介質(zhì)相對介電常數(shù)模型Fig.1 Random media models of relative dielectric permittivity in different auto-correlation length media produced with intermixed ellipsoidal autocorrelation function(a)a=1,b=1;(b)a=10,b=10;(c)a=100,b=20;(d)a=∞,b=20

成像條件的應(yīng)用是影響逆時(shí)偏移成像效果的另一個(gè)重要因素。目前，GPR逆時(shí)偏移中應(yīng)用最廣泛的成像條件是互相關(guān)成像條件[25]?；ハ嚓P(guān)成像條件以時(shí)間一致性原理為基礎(chǔ)，通過計(jì)算源點(diǎn)正傳電磁波場和接收點(diǎn)反傳電磁波場的零延時(shí)互相關(guān)，從而獲得成像剖面。標(biāo)準(zhǔn)的零延遲互相關(guān)成像條件公式為：

(10)

為此，Kaelin等[26]采用源點(diǎn)正傳電磁波場或接收點(diǎn)反傳電磁波場對零延時(shí)互相關(guān)成像條件的計(jì)算結(jié)果進(jìn)行歸一化，提出了歸一化互相關(guān)成像條件，其表達(dá)式為：

(11)

由式(11)可知：歸一化互相關(guān)成像條件不但能削弱近地表源點(diǎn)強(qiáng)振幅影響，而且可提升深部反射界面的能量。目前，歸一化互相關(guān)成像條件在彈性波、聲波和電磁波逆時(shí)偏移中得到廣泛應(yīng)用[27, 28, 21]。筆者采用該成像條件開展GPR逆時(shí)偏移計(jì)算，開發(fā)的隨機(jī)介質(zhì)GPR逆時(shí)偏移成像算法如圖2所示。逆時(shí)偏移的第一步是給定偏移模型和激勵(lì)源，利用FDTD模擬激勵(lì)源在模型中傳播的正傳波場，計(jì)算至最大采樣時(shí)間tmax；第二步是將接收點(diǎn)所接受到的數(shù)據(jù)借助FDTD沿時(shí)間軸進(jìn)行電磁波場的逆時(shí)延拓至零時(shí)刻；第三步，把相同時(shí)刻相同位置上的電磁波場應(yīng)用相關(guān)成像條件進(jìn)行成像，本文用歸一化互相關(guān)成像條件進(jìn)行偏移成像。對單個(gè)源點(diǎn)進(jìn)行歸一化互相關(guān)成像后，需要對偏移結(jié)果進(jìn)行拉普拉斯濾波去除剖面頂部的低頻噪聲，并將多個(gè)源點(diǎn)偏移結(jié)果進(jìn)行疊加，獲得最終的偏移結(jié)果。

圖2 隨機(jī)介質(zhì)GPR逆時(shí)偏移流程圖Fig.2 Flow chart of prestack RTM algorithm in random media

圖3 圓狀體GPR模型及其逆時(shí)偏移結(jié)果Fig.3 Schematic diagram of circle model and its migration reverse time result(a)圓柱體GPR模型;(b)多偏移距GPR正演剖面; (c)逆時(shí)偏移剖面;(d)單道逆時(shí)偏移結(jié)果

3 數(shù)值算例

依據(jù)上述方法原理，編制了相應(yīng)的隨機(jī)介質(zhì)GPR疊前逆時(shí)偏移程序，利用該程序?qū)蓚€(gè)典型隨機(jī)介質(zhì)GPR模型的多偏移距正演數(shù)據(jù)進(jìn)行逆時(shí)偏移成像，并與相應(yīng)背景介質(zhì)的GPR模型的逆時(shí)偏移結(jié)果進(jìn)行對比分析。

圖4 不同自相關(guān)長度的隨機(jī)介質(zhì)圓柱體模型Fig.4 Random media circle models of relative dielectric permittivity with different auto-correlation length(a)a=1,b=1;(b)a=1,b=10; (c)a=15,b=15;(d)a=15,b=30

3.1 模型一

圖3(a)是大小為1.8 m×1.8 m的圓形空洞模型示意圖，模型介質(zhì)的相對介電常數(shù)為5，電導(dǎo)率為0.001 S/m；中間埋有一個(gè)半徑為0.1 m的圓形空洞，其圓心位置為(0.9 m, 0.9 m)。利用FDTD法對該模型進(jìn)行多偏移距正演和逆時(shí)偏移時(shí)，發(fā)射源信號采用中心頻率為500 MHz的雷克子波，采樣時(shí)間間隔為0.01 ns，時(shí)窗長度為25 ns；地表布設(shè)15個(gè)發(fā)射天線，發(fā)射天線的間距為0.1 m，第一個(gè)發(fā)射天線放置在原點(diǎn)，其余發(fā)射天線依次往右布置；每個(gè)發(fā)射天線右側(cè)布置50道接收天線，接收點(diǎn)的間距為0.01 m，水平方向上覆蓋了0.5 m。圖3(b)為模型一的多偏移距GPR正演剖面，由圖可知：圓形空洞產(chǎn)生的雙曲線反射波清晰可見，但不能準(zhǔn)確定位異常體的真實(shí)空間位置。利用二維GPR逆時(shí)偏移成像算法對圖3(b)中的正演數(shù)據(jù)進(jìn)行逆時(shí)偏移，獲得了成像剖面如圖3(c)所示。偏移剖面中，繞射波能量幾乎都得到收斂，反射波能量歸位到空洞的真實(shí)空間位置，且空洞頂部和底部真實(shí)位置與圖3(c)中水平位置0.9 m處的歸一化單道波形能量位置相符，如圖3(d)，所示驗(yàn)證了本文構(gòu)建的GPR逆時(shí)偏移算法的正確性和有效性。

為了更好地分析二維GPR逆時(shí)偏移成像算法對實(shí)際地下隨機(jī)分布介質(zhì)的成像能力，將模型一的背景介質(zhì)設(shè)置為不同自相關(guān)長度的隨機(jī)介質(zhì)，建立的模型a、模型b、模型c、模型d,如圖4所示。由圖4可見，隨機(jī)介質(zhì)中自相關(guān)長度是影響介質(zhì)非均勻尺度的重要因素，自相關(guān)長度越大，非均勻尺度越大，自相關(guān)長度越小，非均勻尺度越小，取不同的自相關(guān)長度可以很好地表示地下不同非均勻尺度的隨機(jī)介質(zhì)模型。

圖5為利用FDTD對圖4所示的隨機(jī)介質(zhì)模型進(jìn)行多偏移距正演計(jì)算，獲得GPR正演剖面。由圖5可見：自相關(guān)長度a、b均為1時(shí)，圓形空洞產(chǎn)生的反射波清晰，如圖5(a)所示；自相關(guān)長度越大，雙曲線反射波能量越微弱，如圖5(b)和圖5(c)所示；當(dāng)自相關(guān)長度a=15,b=30,如圖5(d)所示，雙曲線反射波形態(tài)失真，不易被識別。對比圖5和圖3(b)可知：電磁波在隨機(jī)介質(zhì)中傳播時(shí)散射強(qiáng)烈，散射波相互疊加，形成了明顯的隨機(jī)擾動，致使反射波扭曲變形，不連續(xù)，不易辨識，降低了GPR回波的信噪比和分辨率。

圖5 圖4中的隨機(jī)介質(zhì)GPR模型的多偏移距正演剖面Fig.5 Simulated multi-offset GPR profiles of media GPR models with different autocorrelation length as in Fig 4(a)a=1,b=1;(b)a=1,b=10; (c)a=15,b=15;(d)a=15,b=30

圖6 利用圖5中的多偏移距GPR正演數(shù)據(jù)進(jìn)行逆時(shí)偏移獲得的成像剖面Fig.6 Imaging results of simulated multi-offset GPR dataset by using the GPR reverse time migration of random media(a)a=1,b=1;(b)a=1,b=10; (c)a=15,b=15;(d)a=15,b=30

圖7 圖6中地表中心位置(0.9 m)處的單道波形對比Fig.7 Comparison of single waveforms extracted from the center surface position of 0.9 m in Fig 6(a)0 m～1.8 m;(b) 0.6 m～1.2 m

圖8 復(fù)雜GPR模型二示意圖Fig.8 Schematic diagram of complex GPR model(a)背景介質(zhì)為均勻介質(zhì);(b)背景介質(zhì)為隨機(jī)介質(zhì)

利用隨機(jī)介質(zhì)GPR疊前逆時(shí)偏移成像算法對圖5中的多偏移距正演數(shù)據(jù)進(jìn)行逆時(shí)偏移，獲得的成像剖面如圖6所示。由圖可見，自相關(guān)長度越小，繞射波歸位到其真實(shí)空間位置的效果越好，能量收斂越集中，圓形空洞成像越清晰；自相關(guān)長度越大，繞射波歸位效果越差，圓形空洞成像越模糊；當(dāng)自相關(guān)長度達(dá)到15以上時(shí)，反射波和繞射波能量幾乎不能收斂，異常體的空間分布位置難以被識別。

圖9 復(fù)雜GPR模型二多偏移距正演剖面圖Fig.9 Simulated multi-offset GPR model of complex GPR model II(a)背景介質(zhì)為均勻介質(zhì);(b)背景介質(zhì)為隨機(jī)介質(zhì)

圖10 利用圖9中的多偏移距正演數(shù)據(jù)進(jìn)行逆時(shí)偏移獲得的成像剖面Fig.10 Imaging results of simulated multi-offset GPR dataset by using the GPR reverse time migration of random media(a)背景介質(zhì)為均勻介質(zhì);(b)背景介質(zhì)為隨機(jī)介質(zhì)

圖7(a)為從圖6地表中心0.9 m處提取的歸一化單道波形能量對比圖，對比圖3(d)和圖7(a)可知，隨機(jī)介質(zhì)中的歸一化波形能量連續(xù)出現(xiàn)多個(gè)波峰與波谷，散射波能量收斂差，異常體的空間分布位置難以被識別。而圖3(d)中的散射波能量幾乎都得到收斂。為了更好地對比空洞位置附近的波形能量收斂性，截取了0.6 m～1.2 m范圍內(nèi)的波形對比如圖7(b)所示。由圖7可知：自相關(guān)長度越小，波形能量越強(qiáng)，且繞射波能量越收斂于空洞的真實(shí)位置(如紅色虛線表示)處。

3.2 模型二

模型二是大小為2.0 m×2.0 m的復(fù)雜GPR模型示意圖，如圖8(a)所示。模型被一條起伏界面分成上、下兩層，起伏界面的埋深約0.6 m；上層介質(zhì)的相對介電常數(shù)為5，電導(dǎo)率為0.02 S/m；下層介質(zhì)的相對介電常數(shù)為10，電導(dǎo)率為0.001 S/m。下層介質(zhì)的左右兩邊分別埋有一個(gè)圓形空洞異常體，其半徑分別為0.15 m和0.05m，圓心位置分別為(0.50 m, 1.2 m)和 (1.60 m, 0.9 m)。圖8(b)為圖8(a)相應(yīng)的隨機(jī)介質(zhì)模型示意圖，其上層介質(zhì)自相關(guān)長度a=1,b=10，方差為0.1，其下層介質(zhì)自相關(guān)長度a=1,b=1，方差為0.1。利用FDTD和逆時(shí)偏移成像算法對模型二進(jìn)行計(jì)算時(shí)的參數(shù)與計(jì)算模型一時(shí)相同。

圖9為利用FDTD法對模型二進(jìn)行多偏移距正演計(jì)算獲得的GPR正演剖面。由圖9(a)可知，當(dāng)背景為均勻介質(zhì)時(shí)，電磁波未發(fā)生散射，反射波形規(guī)則、有序，起伏界面清晰可見；當(dāng)背景介質(zhì)為隨機(jī)介質(zhì)時(shí)，電磁波散射強(qiáng)烈、隨機(jī)無序傳播，散射波相互疊加，如圖9(b)所示。由于上層隨機(jī)介質(zhì)的自相關(guān)長度大，電磁波散射強(qiáng)烈，使得起伏界面的反射波能量異常微弱，無法被識別；而下層隨機(jī)介質(zhì)的自相關(guān)長度小，電磁波散射較弱，兩個(gè)圓形空洞的反射波能量較強(qiáng)，能被識別。

利用隨機(jī)介質(zhì)GPR疊前逆時(shí)偏移算法對圖9中的多偏移距正演數(shù)據(jù)進(jìn)行逆時(shí)偏移，獲得的剖面如圖10所示。由圖10可見，當(dāng)背景介質(zhì)為均勻介質(zhì)時(shí)，起伏界面、圓形空洞的成像清晰，反射波和繞射波得到收斂，且歸位到其真實(shí)空間位置，如圖10(a)所示；當(dāng)背景介質(zhì)為隨機(jī)介質(zhì)時(shí)，由于電磁波在其中傳播時(shí)發(fā)生散射，起伏界面的反射波和圓形空洞的雙曲線反射波形態(tài)失真、扭曲，起伏界面和圓形空洞的反射波能量幾乎不能收斂，難以精確定位其空間位置。

通過圖10可知：由于介質(zhì)分布的隨機(jī)性，電磁波在其中傳播時(shí)散射強(qiáng)烈，波形發(fā)生扭曲，并產(chǎn)生隨機(jī)擾動使得電磁波在反傳過程中不能精確到達(dá)真實(shí)空間位置，從而造成異常體的逆時(shí)偏移成像結(jié)果比背景介質(zhì)為均勻介質(zhì)的逆時(shí)偏移結(jié)果差；自相關(guān)長度是影響隨機(jī)介質(zhì)中異常體成像效果的主要因素，自相關(guān)長度越小，異常體的成像越清晰、準(zhǔn)確；自相關(guān)長度越大，異常體的成像效果越差，且不易被識別。

4 結(jié)論

1) 闡述了隨機(jī)過程的譜分解和混合型自相關(guān)函數(shù)理論構(gòu)建GPR隨機(jī)介質(zhì)模型方法及其具體步驟；利用FDTD法和歸一化互相關(guān)成像條件構(gòu)建了GPR隨機(jī)介質(zhì)疊前逆時(shí)偏移成像算法。

2) 不同自相關(guān)長度的GPR隨機(jī)介質(zhì)模型的逆時(shí)偏移剖面對比表明：自相關(guān)長度是影響異常體成像效果的主要因素；與背景介質(zhì)為均勻介質(zhì)的GPR模型的逆時(shí)偏移結(jié)果相比，隨機(jī)介質(zhì)的逆時(shí)偏移成像結(jié)果的空間分辨率更低，低頻噪聲更強(qiáng)。