王洪華， 龔俊波， 梁值歡， 張智， 徐濤

1 桂林理工大學， 地球科學學院， 廣西 桂林 541004 2 中國科學院地質(zhì)與地球物理研究所， 巖石圈演化國家重點實驗室， 北京 100029 3 中國科學院地球科學研究院， 北京 100029

0 引言

探地雷達(Ground Penetrating Radar, GPR)作為一種高頻脈沖電磁探測技術(shù)，以其效率高、分辨率高、無損探測、實時成像等優(yōu)點，而被廣泛應用于工程檢測和淺層勘探領(lǐng)域(盧成明等, 2007; Mochales et al., 2008; 劉瀾波和錢榮毅, 2015; Atef et al., 2016; 廖紅建等, 2016; 郭士禮等, 2019).近年來，工程檢測和淺層勘探的日益細化，給GPR數(shù)據(jù)的快速高精度處理與成像帶來極大挑戰(zhàn)(程久龍等, 2010; 蘇茂鑫等, 2010; 黃忠來和張建中, 2013; 王敏玲等, 2019a).為此，許多學者根據(jù)電磁波與彈性波傳播規(guī)律的相似性和GPR與地震數(shù)據(jù)接收方式的相似性，將一些成熟的地震數(shù)據(jù)處理與成像方法如動校正疊加(Ebihara et al., 2000; Perroud and Tygel, 2004)、速度譜分析(Grandjean et al., 2000; Booth et al., 2011)、層析成像(Johnson et al., 2007; Chang and Alumbaugh, 2011)、繞射疊加(Feng and Sato, 2004; Aitken and Stewart, 2004)、克?；舴蚱?Moran et al., 2000; Porsani and Sauck, 2007)、逆時偏移(Fisher et al., 1992; Leuschen and Plumb, 2001; Bradford, et al., 2018)等引入到GPR數(shù)據(jù)處理與成像中.其中，逆時偏移因其具有原理簡單、計算效率快、成像精度高等優(yōu)點，而被廣泛應用于淺部精細結(jié)構(gòu)的GPR高精度成像中，取得了良好的成像效果(傅磊等, 2014; Liu et al., 2014,2016,2017; Bradford, 2015; Lu et al., 2016;王敏玲等, 2019b; Zhang et al., 2019).

