馮德山，陳佳維，吳奇

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

2北京大學地球與空間科學學院，北京 100871

3Rice大學地球科學系，德克薩斯州，休斯頓 77005-1892

1 引言

探地雷達（Ground penetrating radar，GPR）是根據(jù)脈沖電磁波在地下不同介質之間的反射及繞射等波動規(guī)律，來探測地下結構和特性的重要淺部地球物理方法（曾昭發(fā)等，2006）.GPR正演是計算地球物理學的一個重要分支，許多學者在此領域開展了大量研究：李展輝等（2009）將錯格時域偽譜法應用于GPR模擬，解決了時域偽譜法的Gibbs現(xiàn)象；馮德山等（2007）應用時域多分辨算法開展GPR數(shù)值模擬，其每個波長僅采樣2個點的精度可達到時域有限差分法（Finite difference time domain，F(xiàn)DTD）每個波長采樣10個點的精度；方宏遠等（2012，2013）利用Hamilton系統(tǒng)的辛分塊龍格庫塔方法開展了道路結構層雷達波傳播特征研究；馮德山等（2013）將只需節(jié)點無需單元的無網格Galerkin算法應用于GPR正演，它具有前處理簡單、解高次連續(xù)等優(yōu)點；Lu等（2005）利用離散時域Galerkin進行色散介質中的GPR正演；沈飚（1994）應用有限單元法（FEM）進行了GPR簡單模型正演；底青云與王妙月（1999）推導了含衰減項的GPR波有限元方程，實現(xiàn)了復雜形體的GPR波的FEM正演；馮德山等（2012）提出結合透射邊界和Sarma邊界的混合邊界，改善了有限單元法GPR正演的邊界反射.自Yee（1966）提出FDTD算法以來，F(xiàn)DTD算法在GPR正演中得到廣泛應用：Giannopoulos（2005）編制了GPR正演軟件“GprMax”，因為是免費軟件，已得到許多用戶青睞；Irving和Knight（2006）基于Matlab開發(fā)了GPR二維正演程序；李靜等（2010）為了提高FDTD算法的模擬精度，采用高階FDTD算法結合單軸各向異性理想匹配層（UPML）邊界條件進行了GPR數(shù)值模擬；馮德山等（2010）應用基于UPML邊界條件的交替方向隱式時域有限差分法（ADI-FDTD）進行GPR正演；田鋼等（2011）將FDTD算法應用于地面以上物體反射干擾特征GPR模擬和分析，為干擾壓制研究提供理論依據(jù)；馮晅等（2011）通過水平正交雙方向同時接收獲取全極化信息，構建了基于FDTD算法的全極化GPR正演模擬方法，進而能更好地分析目標屬性.這些研究指導了異常體的判讀與GPR數(shù)據(jù)解釋，考慮到雷達波在地下介質傳播時，易產生波形畸變、能量衰減的頻散現(xiàn)象.Bano（2004）、Cassidy and Millington（2009）、Teixeira（2008）等強調了GPR正演時頻散特性的重要性；劉四新等（2007）開展了頻散介質中二維FDTD法GPR數(shù)值模擬，對比了頻散介質與非頻散介質中GPR曲線的區(qū)別；方廣有等（1996）開展了淺層地下目標電磁散射問題的FDTD數(shù)值模擬；Li等（2012）開展基于遞歸積分不分裂復頻移PML的FDTD三維GPR正演，該邊界條件對隱失波、低頻波及近掠射波具有好的吸收效果.

由此可見，GPR正演算法種類繁多，各有特色，且經歷了從不考慮介質頻散向考慮頻散的發(fā)展歷程，在這些模擬方法中，F(xiàn)DTD具有概念清晰、編程簡單等特點，成為GPR模擬中應用最廣泛的方法.常規(guī)FDTD開展GPR模擬時，常采用一致網格剖分，為了保證計算精度、避免數(shù)值頻散，當模型中存在局部物性參數(shù)變化劇烈小異常體時，需要采用統(tǒng)一的細網格剖分，精細的網格又要求小的時間步長；然而，除去這些局部區(qū)域，其他大部分區(qū)域波場變化都是平緩的，利用較粗的空間網格與大的時間步長就能進行正演.這樣，在波場緩變區(qū)域存在時間與空間上過采樣問題（黃超和董良國，2009a），造成計算機存儲空間加大和計算效率降低.為了兼顧FDTD模擬效率與精度，避免過采樣問題，Moczo（1989）首次提出網格連續(xù)變化的方法來減小計算量；Zivanovic等（1991）將變網格技術應用到電磁場Maxwell方程求解中，但僅限于空間一個方向，由于統(tǒng)一的小時間步長的存在，計算效率仍有待提高；Falk等（1998）提出一種長短時間步長比為2N倍的局部可變時間步長模擬方法；Tessmer（2000）對Falk的方法進行了改進，使長短時間步長的變化范圍突破了2N的限制；黃超和董良國（2009b）提出一種可變網格與局部可變時間步長的高階差分聲波模擬方法，而后又將該算法引入到交錯網格高階差分彈性波模擬中；張慧和李振春（2011）基于時空雙變網格算法的碳酸鹽巖裂縫型儲層正演模擬.這些雙變步長FDTD算法解決了空間與時間維度上的過采樣問題，但計算過程中時間步長算子作為一個全局變量，局部變化非常困難，給計算機編程增加了不少難度，同時也可能因時、空域的雙插值帶來穩(wěn)定性問題.

