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

①(中南大學地球科學與信息物理學院 長沙 410083)

②(中南大學有色金屬成礦預測與地質環(huán)境監(jiān)測教育部重點實驗室 長沙 410083)

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

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

1 引言

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

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

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

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

2 方法理論

2.1 控制方程

2.2 有限單元法

3 數(shù)值算例

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

3.1 波場參數(shù)分析

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

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

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

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

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

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

3.2 正確性驗證

3.2.1 全空間模型

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

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

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

3.2.2 半空間模型

圖6 主剖面電磁場數(shù)值解與解析解對比圖

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

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

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

3.3 Open MP并行測試

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

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

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

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

3.4 復雜模型算例

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

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

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

表1 整個程序OpenMP并行不同線程計算效率

圖9 復雜模型示意圖

4 結論

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

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

圖10 主剖面電磁場分量實部和虛部

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

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