基于變換光學有限差分探地雷達數(shù)值模擬研究

李靜， 劉津杰， 曾昭發(fā)， 劉鳳山

1 吉林大學地球探測科學與技術(shù)學院， 長春　130026　2 美國特拉華州立大學應(yīng)用數(shù)學研究中心， 多佛　19901

?

基于變換光學有限差分探地雷達數(shù)值模擬研究

李靜1,2， 劉津杰2， 曾昭發(fā)1， 劉鳳山2

1 吉林大學地球探測科學與技術(shù)學院， 長春1300262 美國特拉華州立大學應(yīng)用數(shù)學研究中心， 多佛19901

摘要在采用有限差分方法開展探地雷達復雜目標體精細結(jié)構(gòu)模擬時，為了提高計算精度，常采用非均勻網(wǎng)格對目標區(qū)域劃分小尺寸的網(wǎng)格，以壓制離散網(wǎng)格頻散現(xiàn)象和保證有限差分方法的穩(wěn)定性.常規(guī)非均勻網(wǎng)格和自適應(yīng)亞網(wǎng)格技術(shù)在網(wǎng)格剖分數(shù)量和粗細網(wǎng)格邊界處理上難以達到計算效率和計算精度的均衡.本文根據(jù)隱形斗篷 (invisible cloak) 理論，將基于變換光學(Transformation optics)理論應(yīng)用于有限差分探地雷達數(shù)值計算中.該理論的主要思想是基于目標參數(shù)變化而保持電磁場的傳播不變性，在坐標變換后，Maxwell方程的形式可以維持不變，而使得相對介電常數(shù)與磁導率的表達式變得復雜.通過這種方式可以虛擬地擴大目標體所占的網(wǎng)格節(jié)點數(shù)，減少背景介質(zhì)區(qū)域的網(wǎng)格數(shù)，不增加模型空間的網(wǎng)格總數(shù).另外，這種網(wǎng)格劃分方式不但提高了計算效率，同時也可以克服亞網(wǎng)格技術(shù)邊界反射誤差的影響.本文推導實現(xiàn)了基于變換光學的二維有限差分方法，通過典型探地雷達模型測試，對比分析了該方法與常規(guī)有限差分、變網(wǎng)格有限差分和自適應(yīng)亞網(wǎng)格有限差分的優(yōu)缺點.計算結(jié)果驗證了基于變換光學的有限差分可用于探地雷達目標精細結(jié)構(gòu)模擬，具有較高的計算精度和計算效率.

關(guān)鍵詞探地雷達； 非均勻網(wǎng)格； 亞網(wǎng)格； 變換光學； 有限差分

1引言

有限差分方法是地球物理正演模擬中重要的數(shù)值計算方法之一，它以計算精度高、算法穩(wěn)定而被廣泛應(yīng)用于地震波(彈性波，聲波)、低頻電磁波和高頻探地雷達正演計算 (李靜等，2010).在開展探地雷達復雜目標體模擬時，為了提高數(shù)值計算精度，可以對整個模型區(qū)域采用較小的離散網(wǎng)格，以壓制網(wǎng)格頻散現(xiàn)象和保持有限差分方法的計算穩(wěn)定性.但模型網(wǎng)格剖分越細，占用計算內(nèi)存將增大，計算效率降低，同時長時間的迭代計算將引起誤差積累(馮德山等，2010；Li et al., 2012).采用非均勻網(wǎng)格劃分方法可以在保證計算效率的同時有效提高計算精度.Moczo (1989) 首次提出了有限差分可變網(wǎng)格思想，之后，研究者(趙海波等，2007；Diamanti等，2009； 黃超和董良國，2009a, 2009b；孫林潔等，2011； 張慧和李振春，2011；馮德山等，2014)先后將這種思想分別應(yīng)用到聲波、彈性波和電磁波等地球物理方法有限差分數(shù)值模擬計算中，打破了傳統(tǒng)固定網(wǎng)格方式.由介質(zhì)變化程度來決定離散網(wǎng)格的大小，在目標區(qū)域采用細網(wǎng)格剖分，在背景介質(zhì)區(qū)域采用較大空間網(wǎng)格剖分.這樣的方式不僅減少了計算量，又兼顧了方法穩(wěn)定性與計算精度.常用的非均勻網(wǎng)格包括固定方向變網(wǎng)格有限差分方法(Non-uniform FDTD)和粗細網(wǎng)格交錯自適應(yīng)局部亞網(wǎng)格有限差分方法 (Adaptive Mesh Refinement Sub-grid FDTD).上述方法存在以下問題：變網(wǎng)格有限差分方法采用固定方向的非均勻網(wǎng)格，即需要在一個方向整體進行細網(wǎng)格劃分，采用統(tǒng)一的時間步長，當粗細網(wǎng)格差距較大時需要用很小的時間步長才能滿足穩(wěn)定性要求，這樣在一定程度上浪費了計算資源，而且難以實現(xiàn)對任意區(qū)域的精細目標結(jié)構(gòu)的任意劃分(Wang and Liang, 2007, Pitarka, 1999).粗細網(wǎng)格交錯的亞網(wǎng)格技術(shù)雖然能自適應(yīng)對任意區(qū)域劃分不同的網(wǎng)格尺寸，并且采用不同的時間步長，但粗細網(wǎng)格交界的波場傳遞處理始終是制約該方法發(fā)展的瓶頸，特別是當網(wǎng)格尺寸比較大時，邊界場值迭代往往出現(xiàn)發(fā)散，影響計算精度(謝姣,2010； Xiao et al.,2007).