本文繼承了Wang etal.（2001），Ahmedy和Chen（2004），Diamanti和Giannopoulos（2009）等人FDTD混合算法成果，提出基于Debye模型的復介電常數(shù)開展ADI-FDTD（Alternating-Direction-Implicit）亞網格技術的頻散介質GPR正演，即在物性變化劇烈局部區(qū)域采用細網格剖分ADI-FDTD計算，其他的物性參數(shù)緩變區(qū)域采用粗網格常規(guī)FDTD計算，由于ADI-FDTD突破了CFL條件的約束，可選取與粗網格一致的大時間步長，計算中只需考慮空間不一致網格的插值連接問題，編程計算簡便，從根本上降低內存占用及計算時間，在精確性與時效性之間取得滿意的平衡.

2 頻散介質混合ADI-FDTD亞網格技術原理及差分格式

如圖1a所示，當模擬區(qū)域存在小的異常體時，常規(guī)FDTD算法采用固定的網格步長對整個模型區(qū)域進行精細的網格離散易導致物性參數(shù)緩變區(qū)空間與時間上的過采樣.圖1b為局部亞網格的剖分方法，僅在異常體及其附近區(qū)域采用細網格剖分ADIFDTD計算，其他區(qū)域采用粗網格剖分常規(guī)FDTD計算，因為ADI-FDTD可選取與粗網格中FDTD算法一致的大時間步長，計算中只需考慮空間不一致網格的插值連接問題，時間上完全不用考慮插值問題，編程計算較雙變網格算法更為簡便.

圖1 不同的網格剖分方式（a）精細網格剖分方式；（b）局部亞網格剖分（僅異常體及附近區(qū)域采用細網格）.Fig.1 Different grid splid methods

由電磁波傳播理論可知（Yee，1996），GPR遵循的Maxwell方程在二維情況下變?yōu)?個偏微分方程，其電磁分量僅為Ez、Hx、Hy三個，稱之為TM波.在各向同性介質中D（t）＝ε（r）E（t），介電常數(shù)ε（r）只是空間的函數(shù)，可直接將其代入Maxwell旋度方程的離散式.而在頻散介質中介電常數(shù)是關于空間與頻率的函數(shù)，頻散介質中FDTD遞推公式與各向同性介質中FDTD遞推公式類似，僅僅是電磁分量的推導稍顯不同（Luebbers etal.，1990；Luebbers etal.，1991）.處理頻散介質FDTD算法有許多種（劉少斌等，2010），包括遞歸卷積積分法、輔助方程法、Z變換法、位移算子法等.本文應用輔助差分法（葛德彪和閆玉波，2005）推導頻散介質中的FDTD差分方程，無源情況下頻散介質中Maxwell方程表示為

其中，E為電場強度（V·m-1）；H為磁場強度（A·m-1）；B為磁通密度（Wb·m-2）；σm為電通密度（C·m-2）；μ為磁導率（H·m-1）；ε為介電常數(shù)（F·m-1）.

在粗網格區(qū)域，采用傳統(tǒng)FDTD方法.運用Yee氏網格模型，把連續(xù)變量離散化，Δx、Δy為圖1b中粗網格區(qū)域的空間步長，Δt為差分時間步長，將每個節(jié)點進行編號，例如：nΔt時刻空間坐標位置（iΔx，jΔy）的電磁場值f（iΔx，jΔy，nΔt）按fn（i，j）方式對應起來.考慮到頻散介質中，磁場分量的差分公式與非頻散介質中是一致的，將（1）離散，并整理可得到磁場分量Hx與Hy在t＝nΔt時刻時間迭代FDTD公式：

式（4）、（5）中的電磁參數(shù)中標號m的取值與其式中右端電場或磁場分量的空間位置相同，式中

將式（2）在t＝ （n+1／2）Δt時刻做差分離散：

即

表達介質頻散的模型主要有Lorentz（洛淪茲）模型、Drude（德魯）模型、Debye（德拜）模型（劉少斌等，2010）.下面根據(jù)Debye模型關系式，推導頻散介質中的D和E的微分方程.可將有耗頻散媒質的介電常數(shù)表示為

其中εs為靜態(tài)介電常量；ε∞為無限頻率的介電常量；τ0為介質極化的弛豫時間常數(shù).將式（11）代入到式（3）可得在頻域中有

即

利用輔助差分法中頻域到時域算子轉換關系jω＝?／?t，將（13）式變成時域方程：

