馮德山，陳承申，王洪華

1 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083

2 有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室，長沙 410083

3 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410083

1 引 言

探地雷達(dá)（Ground penetrating radar，GPR）正演模擬一直是該領(lǐng)域理論研究的熱點(diǎn)[1-2］，通過對(duì)典型地質(zhì)模型正演結(jié)果的分析，可以加深對(duì)GPR反射剖面的理解與認(rèn)識(shí)[3-4］，對(duì)GPR資料解釋具有指導(dǎo)作用.在GPR正演模擬中，時(shí)域有限差分法[5-7］（FDTD）因簡單靈活而被廣泛采用，但它不適合復(fù)雜的物性分界面.有限單元法（Finite Element Method，F(xiàn)EM）不需要計(jì)算內(nèi)部邊界條件，適用于物性復(fù)雜問題、求解過程規(guī)范化，具有廣泛的適應(yīng)性.

在地震波FEM模擬[8］得到快速發(fā)展的同時(shí)，一些學(xué)者將其引入到GPR正演模擬中：沈飚[9］進(jìn)行了GPR簡單模型正演；底青云和王妙月[10-11］推導(dǎo)了含衰減項(xiàng)的探地雷達(dá)波有限元方程，實(shí)現(xiàn)了復(fù)雜形體的GPR波的FEM正演模擬和偏移處理；Fanning和Boothby[12］將GPR有限元建模技術(shù)應(yīng)用到拱橋檢測中；謝輝等[13］基于二十節(jié)點(diǎn)等參單元的FEM對(duì)Pulse型GPR進(jìn)行了模擬；Arias等[14］將有限元技術(shù)應(yīng)用到GPR考古成像中；Lu等[15］利用離散時(shí)域Galerkin進(jìn)行色散介質(zhì)中的GPR正演.

GPR正演模擬是用有限空間模擬半無限大空間，在模型邊界上必然產(chǎn)生邊界反射.因此，需要在截?cái)噙吔缣幨褂靡环N特殊計(jì)算方法，在保證邊界處場的計(jì)算精度的同時(shí)大大消除人工截?cái)噙吔绲姆瓷洳?如Cerjan邊界條件[16］，它是在邊界區(qū)引入一個(gè)衰減層，讓波的能量在這個(gè)衰減層里呈現(xiàn)指數(shù)衰減，但是這種吸收邊界條件的選取比較困難，衰減層的厚度、衰減系數(shù)的大小都要憑經(jīng)驗(yàn)去估計(jì).Sarma[17］在邊界區(qū)域建立一個(gè)一定厚度的衰減層，根據(jù)材料性質(zhì)估算恒定比例系數(shù)，使衰減層和介質(zhì)區(qū)存在明顯的物性差異，從而產(chǎn)生反射波；王月英等[18］提出了一種改進(jìn)的Sarma邊界條件，通過在衰減層內(nèi)增加一定厚度的過渡帶，過渡帶內(nèi)的比例系數(shù)呈現(xiàn)線性變化，這就使介質(zhì)區(qū)和衰減層之間的物性差異漸變，大大減弱了交界面處的反射，但也存在對(duì)阻尼比的選取不敏感、過渡帶和吸收層的交界面的角點(diǎn)處產(chǎn)生較強(qiáng)的繞射波的缺點(diǎn)，仍需尋求新的解決辦法來取得更好的吸收效果.

本文通過對(duì)透射邊界條件和改進(jìn)Sarma邊界條件的研究，提出一種混合邊界條件.通過在研究區(qū)域外加載一定厚度的吸收層，以改進(jìn)的Sarma邊界條件對(duì)到達(dá)邊界處的GPR波能量進(jìn)行吸收衰減，然后在最外邊加載透射邊界條件，將衰減過的探地雷達(dá)波透射出去，這種混合邊界條件集成了二者的優(yōu)勢，對(duì)人工截?cái)噙吔绲姆瓷洳ㄈ〉酶玫奈招Ч?

2 GPR波動(dòng)方程推導(dǎo)

Maxwell方程組描述了電磁場的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)規(guī)律[19］，由電磁波理論，高頻電磁波在媒質(zhì)中的傳播規(guī)律也應(yīng)服從Maxwell方程組，其方程組的微分形式可以表示為[20-21］

式中，E為電場強(qiáng)度（V／m），H為磁場強(qiáng)度（A／m），B為磁感應(yīng)強(qiáng)度（T），D為電位移（C／m2），J為電流密度（A／m2），ρ為電荷密度（C／m3）.

本構(gòu)方程：

式中ε為介電常數(shù)（F／m），μ為磁導(dǎo)率（H／m），σ為電導(dǎo)率（S／m）.

把方程式（2）代入方程式（1a），并求旋度，得

考慮到式

由于ε，μ，σ為坐標(biāo)的函數(shù)，它們的時(shí)間導(dǎo)數(shù)可以忽略，用電場E取代（4）式中的A，則（4）式中右邊第二項(xiàng)為0，整理得

假設(shè)將雷達(dá)波激勵(lì)源如電場源SE或磁場源SH，代入到（5）式中，有

同理對(duì)（1b）式中第二式兩邊求旋度，可得下式：

（6）、（7）兩式表明，磁場H和電場E及其分量滿足相同的微分方程.

3 GPR波動(dòng)方程有限元求解

求解GPR波動(dòng)方程其實(shí)質(zhì)是結(jié)合雷達(dá)波所要滿足的初、邊值條件求解偏微分方程.邊界條件分為兩種類型[22］：Dirichlet邊界條件與Neumann邊界條件.變分法推導(dǎo)結(jié)果自動(dòng)滿足Dirichlet邊界條件，即未知函數(shù)在邊界上有已知值；伽遼金（Galerkin）法推導(dǎo)結(jié)果自動(dòng)滿足Neumann邊界條件，即未知函數(shù)的導(dǎo)數(shù)在邊界上有已知值.下面采用伽遼金法推導(dǎo)含衰減項(xiàng)的GPR有限元方程.

第1步 網(wǎng)格剖分

將求解的二維區(qū)域剖分成矩形單元，如圖1所示.

圖1 網(wǎng)格剖分及節(jié)點(diǎn)編號(hào)示意圖Fig.1 Sketch map of mesh division and node number

第2步 線性插值

設(shè)（x，y）是子單元的坐標(biāo)，（x0，y0）是子單元中點(diǎn)的坐標(biāo)，a、b是子單元的兩個(gè)邊長，（ξ，η）為母單元的坐標(biāo)，兩個(gè)單元間的坐標(biāo)變換關(guān)系為

x=x0+（a／2）·ξ， y=y0+（b／2）·η， （8）其中微分關(guān)系為

單元形函數(shù)為雙線性插值形函數(shù)，可寫為

第3步 單元積分

根據(jù)伽遼金法，將（6）式兩邊同時(shí)乘以δE，并求積分，有

則式（12）左邊第一項(xiàng)為

其中Me為單元質(zhì)量矩陣：

式（12）左邊第二項(xiàng)為

其中K′e為單元阻尼矩陣：

式（12）左邊第三項(xiàng)為

其中Ke= （kij）= （kji）.

將（11）式求得ξ或者η的微商，代入到（18）式積分即可得到單元?jiǎng)偠染仃噆ij為

式（12）右邊項(xiàng)為

第4步 總體合成

根據(jù)式（13）、（15）、（17）、（20）得到單元積分：

其中ND是節(jié)點(diǎn)總數(shù).將4×4的系數(shù)矩陣Me、K′e和Ke擴(kuò)展成ND×ND的矩陣M，K′和K，將列向量Se擴(kuò)展成ND維列向量S，即

其中，

由于δET≠0，所以必有M¨E+K′E·+KE=SE.于是可以得到時(shí)空域GPR波的二維有限元方程，這是含有ND個(gè)元的ND個(gè)方程聯(lián)合的常微分方程組，可以寫成如下形式：

第5步 解GPR有限元常微分方程

解常微分方程組（24）或（25）前，要代入邊界條件，邊界條件的加載方法將在下節(jié)討論.先采用中心差分法解方程（24）.對(duì)初始條件離散化得到E（0）=把時(shí)間域 [0，T］分為幾個(gè)相等步長Δt=T／n，把時(shí)刻t的微分方程記為

用差商代替微分：

中心差分法是一種條件穩(wěn)定的計(jì)算方法，當(dāng)時(shí)間步長Δt取的過大時(shí)，計(jì)算結(jié)果就會(huì)出現(xiàn)數(shù)值色散，為此本文采用 Δt≤ （ΔXmin／Vmax），式中時(shí)間步長 Δt與最小單元邊長ΔXmin成正比，與最大媒質(zhì)速度Vmax成反比.差分網(wǎng)格空間步長滿足穩(wěn)定性條件Δx＜ （λmin／10），其中λmin為最小波長長度，Δx為單元的最小尺寸.將式（27）代入到式（26）中得到

求解不含衰減項(xiàng)的GPR波動(dòng)方程時(shí)，由于系數(shù)矩陣A往往是大型的病態(tài)稀疏矩陣，其條件數(shù)很大，如果采取對(duì)其直接求逆，計(jì)算量巨大，其求解結(jié)果不可信，只有當(dāng)系數(shù)矩陣A是對(duì)角陣時(shí)，其逆矩陣是對(duì)角元素求倒數(shù)，這種方法才適用，為此采取集中質(zhì)量矩陣方法來處理M矩陣，即將每一行（或列）的元素都加到對(duì)角線元素上去，則式（14）可轉(zhuǎn)化為單元集中質(zhì)量矩陣：

由于總體集中質(zhì)量矩陣也有對(duì)角陣，所以不含衰減項(xiàng)的GPR有限元方程的差分迭代公式可寫為

在求解含有衰減項(xiàng)的波動(dòng)方程（28）時(shí)，系數(shù)矩陣A是質(zhì)量矩陣M 和阻尼矩陣K′的線性組合，同樣可以采用單元集中質(zhì)量矩陣求解.但本文中選用不完全LU分解預(yù)處理的BICGSTAB（Biconjugate Gradient Stabilized）線性方程組求解算法[23-24］，該算法是一種高效迭代法，特別適用于求解條件數(shù)很大的病態(tài)線性方程組.

4 混合邊界條件理論

目前應(yīng)用于電磁場正演模擬的邊界條件主要集中于FDTD相關(guān)算法，如輻射邊界條件[25-26］、基于單行波方程的吸收邊界條件[27-28］、超吸收邊界條件[29］、完全匹配層（PML）[30-31］吸收邊界條件等.而GPR有限元法正演模擬文獻(xiàn)較少，采用的吸收邊界條件多借鑒與地震波有限元正演模擬的吸收邊界條件.下面主要介紹透射邊界條件與Sarma吸收邊界條件，在此基礎(chǔ)上提出了混合邊界條件.

4.1 透射邊界條件

透射邊界條件是通過直接模擬各種單向波動(dòng)的共同運(yùn)動(dòng)學(xué)特征來建立的一種邊界條件，原理簡單易懂，算法容易實(shí)現(xiàn).首先給出穿過人工邊界上一點(diǎn)并沿邊界外法線方向傳播的單向波動(dòng)的一般表達(dá)式，并假定所有單向波動(dòng)以同一人工波速沿法線方向從邊界透射出去，由此可得透射邊界條件公式[32］：

式中L為微分算子.在每個(gè)單元中有

這里Ei為每個(gè)單元上第i個(gè)節(jié)點(diǎn)的電場值.殘余量R為

利用伽遼金余量法，二維有限元方程可寫為

其中Ae為單元面積，r是殘余量R的向量表示，由式（32）代入式（33）得到.

由于內(nèi)邊界的積分相互抵消，所以只需求解圖2所示的區(qū)域外邊界積分.在均勻各向同性介質(zhì)中，

圖2 區(qū)域外邊界條件Fig.2 Boundary condition outside the area

根據(jù)Claerbout[33］推導(dǎo)的旁軸近似，其下行波方程為

左行波方程為

右行波方程為

由法向?qū)?shù)的定義知

其中，nx和ny為邊界外法線的方向余弦.

將式（31）、（32）、（33）代入式（34）中，并利用高斯定理可得

其中Γl、Γr、Γd分別表示左、右及底邊界，上邊界為自由邊界條件.

將式（35）、（36）、（37）給出的邊界條件代入到式（39），可以得到：

在加載透射邊界條件時(shí)，還需求解外法線的方向余弦，其求解示意如圖3所示.θx為GPR發(fā)射天線激勵(lì)源和邊界線單元x方向的夾角，θy為y方向夾角，無論這兩個(gè)角是正還是負(fù)，其余弦值nx和ny均為正值.

圖3 外法線的方向余弦角示意圖Fig.3 Sketch map of direction Cosine angle of external normal

在求解左邊界的邊界條件時(shí)，如圖4所示，令N1= （1-ξ）／2，N2= （1+ξ）／2，dΓ=l·dξ／2，這里l為邊界上的單元寬度，例如AB、CD邊l=b，BC邊l=a.先分析邊界單元在AB邊的情況：

圖4 邊界上的單元Fig.4 Elements on the boundary

將邊界的單元系數(shù)矩陣定義為Fe，則左邊界的邊界單元系數(shù)矩陣擴(kuò)展為Fle：

同理，可以求得底邊界單元系數(shù)矩陣Fde和右邊界單元系數(shù)矩陣Fre：

故邊界的阻尼矩陣F為

將邊界阻尼加入到GPR有限元方程中去，得到如下式子：

將（45）式透射邊界條件結(jié)合式（24）就可以進(jìn)行GPR有限元正演模擬.

4.2 Sarma吸收邊界條件

Sarma吸收邊界條件是在求解空間外的區(qū)域增加一個(gè)由衰減性媒質(zhì)組成的衰減層，如圖5a所示.對(duì)于衰減性媒質(zhì)，阻尼力作用比較復(fù)雜，導(dǎo)致阻尼矩陣C不容易確定.Rayleigh在1877年給出了一種計(jì)算阻尼矩陣C的經(jīng)典方法，該方法[34］是利用整體質(zhì)量矩陣M和整體剛度矩陣K來計(jì)算阻尼矩陣C：

圖5 Sarma邊界條件及其改進(jìn)模型示意圖（a）Sarma邊界條件；（b）改進(jìn)的Sarma邊界條件；（c）優(yōu)化過渡帶角點(diǎn)的Sarma邊界條件.Fig.5 Sketch map of Sarma boundary condition and improved model（a）Sarma boundary condition；（b）Improved Sarma boundary condition；（c）Sarma boundary condition of optimized transitional layer.

其中，

式中ωi和ωj分別為介質(zhì)的第i和第j個(gè)固有的圓頻率，ξi和ξj為對(duì)應(yīng)的阻尼比.

Caughey[35］提出了一種更為普遍的阻尼矩陣計(jì)算公式：

可以看出（46）式是（47）式的二階形式，即m=2，體現(xiàn)了質(zhì)量矩陣和剛度矩陣的線性組合，并且這種組合非常容易計(jì)算，其單元中的計(jì)算公式如下：

對(duì)于二階Caughey阻尼矩陣的兩個(gè)系數(shù)a0和a1與阻尼比ξ存在著如下的關(guān)系：

其中ω為發(fā)射天線的中心頻率，通常0.05≤ξ≤0.30.將二階Caughey阻尼矩陣引入到GPR吸收邊界的處理中，得到吸收邊界區(qū)內(nèi)的有限元方程為

在衰減層內(nèi)的比例系數(shù)a0和a1可由式（49）得到，這種方法雖然可行，但是衰減層和介質(zhì)區(qū)還是存在著明顯的物性差異，這種差異必然會(huì)在交界面處產(chǎn)生人為的反射波.針對(duì)Sarma邊界條件在交界面處產(chǎn)生的人為反射波，王月英等[18］提出了一種改進(jìn)方法，如圖5b所示：即在衰減層內(nèi)添加一個(gè)過渡帶，過渡帶內(nèi)比例系數(shù)由零逐漸增大到a0和a1，衰減層取一個(gè)半波長的厚度，過渡帶大約占半個(gè)波長.但該改進(jìn)的Sarma邊界條件在過渡帶的四個(gè)角點(diǎn)位置仍存在著差異明顯的四個(gè)小邊界，為此本文采用一圈一圈的環(huán)形過渡帶加載辦法，消除了這種角點(diǎn)的人為差異，稱之為優(yōu)化過渡帶角點(diǎn)的Sarma邊界條件，如圖5c所示.

其中，衰減層內(nèi)總節(jié)點(diǎn)數(shù)為n，過渡帶內(nèi)的節(jié)點(diǎn)數(shù)為m，m＜n（0≤i＜m）為過渡帶內(nèi)節(jié)點(diǎn)i處的比例系數(shù)，比例系數(shù)從交界面向衰減層邊緣方向呈線性增長關(guān)系，這種方法就減緩了交界面處的物性參數(shù)突變，從而減小交界面處產(chǎn)生的人為反射波能量，比例系數(shù)增大到a0和a1就保持恒定，以確保GPR波在衰減層內(nèi)能得到充分吸收.既能使交界面處的反射得到大大降低，又能讓邊界處的GPR波能量被充分吸收.

4.3 混合吸收邊界條件

通過對(duì)透射邊界條件和Sarma邊界條件的深入分析，提出了一種結(jié)合兩種邊界條件的混合邊界條件（如圖6所示），其主要思想是利用Sarma邊界條件對(duì)到達(dá)邊界區(qū)域的電磁波能量衰減功能和透射邊界對(duì)電磁波能量的透射功能，使GPR波經(jīng)過Sarma邊界條件的衰減吸收后再通過透射邊界條件將剩余能量透射出去，集成了二者的優(yōu)勢.處理辦法是在計(jì)算區(qū)域外面加一吸收層，用改進(jìn)的Sarma邊界條件使GPR波能量被吸收減弱，然后在最外邊加載透射邊界條件.通過改變公式（48）中給的阻尼比ξ，可以實(shí)現(xiàn)調(diào)節(jié)Sarma邊界的吸收作用，從而和透射邊界條件相互協(xié)作，達(dá)到最佳吸收效果.

5 不同邊界條件的有限元模擬對(duì)比實(shí)例

圖6 混合邊界條件模型示意圖Fig.6 Sketch map of mixed boundary condition

運(yùn)用Matlab編制的有限元GPR模擬程序?qū)?0.0m×5.0m均勻模型進(jìn)行了模擬，介質(zhì)的介電常數(shù)ε=3.0，電導(dǎo)率σ=0.001S·m-1.網(wǎng)格單元大小為0.1m的正方形，GPR脈沖激勵(lì)源的子波形式為f（t）=t2e-atsinω0t，其中ω0為發(fā)射天線中心頻率為100MHz，電磁波脈沖的衰減速度取決于系數(shù)α，此處取α=0.93ω0.得到了半空間圖7所示的不帶邊界條件的雷達(dá)正演wiggle圖.圖中可見，A、B為左右邊界的反射波，C為底邊界的反射波，D、E為左下角和右下角的兩個(gè)角點(diǎn)的繞射波.這些人為截?cái)噙吔绲姆瓷鋸?qiáng)度很大，對(duì)觀測區(qū)域影響嚴(yán)重，說明了對(duì)截?cái)噙吔绲奶幚韯菰诒匦?

5.1 透射邊界條件中心脈沖激勵(lì)源模型仿真

圖7 不帶邊界條件的均勻模型wiggle圖Fig.7 Wiggle diagram for homogeneous model without boundary conditions

圖8 均勻模型示意圖Fig.8 Sketch map of homogeneous model

為了更好地說明不同吸收邊界條件對(duì)人工截?cái)噙吔绶瓷洳ǖ奈招Ч?，建立圖8所示9.0m×9.0m大小的正方形均勻模型，并把脈沖激勵(lì)源置于模擬區(qū)域的正中心，采樣時(shí)間間隔為0.5ns，其它參數(shù)與無邊界條件情況下一致.通過截取不同時(shí)刻的波場快照?qǐng)D，觀測不同邊界條件對(duì)人工截?cái)噙吔绲奈招Ч?圖9（a，b，c）分別為不帶邊界條件的20、30ns、35ns波場快照，圖中可見在30ns時(shí)刻GPR波前開始傳到區(qū)域邊緣，人為截?cái)噙吔绶瓷洳ㄩ_始形成，從35ns快照?qǐng)D可以看出，人為截?cái)噙吔绶瓷浠氐腉PR波能量很強(qiáng)，嚴(yán)重影響到對(duì)目標(biāo)區(qū)域的研究，同樣說明，截?cái)噙吔绲奶幚硎欠浅１匾?

圖10（a，b，c）分別為模型邊界處加載透射邊界的30、35、40ns波場快照.圖中可見，在30ns時(shí)，波的能量剛剛傳播到邊界區(qū)域，在35ns時(shí)可以看出電磁波大部分的能量都從邊界處透射出去，只有少部分的能量反射回來，對(duì)比圖9c相同時(shí)刻的截?cái)噙吔缣幍膹?qiáng)反射波能量，說明了透射邊界條件取得了顯著效果.但是再分析圖10c中40ns的波場快照?qǐng)D仍然可以看到截?cái)噙吔绲姆瓷浠夭?，說明透射邊界條件仍有待改進(jìn).

5.2 改進(jìn)的Sarma邊界條件中心脈沖激勵(lì)源模型仿真

仍以上例模型和相應(yīng)參數(shù)為例，分別以Sarma邊界條件、改進(jìn)的Sarma邊界條件和優(yōu)化了過渡帶的Sarma邊界條件對(duì)截?cái)噙吔邕M(jìn)行了處理.圖11a為沒有加載過渡帶的Sarma邊界條件的30ns波場快照?qǐng)D，由于沒有加載過渡帶，在介質(zhì)區(qū)和衰減層之間的界面產(chǎn)生了明顯的反射和角點(diǎn)繞射；圖11b為改進(jìn)Sarma邊界條件的30ns波場快照?qǐng)D，由于加載了過渡帶，在介質(zhì)區(qū)和邊界吸收層之間的媒質(zhì)特性是漸變的，明顯的減弱了介質(zhì)之間的差異，降低了介質(zhì)區(qū)和邊界吸收層之間界面的反射波能量，較好地處理了介質(zhì)區(qū)和衰減層之間物性差異引起的反射波，但在過渡帶的四個(gè)角點(diǎn)位置存在著差異明顯的四個(gè)小邊界.為了消除這種因過渡帶人為加載方式引起的反射，本文采用一圈一圈的環(huán)形過渡帶加載辦法，圖11c為優(yōu)化過渡帶Sarma邊界條件的30ns波場快照?qǐng)D.它是在改進(jìn)的Sarma邊界條件基礎(chǔ)上，對(duì)過渡帶的加載方法進(jìn)行了優(yōu)化，目的是減少人為的角點(diǎn)引起的異常，雖然效果不是很明顯，但是從側(cè)面反映出改進(jìn)的Sarma邊界條件的過渡帶的角點(diǎn)異常主要是由角點(diǎn)引起，過渡帶的加載方法并非主要影響因素，重點(diǎn)工作應(yīng)該放在研究過渡帶的角點(diǎn)異常消除上.

5.3 混合吸收邊界條件中心中心脈沖激勵(lì)源模型仿真

如圖12a為加載了透射邊界條件的35ns快照?qǐng)D，圖12b為加載了優(yōu)化過渡帶的Sarma邊界條件的35ns快照?qǐng)D，圖12c為加載混合邊界條件的35ns快照?qǐng)D.圖中可見，圖12b優(yōu)化過渡帶的Sarma邊界條件由于過渡帶的四個(gè)角點(diǎn)區(qū)的繞射波能量比較強(qiáng)，其對(duì)GPR波的吸收效果并不是很理想，而圖12c中的混合邊界條件由于集成了Sarma邊界條件和透射邊界條件的優(yōu)點(diǎn)，先衰減了到達(dá)邊界區(qū)域的GPR波能量，然后再將其透射出去，觀察相同時(shí)刻快照?qǐng)D，可以看出，與圖9c無邊界條件的波場快照?qǐng)D對(duì)比，各種邊界條件對(duì)截?cái)噙吔缣幍姆瓷洳ǘ既〉昧肆己眯Ч?，但混合邊界條件明顯優(yōu)于單純的透射邊界條件或Sarma邊界條件.