然而，電磁波與彈性波在地下介質(zhì)傳播規(guī)律存在顯著差異：二階彈性波動方程只涉及位移關(guān)于時間的二次導數(shù)(波動項)，然而地下介質(zhì)電導率的存在，使得二階電磁波方程還涉及電磁場關(guān)于時間的一次導數(shù)(衰減項)(Di and Wang, 2004).與彈性波在地下介質(zhì)中傳播時只表現(xiàn)波動特性相比，電磁波既有波動特性，也有表現(xiàn)地下介質(zhì)吸收的衰減特性.高頻電磁波在地下介質(zhì)中傳播時，傳播速度和衰減系數(shù)是關(guān)于介質(zhì)電導率的函數(shù)，電導率越高，傳播速度越小，衰減系數(shù)越大，能量更易衰減(Bergmann et al., 1998; Neto and Mediros, 2006; 張先武等, 2014; 王洪華等, 2018).因此，在GPR逆時偏移中，考慮地下高電導率介質(zhì)對電磁波的強吸收衰減作用，對提高高電導率區(qū)域的成像質(zhì)量顯得尤為必要.常規(guī)GPR逆時偏移在計算反傳電磁波場時大都沿用正傳電磁波方程，正傳電磁波場在高電導率介質(zhì)中衰減的同時，反傳電磁波場會再次衰減，能量非常微弱，難以實現(xiàn)高電導率區(qū)域的清晰準確成像(朱尉強和黃清華, 2016; 王敏玲等, 2019a).如何在計算反傳電磁波場的同時，對電磁波正傳時衰減的能量進行精確補償，以提高高電導率區(qū)域的成像質(zhì)量是目前該領(lǐng)域的研究熱點.Sena等(2006)在裂步-傅里葉偏移算法的基礎(chǔ)上，通過在頻率域電磁波場外推的同時進行反濾波處理，補償衰減的電磁波場，有效提高了高衰減區(qū)域的成像質(zhì)量.其后，Oden等(2007)將上述方法應用于GPR頻率-波數(shù)偏移算法中，數(shù)值試驗論證了該算法的有效性.然而，反濾波方法大都基于一維衰減模型，難以適用于復雜地質(zhì)結(jié)構(gòu)的GPR逆時偏移成像(Zhu et al., 2016; 朱尉強和黃清華, 2016).近年來，針對黏彈性介質(zhì)高精度成像提出的衰減補償逆時偏移算法為高電導率介質(zhì)區(qū)域的GPR高精度成像提供了一種可行有效的方案(Zhu, 2014; Zhu and Harris, 2014).該算法通過在彈性波場反傳過程中，人為改變黏彈性波動方程中衰減項的正負號，以保持逆時外推的時間對稱性和反傳不變性，精確補償正傳時衰減的彈性波場能量，提高黏彈性介質(zhì)的成像質(zhì)量(Zhu et al., 2014; Zhu and Harris, 2015; Zhu, 2016).目前，該算法在黏彈性介質(zhì)高精度逆時偏移中得到廣泛應用(李振春等, 2014; Sun et al., 2016; 田坤等, 2017; 劉財?shù)? 2018; 豆輝和徐逸鶴, 2019)，取得了良好的成像效果.從數(shù)學上看，二階電磁波方程與二階黏彈性波動方程形式類似，都涉及衰減項，既表現(xiàn)了波動特性，也表現(xiàn)了介質(zhì)的吸收衰減特性.為此，Zhu等(2016)、朱尉強和黃清華(2016)分別根據(jù)兩者的相似性，成功將黏彈性波衰減補償?shù)哪鏁r偏移算法應用于二維高電導率介質(zhì)結(jié)構(gòu)的GPR成像中，詳細推導補償電磁波正傳時衰減能量的反傳電磁波方程，并用數(shù)值試驗論證了該算法應用于提高高電導率介質(zhì)區(qū)域的成像效果的可行性和有效性.

考慮到實際GPR高頻電磁波是在地下三維空間輻射傳播，二維逆時偏移難以實現(xiàn)反射波的準確歸位和繞射波的完全收斂，成像精度降低(Liu et al., 2017, 2018; 張崇明等, 2019; Zhu et al., 2020).本文在Zhu等(2016)、朱尉強和黃清華(2016)的基礎(chǔ)上，開展基于電磁波衰減補償?shù)娜SGPR逆時偏移成像研究.其中，三維時域有限差分法用于計算正傳和反傳電磁波場，并通過改變反傳電磁波方程中衰減項的正負號，以補償電磁波正傳時衰減的能量；零時刻成像條件用于獲得三維逆時偏移結(jié)果.數(shù)值試驗論證了本文構(gòu)建的基于電磁波衰減補償?shù)娜SGPR逆時偏移算法在高電導率區(qū)域成像分辨率和抗干擾能力方面的優(yōu)勢.

1 方法原理

1.1 電磁波的衰減特征

根據(jù)電磁波場理論，忽略激勵源的影響，三維GPR電磁波方程可表示為(馮德山等，2017)：

(1)

假定電場為時諧場，式(1)兩邊都進行傅里葉變換并整理，可推導電磁波復傳播速度為(Carcione, 2014)：

(2)

式中，i為虛數(shù)單位，將式(2)展開，可推導電磁波在地下介質(zhì)中傳播的速度和衰減系數(shù)表達式為(Zhu et al., 2016)：

(3)

(4)

其中，sgn(x)為符號函數(shù)，當x>0時，sgn(x)=1；當x<0時sgn(x)=-1.

由式(3)和(4)可知，電磁波速度和衰減系數(shù)是關(guān)于介質(zhì)電導率的函數(shù).圖1為均勻介質(zhì)(相對介電常數(shù)為6，頻率為400 MHz)中電磁波速度和衰減系數(shù)隨電導率變化曲線，其中實線和虛線分別是電導率取正值和負值所得.由圖可知，與電導率為0.0001 S·m-1時的電磁波速度相比，電導率為0.01 S·m-1時的電磁波速度變化約為0.1%，受電導率變化影響較??；當電導率從0.0001 S·m-1增大到0.01 S·m-1時，衰減系數(shù)從0.008增大到0.8，受電導率變化影響較大.由此可見：電導率是影響電磁波能量衰減的關(guān)鍵參數(shù)，特別是在高電導率區(qū)域中電磁波能量衰減更強.因此，對在高電導率區(qū)域采集的GPR數(shù)據(jù)進行逆時偏移時，補償電磁波衰減的能量顯得尤為必要.