式（14）為與Debye關系式（11）對應的時域中Dz和Ez的微分方程.為推導由Dz到Ez的FDTD差分公式，將（14）式在t＝ （n+1／2）Δt時刻離散得到：

即

由此可得到時域中由Dz到Ez的轉換差分公式（李慶偉，2008）：

公式（17）作空間離散后可得到Ez分量在時域的推進式.所以頻散介質FDTD的隨時間推進的步驟為

（1）根據(jù)常規(guī)的FDTD差分公式（4）—（5），由E計算Hn+1／2的各個分量；

（2）根據(jù)（10）式由H計算Dn+1的各分量；

（3）根據(jù)（17）式由D計算En+1的各分量；

（4）回到步驟（1）.

在異常體及其附近細網格剖分區(qū)域中，使用ADI-FDTD方法將傳統(tǒng)的一個時間步分成兩個子時間步，分別采用前向和后向差分，兼具隱式差分格式的無條件穩(wěn)定性和顯式差分格式計算相對簡單的優(yōu)點（Namiki，1999）.ADI-FDTD方法中兩個時間步的差分公式推導類似，以n→n+1／2步為例，n+1／2→n+1應用同樣的方法導出.再分析Maxwell旋度方程（1）與（2），式中時間偏微分項仍采用中心差分格式，Δx′、Δy′為圖1b中細網格區(qū)域的空間步長，對x方向導數(shù)在子時間步的末時刻 （n+1／2）Δt取隱式差分離散，對y方向導數(shù)在初時刻nΔt作顯式差分離散，Hx分量在n和n+1／2時刻均取值，并定義：

同時，定義標號 （m）＝ （i，j+1／2），則式（1）中的Hx分量差分格式為

同理，如果定義標號 （m）＝ （i+1／2，j），則式（1）的Hy分量差分格式為

再分析（9）式，x方向導數(shù)的差分離散取隱式，y方向導數(shù)的差分離散取顯式，則式（9）的差分格式為

式（14）為與Debye關系式（11）對應的時域中Dz和Ez的微分方程.為推導由Dz到Ez的ADIFDTD差分公式，將（14）式在n→n+1／2步中離散得到：

由此可得到時域中由D到E的轉換關系：

其中

式（21）—（23）及（25）四個式子可用于子時間步n→n+1／2電磁場的時域推時計算，其中式（21）在計算n+1／2時刻場值時其右端只涉及n時刻場值，稱為顯式格式，但式（22）等式兩邊包含同時刻和場分量，式（23）等式兩邊包含同時刻和場分量，式（25）等式兩邊包含同時刻和場分量，不能構成顯式時間推進計算，邏輯上講時間步進是無法完成的，稱為隱式格式.為此，將式（22）代入到式（23）中消除磁場分量的，再將得到的消除了的式（23）代入到式（25）消去，得：

多元方程組式（28）的系數(shù)矩陣可以改組成實三對角條帶矩陣，結合吸收邊界條件就可求解（鄭奎松，2005）.

3 粗細網格場值交換格式

如圖2所示，ADI-FDTD亞網格算法空間網格的布置與常規(guī)的FDTD算法是一致的，電場位于網格的中心垂直向里，磁場分別位于網格的邊上，粗細網格比為3∶1，圖2中紅色箭頭表示粗網格磁場值、灰色小箭頭表示需要插值計算的細網格磁場值，粗網格中應用FDTD算法，細網格中應用ADI-FDTD算法，進行電磁場時間步推進時，在大小網格的交界處，細網格的場值更新的部分數(shù)據(jù)在粗網格上是缺失的，這就需要用到插值技術，粗細網格之間過渡與銜接、場值的交換與更新是本算法的重要內容.

粗細網格的場值交換格式有多種，可以采取粗細網格邊界純磁場交換格式、純電場交換格式、電磁場結合交換格式（Diamanti，2008）.此外，在選取場值交換時也有細微的區(qū)別，可以僅選取傳遞粗細網格邊界上的磁場值，也可以選取粗細網格邊界附近的一個粗網格上的電磁場值.本文采用電磁場結合的交換格式，其電磁場交換格式如圖3所示，圖中粗網格FDTD計算出的磁場值Hx、Hy置于每個網格邊上的中間位置，傳遞進細網格hx、hy，紅色的細箭頭表示需要插值計算的磁場值，ADI-FDTD細網格計算完成以后，再將圖3中粗細網格邊界向細網格內的0.5個粗網格處紅色標記的電場值ez傳回給FDTD粗網格Ez，該算法的詳細計算過程見流程圖4.根據(jù)Chevalier etal.（1997）與Diamanti（2008）的研究，應用亞網格算法進行電磁場計算中，粗細網格的比最好不超過11，為了避免電磁場的突然改變，物性介質分界面不要置于粗細網格的邊界，而異常體的每邊距離粗細網格的分界面最好有3個粗網格的距離.此外，計算過程中為了粗網格的FDTD與細網格的ADI-FDTD計算時間上的同步，將ADIFDTD的時間步長設為粗網格FDTD時間步長的一半，這樣在每計算一個FDTD時間步長，ADIFDTD算法調用兩次，兩次加起來剛好與FDTD的一個時間步長相等，這樣設置比將FDTD與ADIFDTD采用同樣的時間步長精度更高.

