隨機介質中的探地雷達疊前逆時偏移成像

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

(1. 桂林理工大學 地球科學學院，桂林 541004;2. 廣西隱伏金屬礦產勘查重點實驗室，桂林 541004)

0 引言

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

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

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

1 隨機介質GPR模型建立方法

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

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

(1)

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

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

(2)

假設φ是具有零均值、一定方差及自相關函數的空間平穩隨機過程[16]，則φ滿足[17]：

〈φ(x,z)〉=0

(3)

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

(4)

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

1) 選擇自相關函數。目前，描述隨機介質的自相關函數有多種。其中，混合型自相關函數能較好地描述具有多尺度平滑的隨機介質而在彈性波、電磁波正演模擬領域得到廣泛應用，其表達式為：

(5)

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

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

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

4) 計算隨機功率譜函數ψ(kx,kz)，其表達式為：

exp[iθ(kx,kz)]

(6)

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

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

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

μ=〈v(x,z)〉

(7)

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

(8)

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

(9)

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

2 GPR逆時偏移成像方法

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

圖1 不同相對常數的橢圓自相關函數產生的隨機介質相對介電常數模型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

成像條件的應用是影響逆時偏移成像效果的另一個重要因素。目前，GPR逆時偏移中應用最廣泛的成像條件是互相關成像條件[25]?；ハ嚓P成像條件以時間一致性原理為基礎，通過計算源點正傳電磁波場和接收點反傳電磁波場的零延時互相關，從而獲得成像剖面。標準的零延遲互相關成像條件公式為：

(10)

為此，Kaelin等[26]采用源點正傳電磁波場或接收點反傳電磁波場對零延時互相關成像條件的計算結果進行歸一化，提出了歸一化互相關成像條件，其表達式為：

(11)

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

圖2 隨機介質GPR逆時偏移流程圖Fig.2 Flow chart of prestack RTM algorithm in random media

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

3 數值算例

依據上述方法原理，編制了相應的隨機介質GPR疊前逆時偏移程序，利用該程序對兩個典型隨機介質GPR模型的多偏移距正演數據進行逆時偏移成像，并與相應背景介質的GPR模型的逆時偏移結果進行對比分析。

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

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

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

圖5 圖4中的隨機介質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正演數據進行逆時偏移獲得的成像剖面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 復雜GPR模型二示意圖Fig.8 Schematic diagram of complex GPR model(a)背景介質為均勻介質;(b)背景介質為隨機介質

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

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

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

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

3.2 模型二

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

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

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

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

4 結論

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

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