探地雷達(dá)頻率域2.5維正演

戴世坤 歐陽(yáng)振崇* 周印明 張錢江 李 昆趙東東 陳輕蕊 凌嘉宣

①(中南大學(xué)地球科學(xué)與信息物理學(xué)院 長(zhǎng)沙 410083)

②(中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室 長(zhǎng)沙 410083)

③(東方地球物理公司綜合物化探處 涿州 072751)

④(桂林理工大學(xué)地球科學(xué)學(xué)院 桂林 541006)

1 引言

探地雷達(dá)(Ground Penetrating Radar, GPR)是一種利用高頻(106～109Hz)電磁波來(lái)確定介質(zhì)內(nèi)部物質(zhì)分布規(guī)律的地球物理方法，因其具有抗干擾能力強(qiáng)、探測(cè)分辨率高、場(chǎng)地適應(yīng)能力強(qiáng)、操作簡(jiǎn)單且為無(wú)損探測(cè)等優(yōu)勢(shì)而被廣泛地應(yīng)用于工程勘察及地質(zhì)調(diào)查。GPR正演是雷達(dá)理論研究的重點(diǎn)之一，隨著勘探要求更加精細(xì)化和高效，大范圍高效、高精度的數(shù)值模擬和反演成像成為現(xiàn)在GPR技術(shù)的目標(biāo)[1,2]。

波動(dòng)方程正演模擬方法，因其能包含雷達(dá)波的運(yùn)動(dòng)學(xué)特征和動(dòng)力學(xué)特征被廣泛應(yīng)用于GPR正演模擬中，主要包括時(shí)域有限差分法(Finite-Difference Time Domain method, FDTD)和有限單元法(Finite Element Method, FEM)兩種。在FDTD法的應(yīng)用方面，自1996年Yee[3]提出著名的Yee氏網(wǎng)格后，F(xiàn)DTD法被廣泛應(yīng)用于GPR數(shù)值模擬中。Bergmann等人[4]應(yīng)用FDTD法開(kāi)展了不均勻非線性和衰減介質(zhì)的GPR正演模擬；Irving等人[5]采用PML吸收邊界研究了TE和TM模式的2維GPR數(shù)值模擬；劉四新等人[6]對(duì)比了GPR2維頻散介質(zhì)與非頻散介質(zhì)中GPR信號(hào)的區(qū)別；馮德山等人[7]研究了FDTD數(shù)值模擬中不分裂卷積完全匹配層對(duì)倏逝波的吸收效果。

FDTD法原理簡(jiǎn)單，易于編程實(shí)現(xiàn)，但要求模型規(guī)則剖分，對(duì)復(fù)雜問(wèn)題適應(yīng)性差。FEM因具有網(wǎng)格剖分靈活，適用于物性參數(shù)分布復(fù)雜和幾何特征不規(guī)則模型的優(yōu)勢(shì)，被引入到GPR數(shù)值模擬領(lǐng)域。底青云等人[8]推導(dǎo)GPR有限元方程，開(kāi)展了一系列典型模型的正演模擬；杜華坤等人[9]基于優(yōu)化的Delaunay三角剖分，采用線性插值進(jìn)行了FEM的GPR2維正演，提高了FEM對(duì)復(fù)雜模型模擬的適應(yīng)性和數(shù)值模擬精度；石明等人[10]采用Delaunay非結(jié)構(gòu)化網(wǎng)格有限元法進(jìn)行了各向異性介質(zhì)GPR有限元正演；王洪華等人[11]實(shí)現(xiàn)了PML邊界條件在2階電磁波動(dòng)方程GPR時(shí)域有限元模擬中的應(yīng)用，驗(yàn)證了PML邊界條件在復(fù)雜地電結(jié)構(gòu)電磁波傳播模擬中具有良好的吸收效果。

近年來(lái)GPR數(shù)值模擬大都在時(shí)間域內(nèi)研究，時(shí)間域GPR波的傳播滿足的波動(dòng)方程模擬結(jié)果符合雷達(dá)波傳播的運(yùn)動(dòng)學(xué)特征，但在表現(xiàn)波的動(dòng)力學(xué)特征方面存在不足，文獻(xiàn)[12]指出頻率域波形反演在地震中的重要地位，首次研究了頻率域波形反演，對(duì)頻域波形反演存在的問(wèn)題進(jìn)行了初次探討，本文從頻率域出發(fā)，研究了2.5維GPR數(shù)值模擬方法，這樣可以很好地保留雷達(dá)波傳播的運(yùn)動(dòng)學(xué)特征和動(dòng)力學(xué)特征，準(zhǔn)確地研究雷達(dá)波在頻率域的傳播特性，為GPR全波形反演提供重要基礎(chǔ)。本文采用有限單元法，推導(dǎo)基于行波分解的吸收邊界條件的GPR有限元方程，實(shí)現(xiàn)頻率域2.5D探地雷達(dá)正演模擬。重點(diǎn)分析和總結(jié)在雷達(dá)頻率不同相對(duì)介電常數(shù)和不同收發(fā)距波數(shù)選取的規(guī)律；設(shè)計(jì)均勻全空間和半空間模型，將數(shù)值解與解析解對(duì)比驗(yàn)證了算法的正確性；另外，雷達(dá)2.5維頻率域正演在不同波數(shù)之間的計(jì)算具有高度并行性，通過(guò)設(shè)計(jì)Open MP并行，探究不同線程下算法并行的效率，驗(yàn)證了算法的高效性。