6 探地雷達(dá)正演模擬實(shí)例

6.1 起伏界面模型

如圖13所示，模型為一個(gè)10.0m×5.0m的矩形區(qū)域，分為上下兩層介質(zhì)，中間是起伏界面，上層介質(zhì)的相對(duì)介電常數(shù)ε1為3.0，電導(dǎo)率σ1為0.001S·m-1，下層介質(zhì)的相對(duì)介電常數(shù)ε2為20.0，電導(dǎo)率σ2為0.01S·m-1，整個(gè)區(qū)域由單元邊長0.1m的正方形，網(wǎng)格剖分為100×50的網(wǎng)格空間，GPR波脈沖激勵(lì)源的中心頻率為100MHz，采樣時(shí)窗長度為100ns，采樣時(shí)間間隔為0.5ns.應(yīng)用本文的基于混合邊界條件的FEM算法對(duì)起伏界面模型進(jìn)行了正演模擬，得到了圖14所示的GPR正演模擬wiggle圖，圖中可見，在30ns附近有一條能量較強(qiáng)的反射界面，通過計(jì)算可以得出它與圖13模型中的上、下兩層分界面位置相吻合，但是在起伏比較大的地方，仍存在角點(diǎn)繞射現(xiàn)象.

6.2 三角形與階梯模型

