韓曉冰 郭濤 李航 梁冰洋 周遠國

（西安科技大學(xué)通信與信息工程學(xué)院,西安 710600）

引言

地質(zhì)雷達(ground penetrating radar，GPR)數(shù)值模擬是用電磁計算來研究地下介質(zhì)分布規(guī)律的一種方法，通過將數(shù)值計算獲得的模擬結(jié)果與實際研究區(qū)域的檢測結(jié)果進行對比，能夠有效地分析GPR 技術(shù)的實際應(yīng)用效果，可為實際應(yīng)用中檢測數(shù)據(jù)的地質(zhì)解釋奠定理論基礎(chǔ)，該方法逐漸成為近年來地球物理勘探領(lǐng)域的研究熱點[1-3].常用的數(shù)值模擬方法包括有限差分法[4]、有限元法(finite element method,FEM)[5-6]等.在GPR 建模過程中，為了減少計算量，通常會采用完美匹配層(perfectly matched layer，PML)[7-8]對模型進行截斷，PML 的吸收效果直接影響整個模型的計算精度.

然而，PML 邊界條件在數(shù)值模擬過程中，仍存在一些復(fù)雜問題，例如，對倏逝波都不能很好地吸收[9]，因此在文獻[10]中提出了基于單軸完美匹配層(uniaxial perfectly matched layer，UPML).UPML 不需要對吸收介質(zhì)中的波場進行分裂，因此，計算簡潔、內(nèi)存消耗降低.并且在衰減系數(shù)中引入了吸收因子，使得UPML 邊界能夠較好地衰減倏失波、低頻波等干擾波.文獻[11]利用高階時域有限差分(finitedifference time-domain,FDTD)法結(jié)合UPML 邊界條件實現(xiàn)GPR 正演數(shù)值模擬，優(yōu)化迭代UMPL 吸收邊界，完善了三維模擬在計算效率與吸收效果的雙重要求.文獻[12]采用了基于交替方向隱式FDTD(alternating direction implicit-FDTD,ADI-FDTD)法，引入UMPL 消除了截斷邊界處的強反射現(xiàn)象.但UMPL 媒質(zhì)的角形區(qū)域仍需頻時轉(zhuǎn)換，迭代求解大型線性方程組，增加了內(nèi)存消耗.

雖然UPML 對于倏失波、低頻波等干擾波具有一定的吸收效果，但吸收效果依然不是特別完美.復(fù)頻移PML(complex frequency shifted-PML，CFS-PML)技術(shù)主要利用坐標(biāo)展寬來增強倏失波在場中的衰減，其可以吸收低頻信號在邊界表面的反射，同時對大角度掠射波有很好的吸收效果[13](主要通過頻移參數(shù) α降低掠射波，入射到PML 界面時產(chǎn)生的非均勻波).文獻[14]利用雙線性變換技術(shù)增加CFS 因子，以避免場分裂情況，并且實現(xiàn)了柱坐標(biāo)系下的CFS-PML吸收邊界，與傳統(tǒng)UMPL 邊界相比，CFS-PML 邊界條件具有更好的吸收性能.文獻[15]實現(xiàn)了CFSPML 在三維海洋可控源電磁建模中的應(yīng)用，有效地抑制住了人工邊界效應(yīng)，與Dirichlet 邊界相比，降低了內(nèi)存消耗，提高了計算效率.文獻[16]改進了CFSPML 邊界條件，并成功將其應(yīng)用于瞬變電磁法的三維有限差分虛擬波域(finite difference-fictitious wave domain,FD-FWD)建模中.采用均勻半空間模型，數(shù)值結(jié)果表明在FD-FWD 中，CFS-PML 邊界條件具有比PML 更好的吸收效果.