如何克服以上問題，是非均勻網(wǎng)格有限差分方法發(fā)展的關(guān)鍵技術(shù)之一.本文介紹了一種新的非均勻網(wǎng)格有限差分方法，該方法基于麥克斯韋方程的坐標變換不變性，通過坐標變換增加目標區(qū)域參數(shù)所占的網(wǎng)格數(shù)目，減少背景區(qū)域網(wǎng)格數(shù)據(jù)以保證整體網(wǎng)格數(shù)目不變，達到提高計算精度的目的，這種方法被稱為變換光學(transformation optics)有限差分方法.變換光學由Pendry等(2006) 在Science雜志上發(fā)表了一個符合 Maxwell 方程的隱形斗篷理論(Cloak Theory).該理論提出至今，已經(jīng)成了一個十分熱門的研究領(lǐng)域，可應(yīng)用于光線隱身設(shè)備設(shè)計、超透鏡以及聲波隱身設(shè)備等的設(shè)計和制造.Pendry等(1996)提出了基于物質(zhì)結(jié)構(gòu)參數(shù)變化而電磁場的傳播不變性，在任意坐標變換之下，Maxwell 方程式可以維持不變，而相對介電常數(shù)與磁導率的表達式變得更復雜.對于這個坐標變換可以采用以下兩種觀點做解釋(Chen et al, 2010)：第一種是空間電磁場和介質(zhì)都沒有變，只是坐標系改變 (笛卡兒坐標換成球坐標)，這是介質(zhì)相對于坐標系的改變；第二種是空間、介質(zhì)與電磁場都改變，但坐標系沒變.變換光學方法的設(shè)計關(guān)鍵，即在先采用第一種觀點做坐標變換，再將變換結(jié)果用第二種觀點解釋.變換之后的介質(zhì)與原來的介質(zhì)之間存在一一對應(yīng)關(guān)系.變換后的相對介電常數(shù)與磁導率具有跟目標本身一樣的電磁特性，只要連續(xù)地改變不同位置的目標參數(shù)，就可以建構(gòu)出任何可能的電磁介質(zhì)，實現(xiàn)坐標變換介質(zhì) (transformed medium) .這種方法仍然沒有明確定義而暫時被稱作變換光學(Transform Optics).

基于上述變換光學思想，本文推導實現(xiàn)了基于變換光學的二維有限差分方法并應(yīng)用于探地雷達目標數(shù)值計算.針對目標體的精細結(jié)構(gòu)模型，相比常規(guī)固定方向變網(wǎng)格、粗細網(wǎng)格交錯局部亞網(wǎng)格，變換光學有限差分方法不僅很好地提高了模型的計算精度，在計算效率方面也得到了很好的改善.

2固定方向變網(wǎng)格有限差分方法

固定方向變網(wǎng)格方法是沿固定方向?qū)δ繕藚^(qū)域劃分較細的網(wǎng)格，而在其他周圍背景介質(zhì)位置以及介質(zhì)變化緩慢的區(qū)域使用較粗網(wǎng)格.圖1是基于垂直方向劃分非均勻網(wǎng)格的示意圖，通過非均勻網(wǎng)格的方法通過減少數(shù)值色散，改善了精度.這種非均勻網(wǎng)格差分方法采用統(tǒng)一的時間迭代步長，需要滿足最小網(wǎng)格尺寸的穩(wěn)定性條件.在粗細網(wǎng)格內(nèi)部仍然采用標準的交錯網(wǎng)格迭代形式，迭代公式與常規(guī)二階差分格式相同 (王紅，2013).以y方向劃分非均勻網(wǎng)格為例，在x方向的空間步長分別為Δx和δx，在y方向的網(wǎng)格尺寸統(tǒng)一為Δy，根據(jù)完美磁導邊界的電磁場分布形式，只需要計算在網(wǎng)格交界面的電場分量Ez(Taflove, 1995)：