分析圖2、3可知，設置粗細網格比為3∶1的奇數(shù)倍，可將粗網格的電磁場值直接賦予細網格的磁場值，方便粗細網格間的電磁場值的數(shù)據(jù)交換.當采用圖5所示的粗細網格比為4∶1的偶數(shù)倍時，則粗網格的磁場值并不在某一小格的邊上，需要進行插值計算，而ADI-FDTD細網格計算完成以后，電場值ez傳遞回FDTD粗網格Ez時（圖中的綠色大圓點），需要周圍的4個相鄰的ez場值的平均，給計算帶來了一定的不便，但理論上是可行的，因此，粗細網格比理論上可以為任意倍數(shù).而圖3中的粗細網格界面處的磁場值（紅色小箭頭表示的），可以采用簡單的線性插值計算，如圖6所示，以y方向磁場

圖2 ADI-FDTD亞網格技術粗細網格及場值分布示意圖Fig.2 The field distributing sketch map of coarse grid and fine grid on ADI-FDTD subgridding method

圖3 粗細網格場值交換格式Fig.3 The tranfer mode of field value between on coarse grid and fine grid

為例：

本文僅需式（33）中2個插值公式即可.如果采用粗細網格邊界1個粗網格內的場值交換方式，還需計算：

4 正演模擬實例

4.1 模型一

圖4 ADI-FDTD亞網格技術計算流程圖Fig.4 The flow chart of ADI-FDTD subgridding formulation computation

圖5 粗細網格比為4倍時的網格及場值分布示意圖Fig.5 Sketch map of coarse grid to fine grid ratio equal to 4

圖6 磁場線性插值計算示意圖Fig.6 Sketch map of line interpolation

圖7 薄矩狀異常體地電模型及亞網格剖分示意圖Fig.7 The sketch map of thin rectangle anomalous body model and subgridding splitting

模型一如圖7所示，模型區(qū)域3.0m×2.0m，背景介質為混凝土，其介電常數(shù)6.0，電導率0.0005S·m-1，混凝土介質中有1個埋深0.5m，大小0.4m×0.1m薄矩狀素填土缺陷，其介電常數(shù)10.0，電導率0.002S·m-1，對角點坐標（1.3m，0.6m）、（1.7m，0.5m）.正演中UPML邊界設為8層粗網格，Ricker子波距邊界為10層網格，主頻900MHz.為了說明ADI-FDTD亞網格技術的計算效率及精度，將該算法與均勻剖分的粗網格FDTD、細網格FDTD算法進行對比.計算環(huán)境為Intel（R）CoreTMi5CPU，P8600＠2.67GHz，1.17GHz+2.86GB的內存物理地址擴展，Window XP操作系統(tǒng)的聯(lián)想（IBM）X201筆記本.一致粗網格模擬步長Δx＝Δy＝0.009m，粗網格數(shù)332×221，占用內存1.84Mb，時間步長Δt＝2.123×10-11s，時窗長度取12ns，每隔0.05m采集一道雷達數(shù)據(jù)，計算50道雷達數(shù)據(jù)，共計花費時間1′24″.空間一致細網格模擬步長Δx′＝Δy′＝0.003m，細網格數(shù)999×666，占用內存11.92Mb，時間步長Δt′＝0.826×10-11s，時窗長度同樣設為12ns，模擬50道雷達數(shù)據(jù)，共計花費時間21′11″.應用異常體及其周圍區(qū)域細網格剖分，其他區(qū)域采用粗網格剖分的ADI-FDTD亞網格技術計算，粗細網格比為3∶1，細網格區(qū)域大小設為0.81m×0.42m，細網格數(shù)270×140，其余區(qū)域設為粗網格，粗細網格時間步長統(tǒng)一取Δt＝2.123×10-11，由于細網格區(qū)域的ADI-FDTD算法分2步進行，取它的時間步長為FDTD時間步長的一半，該算法占用內存2.38Mb，共計花費時間2′46″.3種算法的計算資源占用情況詳見表1.

表1 模型一3種算法計算資源占用詳情Table 1 The occupation detail of calulate resource with three methods for model 1

圖8為3種方法模擬所得的單道雷達記錄，圖8a中紅色曲線代表粗網格方法，綠色的星號代表細網格方法，黑色曲線代表ADI-FDTD亞網格技術，圖中可見，3種算法結果基本吻合，一致性比較好，只是在直達波的位置，粗網格的Ez負場值較另外兩種方法稍大，8.2ns與10.4ns附近的正負波峰很好地對應了薄矩狀異常體的上下界面；圖8b為截取圖8a中7～12ns時刻放大后的雷達記錄，圖中黑色的曲線（ADI-FDTD亞網格技術）與綠色的曲線（細網格）吻合得更好，紅色曲線（粗網格）與另外兩條曲線稍有偏差.如果認為細網格與計算結果與解析解接近，說明ADI-FDTD亞網格技術在使用較少內存、較短運行時間的前提下，計算結果具有較高的精度，在運算時間與計算精度方面達到了較好的平衡.