GPR 數(shù)值模擬時，吸收邊界產(chǎn)生的人為反射波會影響整個模型的計算精度.為進一步改善GPR 三維數(shù)值模擬邊界條件的吸收效果，同時優(yōu)化計算效率，本文引入高階有限元法(higher order finite element method，HO-FEM).HO-FEM 是FEM 的一種改進形式，不同之處在于其采用高階正交多項式作為基函數(shù)[17].一方面提高了算法的效率，另一方面提高了離散網(wǎng)格內(nèi)的建模精度，降低建模時對網(wǎng)格的依賴，且兼顧了常規(guī)FEM 對復(fù)雜地質(zhì)模型模擬的靈活性[18]，該方法可使計算精度較高且計算量大幅減少[19].相較于HO-FEM 利用PML 作為吸收邊界[20]，CFSPML 能更有效地處理倏失波衰減的問題.HO-FEM算法雖然已經(jīng)應(yīng)用于多種電磁場數(shù)值模擬中[21-22]，但在GPR 三維數(shù)值模擬中的研究尚未涉及.因此，本文做了以下兩部分工作：1)將復(fù)頻移坐標(biāo)拉伸系統(tǒng)引入到PML 域內(nèi)，使得頻移參數(shù) α對倏失波進行有效衰減；2)將CFS-PML 算法與HO-FEM 法進行結(jié)合，以此來求解三維GPR 數(shù)值模擬問題，提高計算精度與效率.

1 理論方法

1.1 HO-FEM 方法

從微分形式麥克斯韋方程組出發(fā)，利用時變因子 eiwt，得到如下控制方程：

式中：E為 電場強度，V/m；ω為角頻率；B為磁通量密度，Wb/m2；H為磁場強度，A/m ；D為電通量密度，C/m2；J為電流密度，A/m2；μ 為磁導(dǎo)率；ε為介電常數(shù).對式(1)中第二個方程取旋度，得到關(guān)于磁場矢量赫姆霍茲方程

式中：εr為相對介電常數(shù)；為真空中的波數(shù)；μr為相對介磁常數(shù).與FEM 相同，對整個計算域進行區(qū)域離散，得到很多子區(qū)域，稱之為單元.采用伽遼金加權(quán)殘差法處理式(2)，令e單元內(nèi)的加權(quán)殘差積分為零，即

式中：Φi為加權(quán)基函數(shù)；Ωe是單元計算域.對式(3)展開如下：

式中：Nb表示單元外邊界棱邊數(shù).使用伽遼金加權(quán)殘差法時，規(guī)定采用的加權(quán)基函數(shù)和近似解展開函數(shù)相同，這里采用高斯-洛巴托-勒讓德(Gauss-Lobatto-Legendre,GLL)高階插值多項式作為加權(quán)基函數(shù)，一維參考空間上的N(N≥2)階GLL 基函數(shù)

式中：PN(ξj)為Legendre 多項式；ξj為GLL 插值節(jié)點；j=0,1,···,N.

將所有單元進行集成，離散系統(tǒng)矩陣如下：

由式(1)控制方程推導(dǎo)得到磁場矢量赫姆霍茲方程，結(jié)合式(3)伽遼金加權(quán)殘差法，引入式(5)中的GLL 高階插值多項式，得到式(6) HO-FEM 的離散系統(tǒng)矩陣形式，最終我們將HO-FEM 算法應(yīng)用于整個計算域內(nèi).

1.2 HO-FEM-CFS-PML 計算方法

CFS-PML[23]的主要思想是在人工邊界內(nèi)水平和垂直方向引入一個復(fù)雜的坐標(biāo)拉伸系統(tǒng).例如，x方向拉伸因子e(x)為

在電磁建模中主 要考慮到 參數(shù)k(x)和dx(x)，且k(x)和dx(x)均取正數(shù)，在本文中k(x)取值為25，與 文獻[16]一致.

α(x)參數(shù)的本質(zhì)是吸收倏逝波并改善掠射角處的吸收性能[24]，表達式為

式中：LPML為 每側(cè)PML 區(qū)域厚度；pd取值一般為2 或3[16]；f0為電磁波主頻率.

dx(x)作為x方向的衰減阻尼因子，表達式為

式中：cp為 電磁波傳播速度；R是反射系數(shù).式(10)中的1 為計算域x方向長度縮放為1.從式(10)看出，上述公式中dx(x)須在PML 域內(nèi)取值.對于三維模型而言，PML 域內(nèi)擁有三種媒質(zhì)分區(qū)，分別是平面區(qū)、棱邊區(qū)和角頂區(qū)，每個分區(qū)都應(yīng)有d(x)作為衰減阻尼因子進行電磁波的衰減.

在HO-FEM 中引入CFS-PML 后，赫姆霍茲方程可以表示為

式(12)只針對PML 域內(nèi)赫姆霍茲方程式.為使PML內(nèi)邊界與計算域外邊界形成阻抗匹配，引入CFSPML 方法后三維空間中的相對介電張量以及相對介磁張量為：

結(jié)合HO-FEM 方法，采用伽遼金加權(quán)殘差法處理，將式(6)中的質(zhì)量矩陣和剛度矩陣修正為：