2 方法理論

2.1 控制方程

2.2 有限單元法

3 數(shù)值算例

本文算法的程序代碼采用Fortran語(yǔ)言編寫(xiě)，測(cè)試平臺(tái)配置為CPU-Inter(R) Core(TM) i9-7980XE，主頻2.60 GHz(18核，36線程)，內(nèi)存64 GB,64位操作系統(tǒng)。

3.1 波場(chǎng)參數(shù)分析

圖1 dx=0.1 m處不同相對(duì)介電常數(shù)的電磁場(chǎng)譜隨波數(shù)變化特征

圖2 dx=0.5 m處不同相對(duì)介電常數(shù)電磁場(chǎng)譜隨波數(shù)變化特征

圖3 dx=1 m處不同相對(duì)介電常數(shù)電磁場(chǎng)譜隨波數(shù)變化特征

圖4 dx=5 m處不同相對(duì)介電常數(shù)電磁場(chǎng)譜隨波數(shù)變化特征

圖5 dx=10 m處不同相對(duì)介電常數(shù)電磁場(chǎng)譜隨波數(shù)變化特征

結(jié)論：頻率為100 MHz，計(jì)算距離范圍為10 m時(shí)，波數(shù)最大值選取103即可將譜的能量全包含其中，而且譜的能量主要分布在波數(shù)0～10的范圍內(nèi)，可對(duì)該范圍內(nèi)波數(shù)適當(dāng)加密，使傅里葉逆變換精度更高。

3.2 正確性驗(yàn)證

3.2.1 全空間模型

設(shè)計(jì)均勻全空間模型，σ =0.001 S/m , εr=1.0，節(jié)點(diǎn)301×301，電流1000 A，偶極距dl為1 m，源中心在原點(diǎn)，x方向源網(wǎng)格均0.02 m，源以外網(wǎng)格采用非均勻剖分最大0.05 m，邊界取15個(gè)網(wǎng)格間距由0.05 m以1.5倍遞增至1 m。x, z方向網(wǎng)格剖分相同，f=100 MHz，正波數(shù)范圍(0.01, 1000)，波數(shù)55個(gè)。

圖6中主剖面上電磁場(chǎng)數(shù)值解與解析解的形態(tài)、數(shù)值都相同，主剖面測(cè)線z=0.5 m處各節(jié)點(diǎn)的相對(duì)誤差均小于1%，驗(yàn)證了本文算法的正確性。

圖7為采用文獻(xiàn)[15]中的算法計(jì)算主剖面電磁場(chǎng)的數(shù)值解與解析解對(duì)比圖，由圖可知，文獻(xiàn)[15]算法在源附近的誤差較大，表明本文算法在源的處理上要優(yōu)于文獻(xiàn)[15]中的算法，計(jì)算精度更高。

3.2.2 半空間模型

圖6 主剖面電磁場(chǎng)數(shù)值解與解析解對(duì)比圖

圖7 主剖面電磁場(chǎng)其他算法數(shù)值解與解析解對(duì)比圖

設(shè)計(jì)均勻半空間模型，空氣層電導(dǎo)率σ =10?12S/m ，地下介質(zhì)電導(dǎo)率σ =10?3S/m，相對(duì)介電常數(shù) εr=1.0，頻率f=100 MHz，網(wǎng)格節(jié)點(diǎn)601×601，x方向的源布設(shè)于地下0.5 m，源的其他參數(shù)和網(wǎng)格剖分與全空間相同，選取正波數(shù)范圍(0.0001, 1000)，波數(shù)277個(gè)。

圖8為主剖面(y = 0 m)數(shù)值解與解析解對(duì)比，第3列為剖面測(cè)線z=0.5 m的相對(duì)誤差，由圖可知，數(shù)值解與解析解的形態(tài)、數(shù)量級(jí)相同，測(cè)線各節(jié)點(diǎn)的相對(duì)誤差均小于1.5%，同樣驗(yàn)證了本文算法的準(zhǔn)確性。

3.3 Open MP并行測(cè)試

