探地雷達頻率域2.5維正演

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

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

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

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

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

1 引言

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

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

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

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

2 方法理論

2.1 控制方程

2.2 有限單元法

3 數值算例

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

3.1 波場參數分析

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

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

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

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

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

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

3.2 正確性驗證

3.2.1 全空間模型

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

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

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

3.2.2 半空間模型

圖6 主剖面電磁場數值解與解析解對比圖

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

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

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

3.3 Open MP并行測試

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

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

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

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

3.4 復雜模型算例

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

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

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

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

圖9 復雜模型示意圖

4 結論

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

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

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

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

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