圖13 起伏界面模型示意圖Fig.13 Sketch map of undulating interface

圖14 起伏界面模型探地雷達(dá)正演模擬wiggle剖面圖Fig.14 Wiggle section for forward simulation of GPR

如圖15所示，模型為一個(gè)10.0m×5.0m的矩形區(qū)域，分為上下兩層介質(zhì)，中間是特殊形態(tài)的起伏界面，包括一個(gè)三角形和一個(gè)階梯，上層介質(zhì)的相對(duì)介電常數(shù)ε1為3.0，電導(dǎo)率σ1為0.001S·m-1，下層介質(zhì)的相對(duì)介電常數(shù)ε2為45.0，電導(dǎo)率σ2為0.1S·m-1，整個(gè)區(qū)域由單元邊長0.1m的正方形網(wǎng)格剖分為100×50的網(wǎng)格空間，GPR波源、時(shí)窗長度、采樣時(shí)間間隔均與上面模型一致.應(yīng)用FEM算法對(duì)起伏界面模型進(jìn)行了正演模擬，得到了圖16所示的三角形與階梯模型的GPR正演模擬wiggle圖，分析該剖面圖可以看出，在三角形頂點(diǎn)、階梯形的斷點(diǎn)處都存在較強(qiáng)的角點(diǎn)繞射波，而相應(yīng)的平界面還會(huì)產(chǎn)生較強(qiáng)的反射波，模擬結(jié)果說明了FEM進(jìn)行GPR正演模擬的正確性與可性性.