1.2 逆時偏移方法

GPR逆時偏移原理是在構(gòu)建偏移速度模型的基礎(chǔ)上，將實測GPR信號作為邊界條件在時間軸上進行逆時外推，當逆推至零時刻時應用相關(guān)成像條件獲取成像結(jié)果，從而實現(xiàn)地下精細結(jié)構(gòu)的高精度成像(王敏玲等, 2019b).根據(jù)時間反轉(zhuǎn)原理(Fink, 1992; Zhu et al.,2016; 朱尉強和黃清華, 2016)，電磁波場進行逆時外推滿足方程為

圖1 電磁波速度(a)和衰減系數(shù)(b)隨電導率變化曲線Fig.1 Curves of velocity (a) and attenuation coefficient (b) of electromagnetic waves varying with conductivity

(5)

(6)

為避免電磁波逆時外推時的衰減，補償正傳過程中電磁波衰減的能量，Zhu等(2016)提出了一種改變式(5)中衰減項前的正負號方法，人為保持電磁波場逆時外推的反轉(zhuǎn)不變性，即：

(7)

式(7)與式(1)的形式一致，可精確補償電磁波正傳中衰減的能量(朱尉強和黃清華, 2016).

本文采用三維時域有限差分法進行接收點電磁波場的逆時外推，零時刻成像條件用于獲取成像結(jié)果.當電磁波逆時外推至零時刻時，零時刻的電磁波場即為地下結(jié)構(gòu)的成像結(jié)果(王敏玲等, 2019b)，可表示為：

(8)

2 數(shù)值計算

2.1 方法的可行性和有效性驗證

為驗證本文提出的三維衰減補償電磁波場逆時外推方法的可行性和有效性，建立了一個1.5 m×1.5 m×1.5 m三維均勻模型，其相對介電常數(shù)εr=8.三維FDTD用于模擬計算時的空間步長均為0.01 m，時間步長為0.015 ns，時間長度為24 ns；激勵源是中心頻率為400 MHz的雷克子波.首先，將激勵源放置于模型的正中心(0.75 m, 0.75 m,0.75 m)，分別將均勻模型的電導率σ設置為 0 S·m-1(無損)、0.001 S·m-1、0.015 S·m-1、-0.015 S·m-1，獲得的7 ns時刻Ey分量的波場快照，如圖2所示.由圖可見，介質(zhì)電導率越大，電磁波能量衰減更強、能量越弱.當電導率σ為-0.015 S·m-1時，即利用式(7)進行模擬計算時，保持了時間反轉(zhuǎn)不變性和時間對稱性，電磁波衰減得到有效補償，與圖2c中的電磁波能量相比，能量被有效恢復；且與圖2a無損情況下的電磁波能量相當，說明通過改變電磁波方程中衰減項前的正負號，可有效補償電磁波在高電導率介質(zhì)中衰減的能量.

圖2 不同電導率均勻模型中7 ns時刻Ey分量的波場快照(a) 0 S·m-1； (b) 0.001 S·m-1； (c) 0.015 S·m-1； (d) -0.015 S·m-1.Fig.2 Snapshots of the Ey wave field of homogenous model at 7ns with different conductivity values