通過以上公式，我們將CFS-PML 復(fù)頻移坐標(biāo)拉伸系統(tǒng)引入到HO-FEM 中的PML 域內(nèi)，對坐標(biāo)空間中的電性參數(shù)進行復(fù)頻移拉伸，從而改善介質(zhì)本構(gòu)張量的矩陣特性，增強了整個PML 域?qū)渴挪ǖ奈招阅?本文并未描述參考空間與物理空間的轉(zhuǎn)換，具體可參考文獻[20].

2 數(shù)值算例

在這一節(jié)中，我們設(shè)計了兩個GPR 模型，第一個模型采用三維分層介質(zhì)結(jié)構(gòu)，第二個模型采用三維復(fù)雜異常體介質(zhì)結(jié)構(gòu).將本文的CFS-PML 方法與傳統(tǒng)PML 方法在數(shù)值精度上做對比，為方便起見，兩種算例都采用線源作為發(fā)射源，極化方向都為y方向，并與商業(yè)軟件WCT 進行比較，來驗證該方法的準(zhǔn)確性.同時，對比特定頻率下反射誤差、相對誤差，來驗證本文算法的有效性和精確性.所有算例的仿真都在Intel(R) Xeon(R) Gold 6248R 的CPU 工作平臺進行.

2.1 分層介質(zhì)GPR 模型

算例1 三維分層介質(zhì)GPR 模型如圖1 所示，計算區(qū)域為4 m×4 m×2 m.以o為坐標(biāo)原點，z=1 m 為分界面，規(guī)定三個方向均為正方向.上層空間為空氣，空氣層厚度1 m.下層空間為土壤，并且在分界面下方0.5 m 處設(shè)置水平層厚0.1 m 的礦體.在分界面上方0.125 m 處，沿x方向兩端設(shè)置50 個接收點，每點相距0.08 m，接收點中心位置坐標(biāo)為(2，2，1.125).發(fā)射頻率設(shè)置為40 MHz，電流強度為25 A，長度為1 m的發(fā)射線源沿y軸放置于分界面上方0.01 m 處，其中心位置坐標(biāo)為(2，2，1.01).表1 給出了算例1 具體的電性參數(shù)[25].

圖1 分層介質(zhì)GPR 模型Fig.1 Layered media GPR model

表1 算例1 電性參數(shù)Tab.1 Electrical parameters of case 1

由圖2 可知，兩種邊界條件的數(shù)值模擬結(jié)果均表明HO-FEM 算法應(yīng)用在GPR 數(shù)值模擬時具有較高精度.其中LMGF 求解器為WCT 中分層介質(zhì)格林函數(shù)，主要應(yīng)用于分層介質(zhì)的電磁仿真.

圖2 算例1 各接收點處磁場強度實部和虛部模擬結(jié)果Fig.2 Simulation results of the real and imaginary parts of the magnetic field strength at each reception point(case 1)

本文算法反射誤差理論計算公式采用文獻[26]方法，如式(18)所示：

式中：Hi為 觀測點對應(yīng)的磁場強度；為解析解中每個觀測點對應(yīng)的磁場強度；為解析解中磁場強度的最大值.

觀察圖3 的反射誤差可知，由于頻移參數(shù) α的影響，CFS-PML 相比于傳統(tǒng)PML 計算的各分量能夠降低10～30 dB.由圖4 可以看到，Hx、Hy、Hz相對誤差分別小于0.72%、0.63%、0.46%，表明CFS-PML 吸收邊界有效地抑制了人工邊界效應(yīng)，使得吸收性能更加優(yōu)異.因此，本文所使用的CFS-PML 吸收邊界在整體的接收點處所得的效果要優(yōu)于PML 吸收邊界.

圖3 算例1 各接收點處磁場強度實部和虛部反射誤差對比Fig.3 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 1)

圖4 算例1 各接收點處磁場強度實部和虛部相對誤差對比Fig.4 Comparison of the relative errors of the real and imaginary parts of the magnetic field strength at each reception point(case 1)

與此同時，在本次算例中，CFS-PML 方法與傳統(tǒng)PML 方法的計算耗時分別為18 044 s 和21 734 s，時間節(jié)省約為17%；內(nèi)存消耗分別約為9.8 GB 和11.0 GB，內(nèi)存節(jié)省約為11%.

2.2 起伏地巖GPR 模型

