探地雷達有限元邊界條件的理論公式

1gpr正演模擬實地雷達的正轉模擬（pgr）一直是該領域理論研究的熱點。通過分析典型地質模型的結果，我們可以加深對gdpr反射剖面的理解和理解，并對gdpr數據的解釋起到了指導作用。在epr正轉模擬中，時域有限差分法（fdtd）由于簡單靈活，被廣泛使用，但復雜的物理界面不適合有限單元法。由于有限單元法（fem）不需要計算內部邊界條件，因此不需要計算內部邊界條件。它適用于物理復雜問題的標準化過程，具有廣闊的適應性。在地震波FEM模擬得到快速發展的同時,一些學者將其引入到GPR正演模擬中:沈飚進行了GPR簡單模型正演;底青云和王妙月推導了含衰減項的探地雷達波有限元方程,實現了復雜形體的GPR波的FEM正演模擬和偏移處理;Fanning和Boothby將GPR有限元建模技術應用到拱橋檢測中;謝輝等基于二十節點等參單元的FEM對Pulse型GPR進行了模擬;Arias等將有限元技術應用到GPR考古成像中;Lu等利用離散時域Galerkin進行色散介質中的GPR正演.GPR正演模擬是用有限空間模擬半無限大空間,在模型邊界上必然產生邊界反射.因此,需要在截斷邊界處使用一種特殊計算方法,在保證邊界處場的計算精度的同時大大消除人工截斷邊界的反射波.如Cerjan邊界條件,它是在邊界區引入一個衰減層,讓波的能量在這個衰減層里呈現指數衰減,但是這種吸收邊界條件的選取比較困難,衰減層的厚度、衰減系數的大小都要憑經驗去估計.Sarma在邊界區域建立一個一定厚度的衰減層,根據材料性質估算恒定比例系數,使衰減層和介質區存在明顯的物性差異,從而產生反射波;王月英等提出了一種改進的Sarma邊界條件,通過在衰減層內增加一定厚度的過渡帶,過渡帶內的比例系數呈現線性變化,這就使介質區和衰減層之間的物性差異漸變,大大減弱了交界面處的反射,但也存在對阻尼比的選取不敏感、過渡帶和吸收層的交界面的角點處產生較強的繞射波的缺點,仍需尋求新的解決辦法來取得更好的吸收效果.本文通過對透射邊界條件和改進Sarma邊界條件的研究,提出一種混合邊界條件.通過在研究區域外加載一定厚度的吸收層,以改進的Sarma邊界條件對到達邊界處的GPR波能量進行吸收衰減,然后在最外邊加載透射邊界條件,將衰減過的探地雷達波透射出去,這種混合邊界條件集成了二者的優勢,對人工截斷邊界的反射波取得更好的吸收效果.2陣風場-電導率-磁場強度法Maxwell方程組描述了電磁場的運動學和動力學規律,由電磁波理論,高頻電磁波在媒質中的傳播規律也應服從Maxwell方程組,其方程組的微分形式可以表示為式中,E為電場強度(V/m),H為磁場強度(A/m),B為磁感應強度(T),D為電位移(C/m2),J為電流密度(A/m2),ρ為電荷密度(C/m3).本構方程:式中ε為介電常數(F/m),μ為磁導率(H/m),σ為電導率(S/m).把方程式(2)代入方程式(1a),并求旋度,得由于ε,μ,σ為坐標的函數,它們的時間導數可以忽略,用電場E取代(4)式中的A,則(4)式中右邊第二項為0,整理得假設將雷達波激勵源如電場源SE或磁場源SH,代入到(5)式中,有同理對(1b)式中第二式兩邊求旋度,可得下式:(6)、(7)兩式表明,磁場H和電場E及其分量滿足相同的微分方程.3《gpr波》5,把合成gpr波作為初始條件,確定傳播方程,即求單元傳播到下一步的傳播到下一步的傳播到下一步的傳播到下一步的傳播到下一步的傳播到下一步的求取求解GPR波動方程其實質是結合雷達波所要滿足的初、邊值條件求解偏微分方程.邊界條件分為兩種類型:Dirichlet邊界條件與Neumann邊界條件.變分法推導結果自動滿足Dirichlet邊界條件,即未知函數在邊界上有已知值;伽遼金(Galerkin)法推導結果自動滿足Neumann邊界條件,即未知函數的導數在邊界上有已知值.下面采用伽遼金法推導含衰減項的GPR有限元方程.第1步網格剖分將求解的二維區域剖分成矩形單元,如圖1所示.第2步線性插值設(x,y)是子單元的坐標,(x0,y0)是子單元中點的坐標,a、b是子單元的兩個邊長,(ξ,η)為母單元的坐標,兩個單元間的坐標變換關系為其中微分關系為單元形函數為雙線性插值形函數,可寫為其中ξi,ηi是點i(i=1,2,3,4)的坐標,形函數的分量可寫為第3步單元積分根據伽遼金法,將(6)式兩邊同時乘以δE,并求積分,有其中E=NTEe,N=(N1,N2,N3,N4)T,EeT=(E1,E2,E3,E4),N1,N2,N3,N4是單元上各節點的形函數,E1,E2,E3,E4是各節點的電場值,Ω為單元面積.則式(12)左邊第一項為其中Me為單元質量矩陣:式(12)左邊第二項為其中K′e為單元阻尼矩陣:式(12)左邊第三項為其中Ke=(kij)=(kji).將(11)式求得ξ或者η的微商,代入到(18)式積分即可得到單元剛度矩陣kij為式(12)右邊項為第4步總體合成根據式(13)、(15)、(17)、(20)得到單元積分:將單元列向量其中ND是節點總數.將4×4的系數矩陣Me、K′e和Ke擴展成ND×ND的矩陣M,K′和K,將列向量Se擴展成ND維列向量S,即其中,由于δET≠0,所以必有M¨E+K′E·+KE=SE.于是可以得到時空域GPR波的二維有限元方程,這是含有ND個元的ND個方程聯合的常微分方程組,可以寫成如下形式:式中,M為質量矩陣,K′為阻尼矩陣,K為剛度矩陣,為時間的二階導數.第5步解GPR有限元常微分方程解常微分方程組(24)或(25)前,要代入邊界條件,邊界條件的加載方法將在下節討論.先采用中心差分法解方程(24).對初始條件離散化得到把時間域[0,T]分為幾個相等步長Δt=T/n,把時刻t的微分方程記為用差商代替微分:中心差分法是一種條件穩定的計算方法,當時間步長Δt取的過大時,計算結果就會出現數值色散,為此本文采用Δt≤(ΔXmin/Vmax),式中時間步長Δt與最小單元邊長ΔXmin成正比,與最大媒質速度Vmax成反比.差分網格空間步長滿足穩定性條件Δx<(λmin/10),其中λmin為最小波長長度,Δx為單元的最小尺寸.將式(27)代入到式(26)中得到此式即為GPR有限元正演遞推公式.由于零時刻的均為0,所以根據(28)式可計算出Δt時刻的EΔt,然后依次遞推可得到所有時刻的E值.令則(28)式可簡化為Ax=B形式的線性方程組.對于GPR電磁波,當媒質為良導體或者近似于良導體時,σ較大,的作用不能忽略,這時,GPR波在這種媒質中傳播就會發生頻散和被媒質吸收的情況;當GPR頻率ω很高,即的作用幾乎可以忽略,此時式(26)中的阻尼矩陣K′可以忽略,(26)式變為純波動方程.求解不含衰減項的GPR波動方程時,由于系數矩陣A往往是大型的病態稀疏矩陣,其條件數很大,如果采取對其直接求逆,計算量巨大,其求解結果不可信,只有當系數矩陣A是對角陣時,其逆矩陣是對角元素求倒數,這種方法才適用,為此采取集中質量矩陣方法來處理M矩陣,即將每一行(或列)的元素都加到對角線元素上去,則式(14)可轉化為單元集中質量矩陣:由于總體集中質量矩陣也有對角陣,所以不含衰減項的GPR有限元方程的差分迭代公式可寫為在求解含有衰減項的波動方程(28)時,系數矩陣A是質量矩陣M和阻尼矩陣K′的線性組合,同樣可以采用單元集中質量矩陣求解.但本文中選用不完全LU分解預處理的BICGSTAB(BiconjugateGradientStabilized)線性方程組求解算法,該算法是一種高效迭代法,特別適用于求解條件數很大的病態線性方程組.4吸收邊界條件目前應用于電磁場正演模擬的邊界條件主要集中于FDTD相關算法,如輻射邊界條件、基于單行波方程的吸收邊界條件、超吸收邊界條件、完全匹配層(PML)吸收邊界條件等.而GPR有限元法正演模擬文獻較少,采用的吸收邊界條件多借鑒與地震波有限元正演模擬的吸收邊界條件.下面主要介紹透射邊界條件與Sarma吸收邊界條件,在此基礎上提出了混合邊界條件.4.1邊界單元系數透射邊界條件是通過直接模擬各種單向波動的共同運動學特征來建立的一種邊界條件,原理簡單易懂,算法容易實現.首先給出穿過人工邊界上一點并沿邊界外法線方向傳播的單向波動的一般表達式,并假定所有單向波動以同一人工波速沿法線方向從邊界透射出去,由此可得透射邊界條件公式:式中L為微分算子.在每個單元中有這里Ei為每個單元上第i個節點的電場值.殘余量R為利用伽遼金余量法,二維有限元方程可寫為其中Ae為單元面積,r是殘余量R的向量表示,由式(32)代入式(33)得到.由于內邊界的積分相互抵消,所以只需求解圖2所示的區域外邊界積分.在均勻各向同性介質中,根據Claerbout推導的旁軸近似,其下行波方程為左行波方程為右行波方程為由法向導數的定義知其中,nx和ny為邊界外法線的方向余弦.將式(31)、(32)、(33)代入式(34)中,并利用高斯定理可得其中Γl、Γr、Γd分別表示左、右及底邊界,上邊界為自由邊界條件.將式(35)、(36)、(37)給出的邊界條件代入到式(39),可以得到:其中為電場對時間的一階導數,v為電磁波在媒質中的傳播速度.在加載透射邊界條件時,還需求解外法線的方向余弦,其求解示意如圖3所示.θx為GPR發射天線激勵源和邊界線單元x方向的夾角,θy為y方向夾角,無論這兩個角是正還是負,其余弦值nx和ny均為正值.在求解左邊界的邊界條件時,如圖4所示,令N1=(1-ξ)/2,N2=(1+ξ)/2,dΓ=l·dξ/2,這里l為邊界上的單元寬度,例如AB、CD邊l=b,BC邊l=a.先分析邊界單元在AB邊的情況:將邊界的單元系數矩陣定義為Fe,則左邊界的邊界單元系數矩陣擴展為Fle:同理,可以求得底邊界單元系數矩陣Fde和右邊界單元系數矩陣Fre:故邊界的阻尼矩陣F為將邊界阻尼加入到GPR有限元方程中去,得到如下式子:將(45)式透射邊界條件結合式(24)就可以進行GPR有限元正演模擬.4.2基于sa社會駁岸的民族反射波/文化邊界法Sarma吸收邊界條件是在求解空間外的區域增加一個由衰減性媒質組成的衰減層,如圖5a所示.對于衰減性媒質,阻尼力作用比較復雜,導致阻尼矩陣C不容易確定.Rayleigh在1877年給出了一種計算阻尼矩陣C的經典方法,該方法是利用整體質量矩陣M和整體剛度矩陣K來計算阻尼矩陣C:其中,式中ωi和ωj分別為介質的第i和第j個固有的圓頻率,ξi和ξj為對應的阻尼比.Caughey提出了一種更為普遍的阻尼矩陣計算公式:可以看出(46)式是(47)式的二階形式,即m=2,體現了質量矩陣和剛度矩陣的線性組合,并且這種組合非常容易計算,其單元中的計算公式如下:對于二階Caughey阻尼矩陣的兩個系數a0和a1與阻尼比ξ存在著如下的關系:其中ω為發射天線的中心頻率,通常0.05≤ξ≤0.30.將二階Caughey阻尼矩陣引入到GPR吸收邊界的處理中,得到吸收邊界區內的有限元方程為在衰減層內的比例系數a0和a1可由式(49)得到,這種方法雖然可行,但是衰減層和介質區還是存在著明顯的物性差異,這種差異必然會在交界面處產生人為的反射波.針對Sarma邊界條件在交界面處產生的人為反射波,王月英等提出了一種改進方法,如圖5b所示:即在衰減層內添加一個過渡帶,過渡帶內比例系數a0i和a1i由零逐漸增大到a0和a1,衰減層取一個半波長的厚度,過渡帶大約占半個波長.但該改進的Sarma邊界條件在過渡帶的四個角點位置仍存在著差異明顯的四個小邊界,為此本文采用一圈一圈的環形過渡帶加載辦法,消除了這種角點的人為差異,稱之為優化過渡帶角點的Sarma邊界條件,如圖5c所示.其中,衰減層內總節點數為n,過渡帶內的節點數為m,m<n,a0i和a1i(0≤i<m)為過渡帶內節點i處的比例系數,比例系數從交界面向衰減層邊緣方向呈線性增長關系,這種方法就減緩了交界面處的物性參數突變,從而減小交界面處產生的人為反射波能量,比例系數a0i和a1i增大到a0和a1就保持恒定,以確保GPR波在衰減層內能得到充分吸收.既能使交界面處的反射得到大大降低,又能讓邊界處的GPR波能量被充分吸收.4.3sa尉邊界狀態的優化設計通過對透射邊界條件和Sarma邊界條件的深入分析,提出了一種結合兩種邊界條件的混合邊界條件(如圖6所示),其主要思想是利用Sarma邊界條件對到達邊界區域的電磁波能量衰減功能和透射邊界對電磁波能量的透射功能,使GPR波經過Sarma邊界條件的衰減吸收后再通過透射邊界條件將剩余能量透射出去,集成了二者的優勢.處理辦法是在計算區域外面加一吸收層,用改進的Sarma邊界條件使GPR波能量被吸收減弱,然后在最外邊加載透射邊界條件.通過改變公式(48)中給的阻尼比ξ,可以實現調節Sarma邊界的吸收作用,從而和透射邊界條件相互協作,達到最佳吸收效果.5不帶邊界條件的雷達正演運用Matlab編制的有限元GPR模擬程序對10.0m×5.0m均勻模型進行了模擬,介質的介電常數ε=3.0,電導率σ=0.001S·m-1.網格單元大小為0.1m的正方形,GPR脈沖激勵源的子波形式為f(t)=t2e-atsinω0t,其中ω0為發射天線中心頻率為100MHz,電磁波脈沖的衰減速度取決于系數α,此處取α=0.93ω0.得到了半空間圖7所示的不帶邊界條件的雷達正演wiggle圖.圖中可見,A、B為左右邊界的反射波,C為底邊界的反射波,D、E為左下角和右下角的兩個角點的繞射波.這些人為截斷邊界的反射強度很大,對觀測區域影響嚴重,說明了對截斷邊界的處理勢在必行.5.1不同邊界條件下的反射波對比為了更好地說明不同吸收邊界條件對人工截斷邊界反射波的吸收效果,建立圖8所示9.0m×9.0m大小的正方形均勻模型,并把脈沖激勵源置于模擬區域的正中心,采樣時間間隔為0.5ns,其它參數與無邊界條件情況下一致.通過截取不同時刻的波場快照圖,觀測不同邊界條件對人工截斷邊界的吸收效果.圖9(a,b,c)分別為不帶邊界條件的20、30ns、35ns波場快照,圖中可見在30ns時刻GPR波前開始傳到區域邊緣,人為截斷邊界反射波開始形成,從35ns快照圖可以看出,人為截斷邊界反射回的GPR波能量很強,嚴重影響到對目標區域的研究,同樣說明,截斷邊界的處理是非常必要的.圖10(a,b,c)分別為模型邊界處加載透射邊界的30、35、40ns波場快照.圖中可見,在30ns時,波的能量剛剛傳播到邊界區域,在35ns時可以看出電磁波大部分的能量都從邊界處透射出去,只有少部分的能量反射回來,對比圖9c相同時刻的截斷邊界處的強反射波能量,說明了透射邊界條件取得了顯著效果.但是再分析圖10c中40ns的波場快照圖仍然可以看到截斷邊界的反射回波,說明透射邊界條件仍有待改進.5.2反射和角點繞射仍以上例模型和相應參數為例,分別以Sarma邊界條件、改進的Sarma邊界條件和優化了過渡帶的Sarma邊界條件對截斷邊界進行了處理.圖11a為沒有加載過渡帶的Sarma邊界條件的30ns波場快照圖,由于沒有加載過渡帶,在介質區和衰減層之間的界面產生了明顯的反射和角點繞射;圖11b為改進Sarma邊界條件的30ns波場快照圖,由于加載了過渡帶,在介質區和邊界吸收層之間的媒質特性是漸變的,明顯的減弱了介質之間的差異,降低了介質區和邊界吸收層之間界面的反射波能量,較好地處理了介質區和衰減層之間物性差異引起的反射波,但在過渡帶的四個角點位置存在著差異明顯的四個小邊界.為了消除這種因過渡帶人為加載方式引起的反射,本文采用一圈一圈的環形過渡帶加載辦法,圖11c為優化過渡帶Sarma邊界條件的30ns波場快照圖.它是在改進的Sarma邊界條件基礎上,對過渡帶的加載方法進行了優化,目的是減少人為的角點引起的異常,雖然效果不是很明顯,但是從側面反映出改進的Sarma邊界條件的過渡帶的角點異常主要是由角點引起,過渡帶的加載方法并非主要影響因素,重點工作應該放在研究過渡帶的角點異常消除上.5.3混合邊界邊界條件如圖12a為加載了透射邊界條件的35ns快照圖,圖12b為加載了優化過渡帶的Sarma邊界條件的35ns快照圖,圖12c為加載混合邊界條件的35ns快照圖.圖中可見,圖12b優化過渡帶的Sarma邊界條件由于過渡帶的四個角點區的繞射波能量比較強,其對GPR波的吸收效果并不是很理想,而圖12c中的混合邊界條件由于集成了Sarma邊界條件和透射邊界條件的優點,先衰減了到達邊界區域的GPR波能量,然后再將其透射出去,觀察相同時刻快照圖,可以看出,與圖9c無邊界條件的波場快照圖對比,各種邊界條件對截斷邊界處的反射波都取得了良好效果,但混合邊界條件明顯優于單純的透射邊界條件或Sarma邊界條件.6雷達檢測模擬的示例6.1混合邊界條件的gpr正演模擬如圖13所示,模型為一個10.0m×5.0m的矩形區域,分為上下兩層介質,中間是起伏界面,上層介質的相對介電常數ε1為3.0,電導率σ1為0.001S·m-1,下層介質的相對介電常數ε2為20.0,電導率σ2為0.01S·m-1,整個區域由單元邊長0.1m的正方形,網格剖分為100×50的網格空間,GPR波脈沖激勵源的中心頻率為100MHz,采樣時窗長度為100ns,采樣時間間隔為0.5ns.應用本文的基于混合邊界條件的FEM算法對起伏界面模型進行了正演模擬,得到了圖14所示的GPR正演模擬wiggle圖,圖中可見,在30ns附近有一條能量較強的反射界面,通過計算可以得出它與圖13模型中的上、下兩層分界面位置相