(1)

圖1　縱向非均勻網(wǎng)格劃分示意圖Fig.1　The map of non-uniform gird along y direction

3亞網(wǎng)格有限差分方法

圖2　二維亞網(wǎng)格分布示意圖Fig.2　The map of 2D sub-grid FDTD

亞網(wǎng)格技術(shù)是在目標體的位置區(qū)域采用細網(wǎng)格尺寸，對于周圍其他位置仍然采用粗網(wǎng)格尺寸，且在不同的網(wǎng)格區(qū)域采用不同的時間步長.通過對目標區(qū)域采用小網(wǎng)格尺寸能保證足夠的計算精度，另外又節(jié)約內(nèi)存和運行時間(Liu et al., 2009)，其基本原理見附錄.在粗網(wǎng)格和細網(wǎng)格區(qū)域內(nèi)部仍然采用普通的均勻網(wǎng)格差分迭代式進行遞推計算.需要注意的是，在

粗細網(wǎng)格交界位置，由于存在網(wǎng)格大小的不連續(xù)性，特別是當粗細網(wǎng)格尺寸比例較大時，不可避免地將帶來非物理性反射，影響亞網(wǎng)格技術(shù)穩(wěn)定性，也是制約亞網(wǎng)格技術(shù)發(fā)展的關(guān)鍵因素之一.在亞網(wǎng)格技術(shù)研究中，如何解決粗細網(wǎng)格交界處的反射誤差是一直以來的研究重點(Diamanti and Giannopoulos, 2009).

4基于變換光學的有限差分方法

基于變換光學的有限差分方法 (TO-FDTD) 應(yīng)用坐標變換來放大局部目標區(qū)域，使得目標區(qū)域包含更多的網(wǎng)格單元.亞網(wǎng)格或非均勻網(wǎng)格FDTD方法要求時間上逐步迭代，并且常常會遇到時間迭代后期不穩(wěn)定的問題.為了避免由變換光學引起的時間步長變小，將介質(zhì)參數(shù)矩陣映射到一種色散介質(zhì)模型中.時間域麥克斯韋方程組可表示為(Liu et al., 2014a)

(2)

將笛卡兒坐標系(x,y,z)進行坐標轉(zhuǎn)換得到另一個坐標系(x′,y′,z′)之后，得到新的麥克斯韋方程組為

(3)

(4)

圖3　基于變換光學有限差分方法變換原理示意圖(a) 傳統(tǒng)的FDTD; (b) 變換光學的FDTD.Fig.3　The illustration of enlarging a small region by using transformation optics(a) Traditional FDTD; (b) Transformation optics FDTD.

三次坐標轉(zhuǎn)換的雅克比矩陣分別為

(1) 由笛卡兒網(wǎng)格(x, y, z)轉(zhuǎn)換到圓柱坐標系(r,θ,z):

(5)

其雅克比矩陣為

(6)

(2) 由(r,θ,z)到(r′,θ′,z)的坐標轉(zhuǎn)換:

(7)

其雅克比矩陣為

(8)

r=f(r′)

(9)

(3) 由圓柱坐標系(r′,θ′,z)轉(zhuǎn)換為笛卡兒網(wǎng)格(x′,y′,z), 有

(10)

雅克比矩陣為

(11)

結(jié)合以上三次坐標轉(zhuǎn)換，得到由(x,y,z)到(x′,y′,z′)的雅克比矩陣：

(12)