7 結(jié)論及建議

（1）推導(dǎo)了GPR有限元波動(dòng)方程，給出了有限單元法求解該泛函變分問題的詳細(xì)解法，并結(jié)合差分離散時(shí)，時(shí)間步長與空間步長須滿足的數(shù)值穩(wěn)定性條件，探討了GPR有限元波動(dòng)方程是否帶衰減項(xiàng)的不同求解算法.

（2）闡述了透射邊界條件和Sarma邊界條件的原理.著重對(duì)Sarma邊界條件進(jìn)行了研究，并通過在Sarma邊界條件衰減層內(nèi)以環(huán)狀方式加入過渡層，壓制了介質(zhì)區(qū)和衰減層交界面處的人為反射，優(yōu)化了Sarma邊界條件.并以二維均勻模型的中心脈沖激勵(lì)源為例，對(duì)比了人工截?cái)噙吔缣幉煌吔鐥l件的處理效果.

（3）提出了一種結(jié)合透射邊界條件和Sarma邊界條件的混合邊界條件，使GPR波經(jīng)過Sarma邊界條件的衰減吸收后再通過透射邊界條件將剩余能量透射出去，集成了兩種邊界條件的優(yōu)勢.最后應(yīng)用基于該混合邊界條件的FEM算法對(duì)兩個(gè)典型地電模型進(jìn)行了正演模擬，模擬結(jié)果表明，混合邊界條件有效地消除了截?cái)噙吔缣幍膹?qiáng)反射，說明了FEM算法進(jìn)行GPR正演模擬的有效性及可行性.

