基于CFS-PML 的高階有限元在地質(zhì)雷達(dá)數(shù)值模擬中的應(yīng)用研究

韓曉冰 郭濤 李航 梁冰洋 周遠(yuǎn)國(guó)

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

引言

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

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

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

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

1 理論方法

1.1 HO-FEM 方法

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

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

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

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

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

式中：PN(ξj)為L(zhǎng)egendre 多項(xiàng)式；ξj為GLL 插值節(jié)點(diǎn)；j=0,1,···,N.

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

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

1.2 HO-FEM-CFS-PML 計(jì)算方法

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

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

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

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

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

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

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

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

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

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

2 數(shù)值算例

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

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

算例1 三維分層介質(zhì)GPR 模型如圖1 所示，計(jì)算區(qū)域?yàn)? m×4 m×2 m.以o為坐標(biāo)原點(diǎn)，z=1 m 為分界面，規(guī)定三個(gè)方向均為正方向.上層空間為空氣，空氣層厚度1 m.下層空間為土壤，并且在分界面下方0.5 m 處設(shè)置水平層厚0.1 m 的礦體.在分界面上方0.125 m 處，沿x方向兩端設(shè)置50 個(gè)接收點(diǎn)，每點(diǎn)相距0.08 m，接收點(diǎn)中心位置坐標(biāo)為(2，2，1.125).發(fā)射頻率設(shè)置為40 MHz，電流強(qiáng)度為25 A，長(zhǎng)度為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ù)值模擬時(shí)具有較高精度.其中LMGF 求解器為WCT 中分層介質(zhì)格林函數(shù)，主要應(yīng)用于分層介質(zhì)的電磁仿真.

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

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

式中：Hi為 觀測(cè)點(diǎn)對(duì)應(yīng)的磁場(chǎng)強(qiáng)度；為解析解中每個(gè)觀測(cè)點(diǎn)對(duì)應(yīng)的磁場(chǎng)強(qiáng)度；為解析解中磁場(chǎng)強(qiáng)度的最大值.

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

圖3 算例1 各接收點(diǎn)處磁場(chǎng)強(qiáng)度實(shí)部和虛部反射誤差對(duì)比Fig.3 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 1)

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

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

2.2 起伏地巖GPR 模型

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

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

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

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

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

圖7 算例2 各接收點(diǎn)處磁場(chǎng)強(qiáng)度實(shí)部和虛部反射誤差對(duì)比Fig.7 Comparison of real and imaginary reflection errors of magnetic field intensity at each reception point(case 2)

圖8 算例2 各接收點(diǎn)處磁場(chǎng)強(qiáng)度實(shí)部和虛部相對(duì)誤差對(duì)比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 算法的電磁場(chǎng)模擬中成功實(shí)現(xiàn)CFS-PML 人工邊界條件，極大地改善了PML 的吸收性能，提高了電磁場(chǎng)GPR 數(shù)值模擬的計(jì)算效率.與商業(yè)軟件WCT 數(shù)值模擬結(jié)果對(duì)比表明，CFSPML 邊界條件的應(yīng)用具有更高的計(jì)算精度和效率，驗(yàn)證了本文所提出算法的準(zhǔn)確性和可靠性.兩個(gè)算例的數(shù)值結(jié)果表明至少能夠降低10 dB 的反射誤差，計(jì)算效率上分別提升17%、30%.

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