圖3為不同電磁波場逆時外推方法在均勻模型正中心位置處接收到的波形對比，其中灰實線為三維FDTD正演在模型正中心接收到的波形；黑點虛線、黑實線和黑虛線分別是模型電導率為0 S·m-1、0.015 S·m-1、-0.015 S·m-1時電磁波場逆時外推接收到的波形，即將模型最外層所有網(wǎng)格點作為接收點接收到的GPR信號進行逆時外推后在模型中心位置處接收的波形.由圖3可知：常規(guī)不考慮電導率的逆時偏移無法對電磁波衰減進行補償，如黑點虛線所示；常規(guī)考慮電導率的逆時偏移比不考慮電導率的電磁波能量衰減更強，成像結(jié)果更差；而本文構(gòu)建的基于電磁波衰減補償?shù)哪鏁r偏移結(jié)果與正演波形相比，能較好地補償由電導率引起的電磁波衰減如黑色虛線所示，驗證了本文構(gòu)建的三維衰減補償電磁波場逆時外推方法的可行性和有效性.

圖3 不同逆時外推電磁波場重構(gòu)方法在均勻模型正中心位置處接收到的波形對比灰實線為正演接收到的波形，黑點虛線、黑實線和黑虛線分別為電導率為0 S·m-1 (無損)、0.015 S·m-1、-0.015 S·m-1時逆時外推接收到的波形.Fig.3 Comparison of reconstructed waveforms at center position of homogenous model by using different reverse time extrapolation methodsThe grey line is the simulated waveform, black dot-dashed line, black line and black dotted line are the reconstructed waves of reverse time extrapolation by using the homogenous model with conductivity of 0 S·m-1, 0.015 S·m-1, and -0.015 S·m-1, respectively.

2.2 空洞模型

為驗證本文構(gòu)建的基于電磁波衰減補償?shù)娜SGPR逆時偏移方法的成像效果，建立了一個大小為0.8 m×2 m×1.1 m的空洞模型，如圖4所示.模型被埋深為0.5 m水平界面分為上下兩層，其相對介電常數(shù)分別為6和8；下層介質(zhì)的左右兩邊分別埋有一個大小為0.1 m×0.1 m×0.1 m的正方體空洞，其中心分別位于(0.4 m,0.5 m,0.8 m)、(0.4 m,1.5 m,0.8 m)，如圖4a所示.模型的背景電導率為0.001 S·m-1，右側(cè)設置了一個大小為0.4 m×0.4 m×0.8 m高電導率區(qū)域，其電導率為0.015 S·m-1，中心位置為(0.4 m,1.5 m,0.7 m)，如圖4b所示.利用三維FDTD進行模擬計算時的參數(shù)與均勻介質(zhì)模型相同，平行Y方向X=0～0.8 m之間等距布設9條測線，測線間距為0.1 m，收發(fā)天線間距為0.06 m；平行X方向Y=0～2 m之間等距布設11條測線，測線間距為0.2 m，收發(fā)天線間距為0.06 m.

圖5a、b分別為空洞模型X方向和Y方向上三維GPR正演切片.由圖5a可見：8 ns處，Y=1.5 m附近出現(xiàn)較強的反射波，這是由于高電導率區(qū)域與背景電導率差異明顯，反射系數(shù)不為零所致.11 ns附近出現(xiàn)上、下層介質(zhì)分界面產(chǎn)生的水平反射波，波形能量強、易識別；受高電導率區(qū)域(Y=1.5 m附近)的影響，電磁波在傳播過程中出現(xiàn)較為明顯地衰減，波形能量較弱，如X=0.3 m、0.4 m、0.5 m位置處的正演切片所示.16 ns開始出現(xiàn)空洞產(chǎn)生的雙曲線繞射波，空洞正上方測線(X=0.4 m)的正演切片中繞射波能量最強，其他測線上的正演切片中繞射波能量隨距離增大而變?nèi)酢⒊霈F(xiàn)時間變長.低電導率區(qū)域(Y=0.5 m附近)中空洞產(chǎn)生的繞射波能量比高電導率區(qū)域(Y=1.5 m附近)能量更強、波形更明顯，這是由于電磁波能量在高電導率區(qū)域衰減更強所致.分析圖5b中Y方向上的正演切片可得到類似的結(jié)論.

圖4 空洞模型的示意圖(a) 相對介電常數(shù)分布； (b) 電導率分布.Fig.4 Schematic diagrams of void GPR model(a) Relative permittivity； (b) Conductivity.

圖5 空洞模型的三維正演剖面(a) X方向； (b) Y方向.Fig.5 3D GPR forward profile of void model(a) X direction； (b) Y direction.