4.2 模型二

模型二設計為了對比頻散介質與非頻散介質的雷達波傳播特征，如圖9所示，該模型為3.0m×2.0m的矩形區(qū)域，左邊非均勻背景媒質的相對介電常數(shù)εr＝6.0，電導率σ1＝0.002s·m-1，右邊頻散媒質采用Debye模型表示，ε∞r＝10.28，ε0r＝19.0，電導率σ1＝0.002S·m-1，馳豫時間τ5ns；在背景媒質中有3個異常體，其中2個圓形異常體具有相同的半徑r＝0.1m，左邊圓形異常體為圓心在（1.2m，0.5m）處的金屬介質，右邊的圓形異常體其圓心在（1.8m，0.5m），相對介電常數(shù)εr＝3.0，電導率σ1＝0.0001S·m-1.在2個圓形異常體的下方為矩形異常體，矩形異常體的長與寬均為0.3m，它距離左邊界1.35m，距離上邊界0.85m，圓形異常體的相對介電常數(shù)εr＝6.0，電導率σ1＝0.005S·m-1.Ricker子波主頻為500MHz，UPML邊界與源的位置與上例一致.

一致粗網格空間步長Δx＝Δy＝0.015m，占用內存0.68Mb，時間步長Δt＝3.538×10-11s，時窗長度取25ns，采樣間隔為0.015m，與上例同等計算條件下計算190道雷達數(shù)據(jù)，共計花費時間2′30″.細網格空間步長Δx′＝Δy′＝0.005m，內存占用5.88Mb，時間步長Δt′＝1.1793×10-11s，時窗長度取25ns，采樣間隔為0.015m，計算190道雷達數(shù)據(jù)，共計花費時間45′36″.應用本文的基于亞網格技術的ADI-FDTD算法，只在異常體及其周圍區(qū)域采用細網格ADI-FDTD計算，其他區(qū)域采用粗網格FDTD計算，粗細網格比為3∶1，細網格區(qū)域大小設為1.0m×1.0m，時間步長設為Δt＝3.538×10-11，ADI-FDTD亞網格技術的算法占用內存為1.55Mb，共計花費時間9′41″.模型二3種算法計算資源占用情況詳見表2.

表2 模型二3種算法計算資源占用詳情Table 2 The occupation detail of calulate resource with three methods for model 2

圖8 模型一3種方法計算的單道雷達數(shù)據(jù)（a）0～12ns的雷達數(shù)據(jù)；（b）7～12ns放大后的雷達數(shù)據(jù).Fig.8 The map of single track radar data of model one with three methods

圖9 模型二示意圖Fig.9 The sketch map of model two

圖10 均勻介質常規(guī)FDTD計算的GPR掃描圖Fig.10 GPR map of homogeneous medium model

圖11 不均勻介質常規(guī)FDTD計算的GPR掃描圖Fig.11 GPR map of nonuniformity model

圖12 去掉模型二中異常體的不均勻介質GPR掃描圖Fig.12 GPR map of remove the anomalous bodies of model 2

圖13 模型二一致細網格GPR掃描圖Fig.13 GPR map of model 2with fine mesh

圖14 模型二一致粗網格GPR掃描圖Fig.14 GPR map of model 2with coarse mesh

圖15 ADI-FDTD亞網格技術GPR掃描圖Fig.15 GPR map of model 2with ADI-FDTD subgridding mesh

模型二的GPR正演剖面比較復雜，為了更好地理解正演圖中的波形特征，采用在模型中分步加入異常體的方法，從簡到繁，便于理解.圖10假定模擬區(qū)域充填圖9中左邊均勻非頻散介質的情況，該均勻介質雷達剖面圖中因沒有異常體的存在故不存在反射波，只有直達波與場源第二次到達時的波峰（簡稱場源二次波），該處顯示的場源二次波及其后面的場源多次波是由于增益分段放大顯示，與作者將數(shù)據(jù)寫入GSSI軟件顯示方式有關；圖11為左右兩邊存在不均勻非頻散介質，左邊介質與圖9中左邊介質一致，右邊介質相對介電常數(shù)為19.0，電導率0.002S·m-1.圖中可見，由于介質分界面的存在，天線在左邊介質時逐漸靠近界面時，分界面反射雙程走時越來越短，最后靠近分界面，雙程走時為零，故圖中反射曲線為一條斜線，而當天線處于右邊介質中時，雷達天線逐漸遠離介質分界面，故分界面反射雙程走時越來越長，由于天線兩邊介質雷達波傳播速度的不同導致了兩條強反射界面斜率的不一致，右邊介質介電常數(shù)大，故傳播速度小，導致斜率大.此外，介質分界面的二次反射波及場源的多次波峰圖中也清晰可辨.