致 謝 中國科學(xué)院地質(zhì)與地球物理研究所底青云教授對(duì)論文的有限元理論、中國石油大學(xué)（北京）資源與信息學(xué)院地質(zhì)與地球物理綜合研究中心的王月英老師對(duì)Sarma邊界條件程序編制方面給予了有益幫助，在此一并致以衷心的謝意.

（References）

[1]Feng D S，Dai Q W.GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.NDT &E International，2011，44（6）：495-504.

[2]Shaari A，Ahmad R S，Chew T H.Effects of antenna-target polarization and target-medium dielectric contrast on GPR signal from non-metal pipes using FDTD simulation.NDT &E International，2010，43（5）：403-408.

[3]James I，Rosemary K.Numerical modeling of groundpenetrating radar in 2-D using MATLAB.Computers &Geosciences，2006，32（9）：1247-1258.

[4]Giannopoulos A.Modelling ground penetrating radar by GprMax.Construction and Building Materials，2005，19（10）：755-762.

[5]劉四新，曾昭發(fā)，徐波.三維頻散介質(zhì)中地質(zhì)雷達(dá)信號(hào)的FDTD數(shù)值模擬.吉林大學(xué)學(xué)報(bào) （地球科學(xué)版），2006，36（1）：123-127.Liu S X，Zeng Z F，Xu B.FDTD simulation for ground penetrating radar signal in 3-Dimensional dispersive medium.Journal of Jilin University （Earth Science Edition）（in Chinese），2006，36（1）：123-127.