本文算法耗時(shí)主要在波數(shù)循環(huán)計(jì)算上，每個(gè)波數(shù)均需求解1次方程組，但各波數(shù)計(jì)算相互獨(dú)立，可采用并行方式提高計(jì)算效率。目前電磁法2.5D正反演應(yīng)用較多的并行方法有MPI(Message Passing Interface)和OpenMP(Open Multi-Processing)。OpenMP使用線程間共享內(nèi)存的方式協(xié)調(diào)并行計(jì)算，對(duì)原串行代碼改動(dòng)小，容易實(shí)現(xiàn)。本文采用OpenMP并行處理不同波數(shù)方程組的求解、傅里葉逆變換和輔助場(chǎng)計(jì)算。

2.5D GPR正演計(jì)算量大，將Pardiso求解器采用OpenMP并行求解，算法效率如表1。模型參數(shù)與3.2.1小節(jié)模型一致，網(wǎng)格節(jié)點(diǎn)301×301，波數(shù)277個(gè)。式中加速比 = 單機(jī)的單線程耗時(shí)/并行計(jì)算耗時(shí)。

表1中，隨著并行線程個(gè)數(shù)增加，加速比增大，算法耗時(shí)減少，改造后的計(jì)算效率明顯提高。2線程耗時(shí)比單線程耗時(shí)長(zhǎng)，可能原因是2線程時(shí)發(fā)生了同步事件，耗時(shí)變長(zhǎng)[16]；16線程耗時(shí)約80 s，比單線程計(jì)算快了3倍；增加到20線程時(shí)，耗時(shí)反倒增大，分析原因是當(dāng)計(jì)算量一定時(shí)，線程數(shù)量增加，用于線程通信/線程調(diào)度的時(shí)間所占比例逐漸增大，計(jì)算效率降低。

結(jié)論：采用OpenMP將Pardiso求解器并行，隨著并行線程個(gè)數(shù)增加，計(jì)算效率提高，但并行線程數(shù)并非越多越好，綜合不同線程算法耗時(shí)和計(jì)算機(jī)資源的占用情況，本節(jié)線程數(shù)為16時(shí)，加速比最大，耗時(shí)最短，是當(dāng)前條件下OpenMP并行選取的最佳線程數(shù)。

3.4 復(fù)雜模型算例

設(shè)計(jì)均勻半空間中存在兩個(gè)異常體，異常體的模型參數(shù)如圖9所示，源和網(wǎng)格等模型參數(shù)均與3.2.2小節(jié)中半空間模型相同。

圖10為主剖面電磁場(chǎng)響應(yīng)特征，圖中可看出異常體的位置和形狀，低阻異常體相對(duì)介電常數(shù)大，電磁波長(zhǎng)小，電磁場(chǎng)波動(dòng)明顯；高阻體相對(duì)介電常數(shù)略小，異常反應(yīng)不如低阻體靈敏，電磁場(chǎng)波動(dòng)小，說(shuō)明波場(chǎng)對(duì)介質(zhì)不均勻性有一定程度的敏感度，說(shuō)明本文提出的GPR正演算法對(duì)于復(fù)雜模型具有很好的適應(yīng)性。

圖8 主剖面半空間電磁場(chǎng)數(shù)值解與解析解的對(duì)比圖

表1 整個(gè)程序OpenMP并行不同線程計(jì)算效率

圖9 復(fù)雜模型示意圖

4 結(jié)論

本文采用有限單元法，推導(dǎo)了基于行波分解的吸收邊界條件的GPR有限元方程，保留雷達(dá)波傳播的運(yùn)動(dòng)學(xué)特征和動(dòng)力學(xué)特征，實(shí)現(xiàn)了頻率域2.5D探地雷達(dá)正演模擬。主要有以下結(jié)論：

(2) 均勻全空間和半空間模型數(shù)值解與解析解除源處之外的剖面節(jié)點(diǎn)誤差均小于1%，表明算法正確。

圖10 主剖面電磁場(chǎng)分量實(shí)部和虛部

(3) 測(cè)試平臺(tái)CPU-Inter(R) Core(TM) i9-7980XE，主頻2.60 GHz(36核)，內(nèi)存64 GB，64位操作系統(tǒng)，節(jié)點(diǎn)301×301、波數(shù)277時(shí)，并行最優(yōu)線程數(shù)是16，此時(shí)算法耗時(shí)80 s，比一個(gè)線程計(jì)算時(shí)的速度快了3倍，驗(yàn)證了算法的高度并行性和高效性。

(4) 傳統(tǒng)探地雷達(dá)的發(fā)射源在空氣中，本文采用接地線源，這樣做更方便能量的輻射，研究雷達(dá)波在頻率域的傳播特性，完善雷達(dá)的電磁理論，為GPR全波形反演提供重要基礎(chǔ)。