變換光學有限差分實現(xiàn)的基本思想可以由如圖4所示的圓形體模型解釋.圖4a是在常規(guī)FDTD相對介電常數(shù)值的分布.模型網(wǎng)格數(shù)據(jù)為140×140，目標體和背景介質(zhì)相對常數(shù)分別為81和6.圓形目標體大約占10個網(wǎng)格.采用變換光學有限差分模式下，圓形體目標相對介電常數(shù)變?yōu)榧s20，體積大小約占20個網(wǎng)格，其代價是周圍背景介質(zhì)是非均勻的漸變參數(shù)分布.為了更為直觀地分析參數(shù)分布特點，垂直抽取圓心所在位置的單道相對介電常數(shù)值對比，如圖4(c、d)所示.圖中標準FDTD數(shù)值分布均勻，目標與背景介質(zhì)區(qū)分明顯，含水圓形體相對介電常數(shù)為81，背景介質(zhì)及空氣層分別為6和1.TO-FDTD方法中，在坐標變換之后，環(huán)狀區(qū)域的相對介電常數(shù)變成半徑的函數(shù)(變換前是常數(shù)).按照上述坐標變換過程使坐標拉伸壓縮.在坐標拉伸的區(qū)域，相對介電常數(shù)變小.在坐標壓縮的區(qū)域，相對介電常數(shù)變大.這就是為什么中間地帶(拉伸)相對介電常數(shù)變小，外環(huán)區(qū)域(壓縮)變大的緣故.

圖4　標準FDTD和變換光學FDTD模型相對介電常數(shù)分布對比(a) 標準FDTD二維矩陣； (b) TO-FDTD二維矩陣； (c) 標準FDTD垂直單道網(wǎng)格； (d) TO-FDTD垂直單道網(wǎng)格.Fig.4　Model dielectric constant parameter comparison of standard FDTD and TO-FDTD(a) The 2D model of standard FDTD； (b) The 2D model of TO FDTD； (c) The single trace of standard FDTD； (d) The single trace of TO FDTD.

5模型測試

圖5　不同網(wǎng)格劃分方法測試模型Fig.5　The 2D test model with different FDTD method

針對常規(guī)非均勻網(wǎng)格、自適應(yīng)亞網(wǎng)格以及基于變換光學的有限差分方法設(shè)計不同的介質(zhì)模型開展探地雷達數(shù)值模擬計算，比較不同方法在計算精度和效率方面的應(yīng)用效果.圖5是兩層介質(zhì)模型，上層介質(zhì)包含兩個大小不同的隨機介質(zhì)圓形體.隨機目標體采用混合型自相關(guān)函數(shù)生成，之所以采用隨機介質(zhì)目標體是為了增加模型復雜程度以測試不同方法的優(yōu)缺點.模型大小為1 m×1 m.最小網(wǎng)格尺寸為0.01 m，上下層介質(zhì)相對介電常數(shù)和電導率分別為12，6和0.01 S·m-1，0.005 S·m-1.其中隨機介質(zhì)目標體位于上層介質(zhì)，圓形目標體半徑分別為0.03 m和0.01 m，隨機圓形目標相對介電常數(shù)變化范圍為30～35，右側(cè)均勻球體為30.采用天線中心頻率為900 MHz的發(fā)射源，收發(fā)天線位于地表分界面.圖6為采用常規(guī)有限差分、漸變非均勻網(wǎng)格和自適應(yīng)亞網(wǎng)格技術(shù)的有限差分數(shù)值模擬結(jié)果.三種方法均能夠?qū)﹄S機圓形體目標及其層界面得到很好地模擬結(jié)果，常規(guī)有限差分結(jié)果中圓形體目標上下界面反射較為模糊，在圓形體反射信號以下存在多次波干擾，非均勻網(wǎng)格有限差分可以很好地壓制干擾波的存在以及網(wǎng)格數(shù)值誤差，模擬結(jié)果有更高的信噪比，目標信號清晰，雜波干擾較少.

但是，圖6b非均勻漸變網(wǎng)格結(jié)果在4 ns和12 ns所示的位置有固定直線反射信號，這是因為采用了橫向網(wǎng)格加密方法，在漸變交界區(qū)域存在異常反射，網(wǎng)格尺寸突變造成.表1是三種方法在相同網(wǎng)格尺寸條件下的計算效率對比，可以看出，基于自適應(yīng)亞網(wǎng)格的計算時間和常規(guī)方法相當，而漸變非均勻網(wǎng)格在相同條件下，計算時間是常規(guī)方法的15倍左右.綜合考慮計算效率和計算精度，漸變非均勻網(wǎng)格計算有不錯的計算精度，但是計算效率較低.自適應(yīng)變換亞網(wǎng)格計算結(jié)果在保證計算精度的同時也能較好地保證計算效率.另外，亞網(wǎng)格技術(shù)可以根據(jù)模型位置任意劃分網(wǎng)格加密區(qū)域，而漸變非均勻網(wǎng)格只能在x或y方向整體加密網(wǎng)格，計算效率較低且對模型劃分沒有自適應(yīng)亞網(wǎng)格方法靈活.