圖12為去掉模型二中3個異常體所得的正演剖面圖，由圖12中可見，由于右邊頻散媒質馳張時間及介電常數(shù)的差異，導致原本同一波形1與2，出現(xiàn)明顯的錯位，再利用該模型進行波場快照，將源與接收位置都置于模型中心，選取圖16中的7ns與9ns波場快照，分析可知，隨著時間的推進，雷達波以球面形式向外擴散能量以指數(shù)方式衰減，離場源中心越遠則能量越弱，由于電磁波在左邊介質傳播得較右邊更快，距離場源更遠，能量應該更弱，但兩幅圖中Ez波場幅值與該結論明顯相反，左邊的波場幅值較右邊的更大，顯然這是由于右邊頻散介質對雷達波吸收較厲害，導致同時刻右邊頻散介質的波場能量較左邊更弱.

圖16 模型2去掉異常體不同時刻波場快照3D顯示（a）7ns波場快照；（b）9ns波場快照.Fig.16 The 3Ddisplay of wavefield snapshot in difference moment of model 2by remove the anomalous bodies

圖17 模型二不同時刻波場快照Fig.17 The wavefield snapshot in difference moment of model 2

圖13 、14、15為分別運用一致細網格FDTD、一致粗網格FDTD及其ADI-FDTD亞網格技術對模型二模擬所得的GPR掃描圖，圖中波形3為左邊圓形金屬異常體上界面的反射波、6為右邊圓形異常體的上界面，4與5分別為矩狀異常體的上、下界面，3種算法對異常體都能有較好的反映，其中一致細網格算法的波場信息最為豐富，對細微結構反映最好，ADI-FDTD亞網格技術其次，粗網格算法在局部波形上甚至出現(xiàn)細微頻散、不光滑現(xiàn)象，而ADI-FDTD亞網格方法盡管也出現(xiàn)了細微頻散現(xiàn)象，但細節(jié)信息較一致粗網格FDTD算法更豐富、局部波形更光滑.綜合考慮占用內存、模擬時間及數(shù)據(jù)精度，基于ADI-FDTD亞網格方法達到了一個較好的均衡.

為了能更直觀地展示雷達波在地下媒質中的傳播過程，將場源與接收點都置于模型中心，然后利用ADI-FDTD亞網格算法模擬模型2并繪制了圖17所示的5ns、6ns、7ns、8ns時的波場快照.圖中顏色的強弱代表電場幅值的大小.圖17a中，GPR波在兩邊介質中以不同速度呈半圓形擴散，當遇到矩狀異常體的4條邊時，波形大部分能量繼續(xù)向外傳播，一部分能量被4條邊反射回來；圖17b中，由于左右兩邊為不同介質，左邊介質中雷達波傳播較快，右邊介質速度傳播較慢，所以兩邊介質中圓的大小顯然不一致，傳播速度較快的左邊介質，雷達波前傳得更快、圓的半徑更大，此時傳播更快的左邊介質中，波前剛好到達了圓形金屬異常體的位置，并開始產生繞射波；圖17c中，GPR波繼續(xù)向外擴散，速度較慢介質中雷達波前到達右邊圓形異常體，剛開始產生繞射，而左邊的圓形金屬異常體、及其中間矩形異常體產生的繞射波與原波前以相反方向呈圓形擴散；圖17d中清晰可見GPR波波前到達的位置，由于介質兩邊雷達波傳播速度的不一致，左邊波前傳播更快、左邊圓形金屬異常體繞射波產生的圓弧更大，而右邊的異常體也產生了明顯的繞射波，只是該繞射波產生的圓弧更小些.圖中的各異常體反射波、繞射波清晰可辨.

5 結論

（1）提出了基于混合ADI-FDTD亞網格技術GPR正演算法，即在物性參數(shù)變化劇烈局部區(qū)域采用細網格剖分ADI-FDTD計算，其他的區(qū)域采用常規(guī)FDTD計算，因為ADI-FDTD突破了CFL條件的約束，可選取與粗網格一致的大時間步長，計算中只需考慮空間不一致網格的場值交換問題，時間步長一致完全不用考慮，編程計算較雙變步長算法簡便，有效地提高了計算效率.

（2）該混合ADI-FDTD亞網格技術能較好地適用于局部細微結構的模擬，在保證計算精度前提下有效減少內存占用，在精確性與時效性之間取得滿意的平衡.該算法被用于頻散介質及非頻散介質的組合模型的正演計算中，雷達正演剖面圖及波場快照結果都表明，考慮介質對雷達波的衰減的正演計算結果更符合雷達波實際地下的傳播特征，能更有效指導工程資料的精確解釋.

Ahmed I，Chen Z Z.2004.A hybrid ADI-FDTD subgridding scheme for efficient electromagnetic computation.International Journal of Numerical Modelling：Electronic Networks，Devices and Fields，17（3）：237-249.

Bano M.2004.Modelling of GPR waves for lossy media obeying a complex power law of frequency for dielectric permittivity.Geophysical Prospecting，52（1）：11-26.