算例2 三維起伏地巖GPR 模型如圖5 所示，計算區(qū)域為4 m×4 m×4 m.以o為坐標(biāo)原點，z=2.5 m為分界面，規(guī)定三個方向均為正方向.上層空間為空氣，空氣層厚度1.5 m，下層空間土壤厚度為1.1 m，地巖隆起部分上方寬度為0.5 m，隆起高度為0.12 m，地巖未隆起的厚度為1.4 m.在模型結(jié)構(gòu)上采用的是不規(guī)則網(wǎng)格剖分，假設(shè)在巖層下方存在一個棱柱形異常體，用以模擬地質(zhì)空洞.發(fā)射中心頻率設(shè)置為62 MHz，電流強度為25 A，長度為1 m 的發(fā)射線源沿y軸放置在分界面上方0.01 m，發(fā)射線源中心位置坐標(biāo)為(2，2，2.51).接收點中心位置坐標(biāo)為(2，2，2.625)，采用50 個接收點，間隔為0.08 m，沿著x方向展開長度為4 m.表2 給出了算例2 具體的電性參數(shù)[25].采用商業(yè)軟件WCT 中CN-FDTD 算法與本文提出的算法進行對比.

表2 算例2 電性參數(shù)Tab.2 Electrical parameters of case 2

圖5 起伏地巖GPR 模型Fig.5 GPR model of uphill terrane

算例2 雖然受地形改變的影響，但通過對比CNFDTD 算法結(jié)果，HO-FEM 算法仍能保持高精度.在測試過程中可以看到，PML 吸收效果并不能達到100%的吸收效率，但從實驗結(jié)果來看，圖6 觀察到磁場強度變化較快的區(qū)域，CFS-PML 仍能很好地吻合CN-FDTD 結(jié)果，表現(xiàn)出較高的精度.圖7 展示了傳統(tǒng)PML 與CFS-PML 在觀測點的反射誤差，從中可以看到CFS-PML 的反射誤差相比PML 較小.三維起伏地巖模型含有復(fù)雜有耗的結(jié)構(gòu)，并且采用的不規(guī)則網(wǎng)格剖分也會使得傳統(tǒng)的PML 邊界產(chǎn)生一定的誤差，然而CFS-PML 中的頻移參數(shù) α使截斷區(qū)域產(chǎn)生頻率相關(guān)的衰減，提高掠射角處吸收性能，使得CFS-PML 在傳統(tǒng)的PML 基礎(chǔ)上反射誤差能夠降低10 dB 左右，并且相對誤差能夠控制在合理的范圍.從圖8 可以看出CFS-PML 邊界條件下Hx、Hy、Hz相對誤差分別小于0.99 %、0.52 %、0.74 %，也體現(xiàn)了較高的精度.與此同時，算例2 中，CFS-PML 方法與傳統(tǒng)PML 方法的計算耗時分別為54 243 s 和77 211 s，時間節(jié)省約為30%；內(nèi)存消耗分別約為14.6 GB 和16.8 GB，內(nèi)存節(jié)省約為14%.

圖6 算例2 各接收點處磁場強度實部和虛部模擬結(jié)果Fig.6 Simulation results of the real and imaginary parts of the magnetic field strength at each reception point(case 2)

圖7 算例2 各接收點處磁場強度實部和虛部反射誤差對比Fig.7 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 2)

圖8 算例2 各接收點處磁場強度實部和虛部相對誤差對比Fig.8 Comparison of the relative errors of the real and imaginary parts of the magnetic field strength at each reception point(case 2)

3 結(jié)論

本文在頻域HO-FEM 算法的電磁場模擬中成功實現(xiàn)CFS-PML 人工邊界條件，極大地改善了PML 的吸收性能，提高了電磁場GPR 數(shù)值模擬的計算效率.與商業(yè)軟件WCT 數(shù)值模擬結(jié)果對比表明，CFSPML 邊界條件的應(yīng)用具有更高的計算精度和效率，驗證了本文所提出算法的準(zhǔn)確性和可靠性.兩個算例的數(shù)值結(jié)果表明至少能夠降低10 dB 的反射誤差，計算效率上分別提升17%、30%.

值得注意的是，本文公式推導(dǎo)僅考慮各向同性介質(zhì)，并未考慮各向異性介質(zhì)材料.此外，本文基于CFS-PML 的HO-FEM 并未采用加速并行等技術(shù)，以后這將成為我們研究的方向.