表1　不同有限差分方法計算時間對比

上述結(jié)果驗證了常規(guī)非均勻網(wǎng)格在計算效率和計算精度上存在的不足.圖7是常規(guī)有限差分、自適應(yīng)亞網(wǎng)格和基于變換光學的有限差分對比.圖7所示為兩層介質(zhì)模型，模型大小為1 m×1 m.最小網(wǎng)格尺寸為0.01 m，上下層介質(zhì)相對介電常數(shù)和電導率分別為6，8和0.01 S·m-1，0.005 S·m-1.其中隨機介質(zhì)目標體位于上層介質(zhì).在上層模型中設(shè)計了微小隨機介質(zhì)圓形體模型，圓形體半徑僅為0.003 m，小于網(wǎng)格劃分的最小尺寸(0.01 m)，如圖7中放大區(qū)域所示.采用天線中心頻率為900 MHz的發(fā)射源，收發(fā)天線位于地表分界面.

圖8所示為三種方法的計算測試結(jié)果，圖8a為常規(guī)有限差分方法模擬結(jié)果，由于目標尺寸小于最小網(wǎng)格，無法獲得任何目標信息.在相同的條件下，采用局部亞網(wǎng)格有限差分，對目標區(qū)域采用Δx/4網(wǎng)格尺寸，圖8b中能看到明顯的目標信號，箭頭所示位置是因為局部粗細網(wǎng)格界面處理造成的計算誤差，圖8c所示為采用基于變換光學的有限差分計算結(jié)果，不僅能夠?qū)δＰ陀泻芎玫姆从常捎诓恍枰獙W(wǎng)格尺寸劃分，即不存在亞網(wǎng)格技術(shù)的邊界誤差，顯得尤為“干凈”.

圖6　隨機介質(zhì)圓形體模型不同方法測試結(jié)果(a) 常規(guī)有限差分方法； (b) 非均勻漸變網(wǎng)格有限差分方法； (c) 基于自適應(yīng)亞網(wǎng)格有限差分方法.Fig.6　The synthetic result of two random circle target model with different method(a) Standard FDTD; (b) Non-uniform FDTD; (c) AMR FDTD.

圖7　局部微小隨機介質(zhì)模型示意圖(虛框所示為局部微小目標放大)Fig.7　The layer model with single local small random target (right figure is the zoom of the dotted box)

圖8　不同有限差分方法測試結(jié)果(a) 常規(guī)有限差分方法； (b) 基于自適應(yīng)亞網(wǎng)格有限差分方法； (c) 基于變換光學的有限差分方法.Fig.8　The synthetic result of local small random target with different method(a) Standard FDTD; (b) AMR FDTD; (c) TO-FDTD.

圖9　局部微小隨機介質(zhì)模型示意圖(右側(cè)小圖為黑框局部放大)Fig.9　The layer model with two local small random target (the right figure is the zoom of the dotted box)

圖9所示的兩個半徑不同的圓形目標體，半徑尺寸分別為0.003 m和0.01 m，其他參數(shù)設(shè)置與上述模型相同.采用相同的計算方式可以得到如圖10所示的計算結(jié)果，對比三種方法的結(jié)果可以看出，基于變換光學的有限差分技術(shù)通過對目標區(qū)域進行變換加密，不僅能夠?qū)ξ⑿∧繕擞星逦淖R別，而且信噪比較高，不存在邊界誤差干擾.層界面反射由于采用了正常的網(wǎng)格模式，存在的繞射干擾和多次波與常規(guī)方法和亞網(wǎng)格技術(shù)幾乎相近.

圖10　不同有限差分方法測試結(jié)果Fig.10　The synthetic result of local small random target with different method

