[發明專利]自適應有限元GPR頻率域正演方法

有效

• 申請號: 201910476650.6
• 申請日: 2019-06-03
• 公開（公告）號: CN110210129B
• 公開（公告）日: 2021-05-11
• 發明（設計）人: 馮德山;丁思元;王珣
• 申請（專利權）人: 中南大學
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06F17/14
• 代理公司: 長沙永星專利商標事務所(普通合伙) 43001
• 代理人: 周詠;米中業

1.一種自適應有限元GPR頻率域正演方法，包括如下步驟：

S1.設定正演方法的參數；

S2.根據步驟S1設定的參數，計算當前網格的數值解；具體為采用如下步驟進行計算：

A.采用頻域有限單元算法進行計算網格的數值解；

根據電磁波傳播理論，雷達波在有耗媒質中的傳播遵循Maxwell方程組，其二維TM波的頻率域表示形式為：

其中E_z為頻域率電場強度V/m，J_z為頻率域電流密度A/m²，μ為磁導率H/m，ε為介電常數F/m，σ為電導率S/m，ω＝2πf為角頻率rad/s，e^iωt為時諧因子，i為虛數單位；

假定介電常數ε和電導率σ是實數，再引入表示波的傳播和衰減的復介電常數，進一步簡化上式可得

考慮到探地雷達探測的目標對象多為非磁性介質，取其磁導率μ＝μ₀μ_r＝μ₀＝4π×10^-1H/m；μ_r＝1；已知自由空間的介電常數為ε₀＝8.854×10^-12F/m，相對介電常數表達式為ε_r＝ε/ε_r，則

設為自由空間波數，將頻率域波數定義為：

進一步改寫成為如下形式：

其中，μ_r＝1根據電磁場理論，將時諧場的各向異性介質的頻率域波動方程表示為：

其中，α_x＝α_y＝1/μ_r，β＝μ_r；

對上式進行積分處理，根據分部積分原理，該方程相應的弱解形式為：

其中，為權函數，而由于計算邊界為完美電導，二維TM波的頻率域表示式的右邊為0，因此，弱解形式可簡化為：

計算區域劃分單元后，積分可以寫成各個單元積分之和：

其中，e表示單元標號，N_e為單元總個數，Ω^e表示在第e個單元上的積分；若采用三角形單元進行網格剖分離散，并將場值函數u在各單元上應用線性插值：

根據Galerkin法，令權函數i＝1,2,3，對于任意三角形單元，有如下單元積分：

假設介質參數在單元內為均勻的，與積分無關，則可改寫為：

其中為單元剛度矩陣，為單元質量矩陣，f_i^e為單元源向量：

再將單元場向量、激勵源及單元系數矩陣擴展成整體矩陣，有如下矩陣表達：

-Ku_u+Mu_u＝f

其中系統矩陣K，M，和源向量f分別由它們的單元級別矩陣K^e，M^e，和源向量f^e組合而成，u_u是模擬波場向量；

B.采用精確完美匹配層處理邊界條件；

引入完美匹配層邊界條件，在完美匹配層吸收邊界區域，式中a_x＝s_y/s_x，a_y＝s_x/s_y，β＝s_xs_y，s_x和s_y為坐標拉伸變量，則上式可具體表示為：

本發明使用精確完美匹配層，其吸收函數在垂直于完美匹配層厚度增加方向的積分無窮大，使得最終得到的完美匹配層邊界條件下的解與無界區域下的真實解一致：

∫σ_i(ξ)dξ＝+∞ i＝x,y

坐標拉伸變量定義為s_x＝1-iσ_x/ω，s_y＝1-iσ_y/ω，其中k為探地雷達波動方程頻率域波數，以x方向為例，電導率σ_x具體計算公式為：

其中L_PML為完美匹配層的厚度，l_x表示完美匹配層內部單元中心到模擬區域邊界的距離，γ為自由變量；

S3.計算局部誤差指示值；

S4.標記需細化的單元集合；

S5.對步驟S4標記的網格進行細化；

S6.在步驟S5進行細化后的網格基礎上，進行數值求解，從而得到最終的自適應有限元探地雷達正演結果。