Cassidy N J，Millington T M.2009.The application of finitedifference time-domain modelling for the assessment of GPR in magnetically lossy materials.Journal of Applied Geophysics，67（4）：296-308.

Chevalier M W，Luebbers R L，Cable V P.1997.FDTD local grid with material traverse.IEEE Transactions on Antennas and Propagation，45（3）：411-421.

Di Q Y，Wang M Y.1999.2Dfinite element modeling for radar wave.Chinese J.Geophys.（in Chinese），42（6）：818-825.

Diamanti N，Giannopoulos A.2009.Implementation of ADI-FDTD subgrids in ground penetrating radar FDTD models.Journal of Applied Geophysics，67（4）：309-317.

Diamanti N.2008.An efficient Ground Penetrating Radar finite-Difference time-Domain Subgridding Scheme and Its Application to the Non-descructive Testing of Masonry arch Bridges［D.Ph.thesis］.Edinburgh：The University of Edinburgh.

Falk J，Tessmer E，Gajewski D.1998.Efficient finite-difference modelling of seismic waves using locally adjustable time steps.Geophysical Prospecting，46（6）：603-616.

Fang G Y，Zhang Z Z，Wang W B.1996.Numerical simulation of time-domain EM scattering by shallow subsurface objects.Journal of Electronics and Information Technology（in Chinese），18（3）：276-282.

Fang H Y，Lin G.2012.Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar.Computers ＆Geosciences，49：323-329.

Fang H Y，Lin G.2013.Simulation of GPR wave propagation in complicated geoelectric model using symplectic method.Chinese J.Geophys.（in Chinese），56（2）：653-659.

Feng D S，Chen C S，Dai Q W.2010.GPR numerical simulation of full wave field based on UPML boundary condition of ADIFDTD.Chinese J.Geophys.（in Chinese），53（10）：2484-2496.

Feng D S，Chen C S，Wang H H.2012.Finite element method GPR forward simulation based on mixed boundary condition.Chinese J.Geophys.（in Chinese），55（11）：3774-3785.

Feng D S，Dai Q W，Weng J B.2007.Application of multiresolution time domain method in three dimensions forward simulation of ground penetrating radar.Journal of Central South University：Science and Technology（in Chinese），38（5）：975-980.

Feng D S，Wang H H，Dai Q W.2013.Forward simulation of Ground Penetrating Radar based on the element-free Galerkin method.Chinese J.Geophys.（in Chinese），56（1）：298-308.

Feng X，Zou L L，Liu C，etal.2011.Forward modeling for full polarimetric ground penetrating radar.Chinese J.Geophys.（in Chinese），54（2）：349-357.

Ge D B，Yan Y B.2005.Finite-Difference Time-Domain Method for Electromagetic Waves（2nd ed）.（in Chinese）.Xi′an：Xidian University Press.

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

Huang C，Dong L G.2009.High-order finite-difference method in seismic wave simulation with variable grids and local timesteps.Chinese J.Geophys.（in Chinese），52（1）：176-186.

Huang C，Dong L G.2009.Staggered-grid high-order finitedifference method in elastic wave simulation with variable grids and local time steps.Chinese J.Geophys.（in Chinese），52（11）：2870-2878.

Irving J，Knight R.2006.Numerical modeling of groundpenetrating radar in 2-D using MATLAB.Computers ＆Geosciences，32（9）：1247-1258.

Li J，Zeng Z F，Huang L，Liu F S.2012.GPR simulation based on complex frequency shifted recursive integration PML boundary of 3Dhigh order FDTD.Computers ＆Geosciences，49：121-130.

Li J，Zeng Z F，Wu F S，etal.2010.Study of three dimension highorder FDTD simulation for GPR.Chinese J.Geophys（in Chinese），53（4）：974-981.

Li Q W.2008.The Propagation of GPR Pulse in the Lossy，Dispersive，and Inhomogeneous Earth Medium（in Chinese）［D.Ph.thesis］.Chengdu：Chengdu University of Technology.

Li Z H，Huang Q H，Wang Y B.2009.A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media.Chinese J.Geophys.（in Chinese），52（7）：1915-1922.

Liu S B，Liu S，Hong W.2010.Dispersive medium FDTD method（in Chinese）.Beijing：Science Press.

Liu S X，Zeng Z F.2007.Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.Chinese J.Geophys.（in Chinese），50（1）：320-326.

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

Luebbers R J，Hunsberger F P，Kunz K S，etal.1990.A frequency-dependent finite-difference time-domain formulation for dispersive materials.IEEE Trans.Electromagn.Compat.，32（3）：222-227.

Luebbers R J，Hunsberger F P，Kunz K S.1991.A frequencydependent finite-difference time-domain formulation for transient propagation in plasma.IEEE Trans.Antennas Propagat，39（1）：29-34.

Moczo P.1989.Finite-difference technique for SH-waves in 2-D media using irregular grids-application to the seismic response problem.Geophysical Journal International，99（2）：321-329.