圖11為充水環(huán)狀模型，通過該模型對比測試不同網(wǎng)格尺寸條件下變換光學和常規(guī)有限差分在計算精度和計算效率上的差異.在兩層介質(zhì)交界處存在一個充水環(huán)狀圓形模型，圓形內(nèi)半徑為0.03 m，外半徑為0.04 m.采用上述相同的計算參數(shù)進行數(shù)值模擬，圖12是選取收發(fā)天線位于地表中心位置的單道雷達波信號，圖中分別包括最小網(wǎng)格尺寸為0.0025 m，0.005 m和0.01 m的常規(guī)有限差分計算結(jié)果和最小網(wǎng)格尺寸為0.01 m和0.015 m的變換光學有限差分計算結(jié)果.對比可以看出，常規(guī)方法在網(wǎng)格尺寸為0.005 m以及0.0025 m條件下計算結(jié)果與變換光學有限差分相符，而最小網(wǎng)格為0.01 m的結(jié)果存在明顯的計算誤差.表2是不同網(wǎng)格尺寸的計算效率對比.相對于常規(guī)有限差分方法，基于變換光學有限差分方法可以在較大網(wǎng)格尺寸條件下獲得與常規(guī)方法小網(wǎng)格相當?shù)慕Y(jié)果，計算效率方面卻有較大的提高.例如常規(guī)方法網(wǎng)格為0.005 m需要52.3 s的計算時間，而變換光學有限差分網(wǎng)格為0.01 m時只需要35.3 s，在計算大規(guī)模長時間模型時，計算效率的優(yōu)勢將更加明顯.

圖11　充水環(huán)狀目標模型示意圖Fig.11　The layer model with filling water ring

圖12　采用基于變換光學有限差分和常規(guī)有限差分單道信號對比.(a) 完整單道信號； (b) 局部目標信號放大顯示Fig.12　The comparison of single GPR signal with standard FDTD and TO-FDTD method (the right figure is the zoom of the dotted box)

表2　不同F(xiàn)DTD網(wǎng)格尺寸計算效率對比

基于變換光學的有限差分方法在探地雷達數(shù)值模擬中不僅能克服亞網(wǎng)格差分存在邊界處理誤差的影響，而且可以改善常規(guī)非均勻網(wǎng)格計算量大的問題，保持計算精度和計算效率的平衡.對于精細目標結(jié)構(gòu)探地雷達模擬具有適用性，例如在探地雷達天線模擬中，饋電點的距離和天線尺寸往往存在較大的比例，特別是低頻天線，饋電點的距離和天線尺寸比例接近1∶100(以饋電點1 cm，天線尺寸1 m為例)，常規(guī)模擬方法為了模擬饋電點的影響，需要小于1 cm的網(wǎng)格尺寸，而采用本文方法，在滿足穩(wěn)定性條件的前提下，最小網(wǎng)格尺寸可選擇大于1 cm.另外，變換光學也可以擴展到聲波及彈性波差分方法模擬，在地震勘探大尺度目標的局部精細結(jié)構(gòu)，隱形斗篷理論用于地震災(zāi)害防御 (Nicoletti, 2014)是當前有限差分研究的熱點問題之一.

6結(jié)論

本文推導了基于變換光學的非均勻網(wǎng)格二維有限差分算法并應(yīng)用于探地雷達正演模擬計算，對比分析了本文方法與常規(guī)FDTD方法，常規(guī)非均勻網(wǎng)格以及自適應(yīng)亞網(wǎng)格(AMR)在計算效率以及計算精度上的差異.正演模擬結(jié)果表明：對于復雜隨機介質(zhì)模型，非均勻網(wǎng)格方法相對于常規(guī)有限差分可以很好地提高計算精度但是計算量也大幅度增加；自適應(yīng)亞網(wǎng)格(AMR)有限差分方法在提高計算精度的同時可以很好地保證計算效率，但是自適應(yīng)網(wǎng)格方法在處理粗細網(wǎng)格交界面時，很難避免邊界反射的影響，這也是制約亞網(wǎng)格技術(shù)發(fā)展的瓶頸.本文提出的基于隱形斗篷理論 (Cloak) 的有限差分方法，通過對目標的屬性參數(shù)的坐標變換，增加目標所占網(wǎng)格的數(shù)目，減少周圍介質(zhì)所占網(wǎng)格數(shù)目，其本質(zhì)上并沒有對網(wǎng)格進行加密，即沒有增加計算量，也不存在邊界反射問題.通過與常規(guī)有限差分方法對比，基于變換光學的有限差分方法在保證計算效率的同時，很好地消除了邊界效應(yīng)的影響，具有很好的計算精度，為探地雷達以及其他地球物理勘探方法的復雜介質(zhì)精細目標結(jié)構(gòu)數(shù)值模擬和參數(shù)反演提供了很好的數(shù)值計算手段.

附錄A亞網(wǎng)格非均勻網(wǎng)格原理

(A1)

對于二維情況，可以寫成

(A2)

(A3)

可以得到

(A4)