利用本文構(gòu)建的基于電磁波衰減補償?shù)娜SGPR逆時偏移算法對圖5所示的三維正演剖面進行逆時偏移成像，并與常規(guī)逆時偏移和介質(zhì)無損情況下的逆時偏移結(jié)果進行對比，獲得結(jié)果如圖6所示.由圖6可知：三種逆時偏移成像剖面中水平界面產(chǎn)生的反射波能量得到準確歸位，空洞產(chǎn)生的繞射波完全收斂，成像結(jié)果清晰準確.但三種逆時偏移方法對高電導率區(qū)域的成像分辨率存在明顯差別：圖6a展示的常規(guī)三維GPR逆時偏移結(jié)果中，由于未考慮電導率對電磁波能量衰減的影響，高電導率區(qū)域處的水平界面與空洞位置處的成像非常模糊、不易被識別；這是由于電磁波在高電導率區(qū)域進行逆時外推時能量再次衰減所致.與圖6a相比，圖6b所示的基于電磁波衰減補償?shù)腉PR三維逆時偏移結(jié)果中高電導率區(qū)域中衰減的電磁波能量得到較好補償，水平界面和空洞的成像能量得到較好的恢復，成像結(jié)果更清晰、準確；且與介質(zhì)無損(電導率為0)情況下的三維GPR逆時偏移結(jié)果圖6c吻合較好.

為更好地分析基于電磁波衰減補償?shù)娜SGPR逆時偏移成像方法對高電導率區(qū)成像的優(yōu)勢，提取X=0.4 m、Y=1.5 m的單道波形對比，如圖7所示.由圖7可見，三種GPR逆時偏移成像結(jié)果中，水平界面和空洞成像位置與真實位置相符；相比介質(zhì)無損情況下三維GPR逆時偏移結(jié)果中的波形能量，常規(guī)GPR逆時偏移結(jié)果波形能量衰減約86%；基于電磁波衰減補償?shù)娜SGPR逆時偏移結(jié)果中的能量得到有效恢復，且與介質(zhì)無損情況下的逆時偏移結(jié)果基本吻合.由此可見：基于電磁波衰減補償?shù)娜SGPR逆時偏移可有效補償電磁波在高電導率介質(zhì)中傳播損失的能量，大大提升了目標體的成像精度和分辨率，其結(jié)果更有利于后續(xù)雷達資料的解釋.

圖6 空洞模型三維GPR正演數(shù)據(jù)的逆時偏移剖面(a) 常規(guī)逆時偏移結(jié)果； (b) 衰減補償逆時偏移結(jié)果； (c) 介質(zhì)無損情況下的逆時偏移結(jié)果.Fig.6 RTM imaging results of 3D forward GPR profile of void model(a) Conventional RTM； (b) Aattenuation compensated RTM； (c) Conventional RTM in lossless media.

2.3 分層界面模型

圖8是一個大小為0.8 m×2 m×1.5 m的分層界面模型，從上至下分為四層，各層的介電參數(shù)和幾何參數(shù)分布如圖所示，在第二層介質(zhì)左邊存在一個局部高電導率區(qū)域，其電導率為0.012 S·m-1.利用三維FDTD對該模型進行模擬的激勵源是中心頻率為600 MHz的雷克子波，其余計算參數(shù)與空洞模型相同.

利用基于電磁波衰減補償?shù)娜SGPR逆時偏移算法對該模型正演數(shù)據(jù)進行逆時偏移成像，并與常規(guī)GPR逆時偏移結(jié)果和介質(zhì)無損情況下的逆時偏移結(jié)果進行對比，如圖9所示.由圖9可見：常規(guī)逆時偏移結(jié)果中由于電磁波在高電導率區(qū)域傳播時衰減較強，附近分界面的成像模糊、能量微弱，特別是最下層的分界面難以被識別；相比常規(guī)逆時偏移結(jié)果，基于電磁波衰減補償?shù)哪鏁r偏移結(jié)果中，高電導率區(qū)域附近的反射界面的成像能量得到有效補償且與介質(zhì)無損情況下的逆時偏移結(jié)果相當，成像結(jié)果清晰可見，且易被識別.