Namiki T.1999.A new FDTD algorithm based on alternatingdirection implicit method.IEEE Transactions on Microwave Theory and Techniques，47（10）：2003-2007.

Shen B.1994.A study on wave eguation theory for ground——penetrating radar and forward mooelling.Computing Techniques for Geophysical and Geochemical Exploration（in Chinese），16（1）：29-33.

Teixeira F L.2008.Time-domain finite-difference and finite-element methods for Maxwell equations in complex media.IEEE Transactions on Antennas and Propagation，56（8）：2150-2166.

Tessmer E.2000.Seismic finite-difference modeling with spatially varying time steps.Geophysics，65（4）：1290-1293.

Tian G，Lin J X，Wang B B，etal.2011.Simulation and analysis reflections interference from above surface objects of ground penetrating radar.Chinese J.Geophys.（in Chinese），54（10）：2639-2651.

Wang B Z，Wang Y J，Yu W H，etal.2001.A hybrid 2-D ADIFDTD subgridding scheme for modeling on-chip interconnects.IEEE Transactions on Advanced Packaging，24（4）：528-533.

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

Zeng Z F，Liu S X，F(xiàn)eng X.2010.Theory and application of Ground Penetrating Radar（in Chinese）.Beijing：Electronic Industry Press.

Zhang H，Li Z C.2011.Forward modeling of carbonate fracture reservoir based on time-space dual variable grid algorithm.Journal of China University of Petroleum：Edition of Natural Science（in Chinese），35（3）：51-57.

Zheng K S.2005.Study of FDTD Parallel Computing and ADIFDTD Method（in Chinese）［D.Ph.thesis］.Xi′an：Institute of Science，Xidian University.

Zivanovic S S，Yee K S，Mei K K.1991.A subgridding method for the time-domain finite-difference method to solve Maxwell′s equations.IEEE Transactions on Microwave Theory and Techniques，39（3）：471-479.

附中文參考文獻

底青云，王妙月.1999.雷達波有限元仿真模擬.地球物理學報，42（6）：818-825.

方廣有，張忠治，汪文秉.1996.淺層地下目標時域電磁散射問題的數(shù)值模擬.電子與信息學報，18（3）：276-282.

方宏遠，林皋.2013.基于辛算法模擬探地雷達在復雜地電模型中的傳播.地球物理學報，56（2）：653-659.

馮德山，陳承申，戴前偉.2010.基于UPML邊界條件的交替方向隱式有限差分法GPR全波場數(shù)值模擬.地球物理學報，53（10）：2484-2496.

馮德山，陳承申，王洪華.2012.基于混合邊界條件的有限單元法GPR正演模擬.地球物理學報，55（11）：3774-3785.

馮德山，戴前偉，甕晶波.2007.時域多分辨率法在探地雷達三維正演模擬中的應用.中南大學學報：自然科學版，38（5）：975-980.

馮德山，王洪華，戴前偉.2013.基于無單元Galerkin法探地雷達正演模擬.地球物理學報，56（1）：298-308.

馮晅，鄒立龍，劉財?shù)?2011.全極化探地雷達正演模擬.地球物理學報，54（2）：349-357.

葛德彪，閆玉波.2005.電磁波時域有限差分方法（第二版）.西安：西安電子科技大學出版社.

黃超，董良國.2009a.可變網格與局部時間步長的高階差分地震波數(shù)值模擬.地球物理學報，52（1）：176-186.

黃超，董良國.2009b.可變網格與局部時間步長的交錯網格高階差分彈性波模擬.地球物理學報，52（11）：2870-2878.

李靜，曾昭發(fā)，吳豐收等.2010.探地雷達三維高階時域有限差分法模擬研究.地球物理學報，53（4）：974-981.

李慶偉.2008.有耗色散地質介質中的GPR脈沖傳播研究［博士論文］.成都：成都理工大學.

李展輝，黃清華，王彥賓.2009.三維錯格時域偽譜法在頻散介質井中雷達模擬中的應用.地球物理學報，52（7）：1915-1922.

劉少斌，劉崧，洪偉.2010.色散介質時域有限差分方法.北京：科學出版社.

劉四新，曾昭發(fā).2007.頻散介質中地質雷達波傳播的數(shù)值模擬.地球物理學報，50（1）：320-326.

沈飚.1994.探地雷達波波動方程研究及其正演模擬.物探化探計算技術，16（1）：29-33.

田鋼，林金鑫，王幫兵等.2011.探地雷達地面以上物體反射干擾特征模擬和分析.地球物理學報，54（10）：2639-2651.

曾昭發(fā)，劉四新，馮晅等.2010.探地雷達原理與應用.北京：電子工業(yè)出版社.

張慧，李振春.2011.基于時空雙變網格算法的碳酸鹽巖裂縫型儲層正演模擬.中國石油大學學報：自然科學版，35（3）：51-57.

鄭奎松.2005.FDTD網絡并行計算及ADI-FDTD方法研究［博士論文］.西安：西安電子科技大學理學院.