式(A4)應(yīng)用于粗細網(wǎng)格邊界.對于亞網(wǎng)格，從當前時刻t=nΔtc算起，需要沿時間軸往前推進nf步.在這里設(shè)nf=4，則需要分別為n+Δtf，n+2Δtf，n+3Δtf，n+4Δtf時刻，其中每一步Ezf的計算都需要用到上式.D可表示為

(A5)

(A6)

References

Chen H Y, Chan C T, Sheng P. 2010. Transformation optics and metamaterials.Nat.Mater., 9(5): 387-396.

Diamanti N, Giannopoulos A. 2009. Implementation of ADI-FDTD subgrids in ground penetrating radar FDTD models.JournalofAppliedGeophysics, 67(4): 309-317.

Feng D S, Chen C S, Dai Q W. 2010. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.ChineseJournalofGeophysics(in Chinese), 53(10): 2484-2496, doi: 10.3969/j.issn.0001-5733.2010.10.022.

Feng D S, Chen J W, Wu Q. 2014. A hybrid ADI-FDTD subgridding scheme for efficient GPR simulation of dispersion media.ChineseJournalofGeophysics(in Chinese), 57(4): 1322-1334, doi: 10.6038/cjg20140429.

Huang C, Dong L G. 2009a. High-order finite-difference method in seismic wave simulation with variable grids and local time-steps.ChineseJournalofGeophysics(in Chinese), 52(1): 176-186.

Huang C, Dong L G. 2009b. Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps.ChineseJournalofGeophysics(in Chinese), 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.

Li J, Zeng Z F, Huang L, et al. 2012. GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD.Comp.Geosci., 49: 121-130.

Li J, Zeng Z F, Wu F S, et al. 2010. Study of three dimension high-order FDTD simulation for GPR.ChineseJournalofGeophysics(in Chinese), 53(4): 974-981, doi: 10.3969/j.issn.0001-5733.2010.04.022.

Liu J J, Brio M, Moloney J V. 2014a. Transformation optics based local mesh refinement for solving Maxwell′s equations.JournalofComputationalPhysics, 258: 359-370.

Liu J J, Brio M, Moloney J V. 2009. Overlapping Yee FDTD method on nonorthogonal grids.J.Sci.Comput., 39(1):129-143.

Moczo P. 1989. Finite difference technique for SH-waves in 2-D media using irregular grids-application to the seismic response problem.Geophys.J.Int., 99(2): 321-329.

Nicoletti O. 2014, Seismic cloaks.NatureMaterials, 13: 428 doi:10.1038/nmat3968Pendry J B, Schurig D, Smith D R. 2006. Controlling electromagnetic fields.Science, 312(5781): 1780-1782.

Pitarka A. 1999. 3D finite-difference modeling of seismic motion using staggered grids with nonuniform spacing.Bull.seis.Soc.Am., 89(1), 54-68.

Sun C Y, Ding Y C. 2012. The precision analysis of wave equation finite difference double variable grid algorithm.OilGeophysicalProspecting(in Chinese), 47(4): 545-551. Sun L J, Yin X Y. 2011. A finite-difference scheme based on PML

boundary condition with high power grid step variation.ChineseJournalofGeophysics(in Chinese), 54(6): 1614-1623, doi: 10.3969/ j.issn.0001-5733. 2011.06.021. Taflove A. 1995. Computational Electrodynamics: the Finite Difference Time Domain Method. USA: Artech House.Xiao K, Pommerenke D J, Drewniak J L. 2007. A three-dimensional FDTD subgridding algorithm with separated temporal and spatial interfaces and related stability analysis.IEEETransactionsonAntennasandPropagation, 55(7): 1981-1990.

Wang H. 2013. Numerical modeling of acoustic wave equation with spatial varying meshes[Master thesis] (in Chinese). Changsha: Central South University.Wang L N, Liang C H. 2006. A new implementation of CFS-PML for ADI-FDTD method.MicrowaveandOpticalTechnologyLetters, 48(10): 1924-1928.

Xie J. 2010. Research on Subgrid Technology for FDTD[Master thesis](in Chinese). Nanjing: Nanjing University of Aeronautics and Astronautics.Zhang H, Li Z C. 2011. Seismic wave simulation method based on dual-variable grid.ChineseJournalofGeophysics(in Chinese), 54(1): 77-86, doi: 10.3969/j.issn.0001-5733.2011.01.009.

附中文參考文獻