[6]Bergmann T，Robertsson J O A，Holliger K.Numerical properties of staggered finite-difference solutions of Maxwell＇s equations for ground-penetrating radar modeling.Geophysical Research Letters，1996，23（1）：45-48.

[7]Carcione J M.Radiation patterns for 2-D GPR forward modeling.Geophysics，1998，63（2）：424-430.

[8]王妙月，郭亞曦，底青云.二維線性流變體波的有限元模擬.地球物理學(xué)報(bào)，1995，36（4）：499-506.Wang M Y，Guo Y X，Di Q Y.2-D finite element modelling for seismic wave in media with linear rheological property.Acta Geophysica Sinica （in Chinese），1995，36（4）：499-506.

[9]沈飚.探地雷達(dá)波波動(dòng)方程研究及其正演模擬.物探化探計(jì)算技術(shù)，1994，16（1）：29-33.Shen B.A study on wave equation theory for groundpenetrating radar and forward modelling.Computing Techniques for Geophysical and Geochemical Exploration（in Chinese），1994，16（1）：29-33.

[10]底青云，王妙月.雷達(dá)波有限元仿真模擬.地球物理學(xué)報(bào)，1999，42（6）：818-825.Di Q Y，Wang M Y.2Dfinite element modeling for radar wave.Chinese J.Geophys.（in Chinese），1999，42（6）：818-825.