圖7 圖6中X=0.4 m、 Y=1.5 m的單道波形對比Fig.7 Comparison of single-trace waveforms at position of X=0.4 m, Y=1.5 m in Fig.6

圖8 分層界面模型示意圖(a) 相對介電常數(shù)分布； (b) 電導率分布.Fig.8 Schematic diagrams of layered interface GPR model(a) Relative permittivity； (b) Conductivity.

圖9 分層界面模型三維GPR正演數(shù)據(jù)的逆時偏移剖面(a) 常規(guī)逆時偏移結(jié)果； (b) 衰減補償逆時偏移結(jié)果； (c) 介質(zhì)無損情況下的逆時偏移結(jié)果.Fig.9 RTM imaging results of the 3D forward GPR profile of layered interface model(a) Conventional RTM algorithm； (b) Attenuation compensated RTM algorithm； (c) The result conventional RTM with lossless media.

為驗證基于電磁波衰減補償?shù)娜SGPR逆時偏移算法的抗干擾能力，對X=0.4 m處的正演剖面圖10a分別施加50%、70%和90%噪聲后的GPR剖面如圖10b、c、d所示.由圖可知，噪聲越強，界面產(chǎn)生的反射波越難以被識別，施加90%噪聲后的GPR剖面中，第三層界面的反射波難以被識別.

圖11是對施加噪聲的GPR剖面進行三維GPR衰減補償逆時偏移和常規(guī)逆時偏移結(jié)果對比，其中圖11 a、c、e是常規(guī)逆時偏移剖面，圖11 b、d、f是衰減補償逆時偏移結(jié)果.由圖可知：施加的噪聲越強，成像結(jié)果越模糊，雜波干擾越強；但不同程度噪聲干擾下的GPR剖面中的反射波均得到準確歸位，且與其真實位置相符，也易被識別.特別是施加90%噪聲情況下，部分有效波已經(jīng)被嚴重污染，但通過逆時偏移仍然能對其進行較好成像.對比常規(guī)GPR逆時偏移和衰減補償GPR逆時偏移結(jié)果可知，本文構(gòu)建的基于電磁波衰減補償?shù)娜SGPR逆時偏移成像算法具有較強的抗干擾能力.值得一提是，基于電磁波衰減補償?shù)哪鏁r偏移在對反射波能量進行補償?shù)耐瑫r，也存在對噪聲信號進行補償?shù)默F(xiàn)象.

3 結(jié)論

(1) 本文從電磁波的衰減特性和逆時偏移原理出發(fā)，通過改變?nèi)S反傳電磁波方程中包含電導率的衰減項的正負號，保持電磁波反傳的時間對稱性和不變性，以精確補償電磁波在正傳中衰減的能量，構(gòu)建了一種基于電磁波衰減補償?shù)娜SGPR逆時偏移算法.其中，三維FDTD用于計算正傳和反傳電磁波場，零時刻成像條件用于獲取地下介質(zhì)的成像結(jié)果.

(2) 兩個典型三維GPR模型的正演剖面的電磁波衰減補償逆時偏移和常規(guī)逆時偏移結(jié)果對比表明：本文構(gòu)建的基于電磁波衰減補償?shù)娜SGPR逆時偏移算法可精確補償電磁波在地下介質(zhì)傳播時衰減的能量，高電導率區(qū)域的成像分辨率更高，抗干擾能力更強，其結(jié)果更有利于指導后續(xù)雷達剖面的解譯.

圖10 施加不同比例噪聲的GPR剖面(a) 未施加噪聲； (b) 50%； (c) 70%； (d) 90%.Fig.10 GPR profile added with noises of different proportions(a) Without noise; (b) 50%; (c) 70%; (d) 90%.

圖11 圖10中GPR剖面的衰減補償逆時偏移和常規(guī)逆時偏移結(jié)果對比(a)、(c)和(e)是常規(guī)逆時偏移結(jié)果； (b)、(d)和(f)是衰減補償逆時偏移結(jié)果.Fig.11 Comparison of attenuation compensated RTM and conventional RTM of the GPR profile shown in Fig.10(a), (c) and (e) Conventional RTM；(b), (d) and (f) Attenuation compensated RTM.