馮德山, 陳承申, 戴前偉. 2010. 基于UPML邊界條件的交替方向隱式有限差分法GPR全波場數(shù)值模擬. 地球物理學報, 53(10): 2484-2496, doi: 10.3969/j.issn.0001-5733.2010.10.022. 馮德山, 陳佳維, 吳奇. 2014. 混合ADI-FDTD亞網(wǎng)格技術(shù)在探地雷達頻散媒質(zhì)中的高效正演. 地球物理學報, 57(4): 1322-1334, doi: 10.6038/cjg20140429.

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

黃超, 董良國. 2009b. 可變網(wǎng)格與局部時間步長的交錯網(wǎng)格高階差分彈性波模擬. 地球物理學報, 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.

李靜, 曾昭發(fā), 吳豐收等. 2010. 探地雷達三維高階時域有限差分法模擬研究. 地球物理學報, 53(4): 974-981, doi: 10.3969/j.issn.0001-5733.2010.04.022.

孫林潔, 印興耀. 2011. 基于PML邊界條件的高倍可變網(wǎng)格有限差分數(shù)值模擬方法. 地球物理學報, 54(6): 1614-1623, doi: 10.3969/j.issn.0001-5733.2011.06.021.

王紅. 2013. 變網(wǎng)格步長聲波方程有限差分數(shù)值模擬[碩士學位論文]. 長沙: 中南大學.

謝姣. 2010. 時域有限差分方法亞網(wǎng)格技術(shù)的研究[碩士學位論文]. 南京: 南京航空航天大學.

張慧, 李振春. 2011. 基于雙變網(wǎng)格算法的地震波正演模擬. 地球物理學報, 54(1): 77-86, doi: 10.3969/j.issn.0001-5733.2011.01.009.

趙海波，王秀明．2007. 一種優(yōu)化的交錯變網(wǎng)格有限差分法及其在井間聲波中的應(yīng)用. 科學通報, 52(12): 1387-1395.

(本文編輯胡素芳)

基金項目國家自然科學青年基金項目(41504083)和中國博士后基金面上項目(2015M571366)聯(lián)合資助.

作者簡介李靜，男，1988年生，博士，吉林大學地球探測科學與技術(shù)學院講師，主要從事探地雷達數(shù)值模擬及電磁法正反演研究. E-mail：inter_lijing@jlu.edu.cn

doi:10.6038/cjg20160630 中圖分類號P631

收稿日期2015-09-28，2016-08-10收修定稿

Study of GPR simulation based on the transformation optics FDTD

LI Jing1,2, LIU Jin-Jie2, ZENG Zhao-Fa1, LIU Feng-Shan2

1CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China2AppliedMathematicsResearchCenter,DelawareStateUniversity,Dover,DE19901,USA

AbstractIn complex target fine structure GPR finite difference time domain (FDTD) simulation, in order to improve the computational accuracy, we commonly use the non-uniform grid to discretize small size grid in the target area to reduce the grid dispersion and keep stability. The conventional non-uniform grid and adaptive mesh refinement (AMR) sub-grid FDTD can not achieve balance between computational time and accuracy. In this paper, according to the invisible cloak theory, we apply the FDTD based on transformation optics (TO) in GPR complex target simulation. The basic theory of TO method is that the medium parameters change and the resulted electromagnetic field is invariant, and that the Maxwell equation can remain unchanged in coordinate transformation and make the expression of dielectric constant and magnetic permeability become more complex. It will increase the number of grid in target area and reduce that of background medium without changing the total grid. This grid mode can not only overcome the reflection error from sub-grid boundary, improve the target computational accuracy, but also save the computational time. This paper derives the basic formula of 2D FDFD with TO theory and compares the advantages and disadvantages of standard FDTD, Non-uniform FDTD, sub-grid FDTD and TO-FDTD with synthetic model. The results demonstrate that the TO-FDTD method can simulate the small target GPR signal response with high computational accuracy and efficiency.KeywordsGround penetrating radar (GPR)； Non-uniform grid； Sub-grid； Transformation optics； FDTD

李靜， 劉津杰， 曾昭發(fā)等. 2016. 基于變換光學有限差分探地雷達數(shù)值模擬研究.地球物理學報，59(6):2280-2289,doi:10.6038/cjg20160630.

Li J, Liu J J, Zeng Z F, et al. 2016. Study of GPR simulation based on the transformation optics FDTD.ChineseJ.Geophys. (in Chinese),59(6):2280-2289,doi:10.6038/cjg20160630.