[11]Di Q Y，Wang M Y.Migration of ground-penetrating radar data with a finite-element method that considers attenuation and dispersion.Geophysics，2004，69（2）：472-477.

[12]Fanning P J，Boothby T E.Three-dimensional modelling and full-scale testing of stone arch bridges.Computer &Structures，2001，79（29）：2645-2662.

[13]謝輝，鐘燕輝，蔡迎春.電磁場有限元法在GPR正演模擬中的應(yīng)用.河南科學(xué)，2003，21（3）：295-298.Xie H，Zhong Y H，Cai Y C.Application of finite element methods for electromagnetic fields in forward model of GPR.Henan Science（in Chinese），2003，21（3）：295-298.

[14]Arias P， Armesto J， Capua D D， et al.Digital photogrammetry，GPR and computational analysis of structural damages in a mediaeval bridge.Engineering Failure Analysis，2007，14（8）：1444-1457.

[15]Lu T，Cai W，Zhang P W.Discontinuous galerkin timedomain method for GPR simulation in dispersive media.IEEE Transactions on Geoscience and Remote Sensing，2005，43（1）：72-80

[16]Cerjan C，Kosloff D，Kosloff R.A nonreflecting boundary condition for discrete acoustic and elastic wave equations.Geophysics，1985，50（4）：705-708.

[17]Sarma G S，Mallick K，Gadhinglajkar V R.Nonreflecting boundary condition in finite-element formulation for an elastic wave equation.Geophysics，1998，63（3）：1006-1016.

[18]王月英，宋建國.波場正演模擬中Sarma邊界條件的改進(jìn).石油物探，2007，46（4）：359-362，389.Wang Y Y，Song J G.Improvement of Sarma boundary condition in wavefield forward modeling.Geophysical Prospecting for Petroleum （in Chinese），2007，46（4）：359-362，389.

[19]Yee K S.Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEE Transactions on Antennas and Propagation，1996，14（3）：302-307.

[20]Taflove A.Re-inventing electromagnetics：supercomputing solution of Maxwell＇s equations via direct time integration on space grids.Computing Systems in Engineering，1992，3（1）：153-168.

[21]Alvarez G B，Loula A F D，Dutra do Carmo E G，et al.A discontinuous finite element formulation for the Helmholtz equation.Computer Methods in Applied Mechanics and Engineering，2006，195（33-36）：4018-4035.

[22]徐世浙.地球物理中的有限單元法.北京：科學(xué)出版社，1994：260-277.Xu S Z.The Finite Element Method in Geophysics（in Chinese）.Beijing：Science Press，1994：260-277.

[23]Van der Vorst H A.Bi-CGSTAB：A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems.SIAM Journal on Scientific and Statistical Computing，1992，13（2）：631-644.

[24]柳建新，蔣鵬飛，童孝忠等.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用.中南大學(xué)學(xué)報(bào) （自然科學(xué)版），2009，40（2）：484-491.Liu J X，Jiang P F，Tong X Z，et al.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling.Journal of Central South University（Science and Technology）（in Chinese），2009，40（2）：484-491.

[25]Kriegsmann G， Morawetz C S.Numerical solutions of exterior problems with the reduced wave equation.Journal of Computational Physics，1978，28（2）：181-197.

[26]Bayliss A，Turkel E.Boundary conditions for the helmholtz equation in duct-like geometries.International Association for Mathematics &Computers in Simulation，1983，455-458.

[27]Engquist B，Majda A.Absorbing boundary conditions for numerical simulation of waves.ProceedingsoftheNational AcademyofSciencesoftheUnitedStatesofAmerica，1977，74（5）：1765-1766.

[28]Mur G.Absorbing boundary conditions for the finitedifference approximation of the time-domain electromagneticfield equations.IEEE Transactions on Electromagnetic Compatibility，1981，23（4）：377-382.

[29]Mei K K，F(xiàn)ang J Y.Superabsorption-a method to improve absorbing boundary conditions.IEEE Transactions on Antennas and Propagation，1992，40（9）：1001-1010.

[30]Gedney S D，Liu G，Roden J A，et al.Perfectly matched layer media with CFS for an unconditionally stable ADI-FDTD method.Antennas and Propagation IEEE Transactions on，2001，49（11）：1554-1559.

[31]Berenger J P.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，1994，114（2）：185-200.

[32]薛東川，王尚旭，焦淑靜.起伏地表復(fù)雜介質(zhì)波動(dòng)方程有限元數(shù)值模擬方法.地球物理學(xué)進(jìn)展，2007，22（2）：522-529.Xue D C，Wang S X，Jiao S J.Wave equation finite-element modeling including rugged topography and complicated medium.Progress in Geophysics（in Chinese），2007，22（2）：522-529.

[33]Claerbout J F.Imaging the Earth′s Interior.Boston：Blackwell Scientific Publications，1985.

[34]Semblat J F，Gandomzadeh A L L.A simple numerical absorbing layer method in elastodynamics.Comptes Rendus Mécanique，2010，338（1）：24-32.

[35]Caughey T K.Classical normal modes in damped linear dynamic systems.Journal of Applied Mechanics，1960，27（2）：269-271.