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

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

?

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

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

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

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

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

1引言

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

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

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

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

固定方向變網(wǎng)格方法是沿固定方向?qū)δ繕?biāo)區(qū)域劃分較細(xì)的網(wǎng)格，而在其他周圍背景介質(zhì)位置以及介質(zhì)變化緩慢的區(qū)域使用較粗網(wǎng)格.圖1是基于垂直方向劃分非均勻網(wǎng)格的示意圖，通過非均勻網(wǎng)格的方法通過減少數(shù)值色散，改善了精度.這種非均勻網(wǎng)格差分方法采用統(tǒng)一的時(shí)間迭代步長(zhǎng)，需要滿足最小網(wǎng)格尺寸的穩(wěn)定性條件.在粗細(xì)網(wǎng)格內(nèi)部仍然采用標(biāo)準(zhǔn)的交錯(cuò)網(wǎng)格迭代形式，迭代公式與常規(guī)二階差分格式相同 (王紅，2013).以y方向劃分非均勻網(wǎng)格為例，在x方向的空間步長(zhǎng)分別為Δx和δx，在y方向的網(wǎng)格尺寸統(tǒng)一為Δy，根據(jù)完美磁導(dǎo)邊界的電磁場(chǎng)分布形式，只需要計(jì)算在網(wǎng)格交界面的電場(chǎ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ù)是在目標(biāo)體的位置區(qū)域采用細(xì)網(wǎng)格尺寸，對(duì)于周圍其他位置仍然采用粗網(wǎng)格尺寸，且在不同的網(wǎng)格區(qū)域采用不同的時(shí)間步長(zhǎng).通過對(duì)目標(biāo)區(qū)域采用小網(wǎng)格尺寸能保證足夠的計(jì)算精度，另外又節(jié)約內(nèi)存和運(yùn)行時(shí)間(Liu et al., 2009)，其基本原理見附錄.在粗網(wǎng)格和細(xì)網(wǎng)格區(qū)域內(nèi)部仍然采用普通的均勻網(wǎng)格差分迭代式進(jìn)行遞推計(jì)算.需要注意的是，在

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

4基于變換光學(xué)的有限差分方法

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

(2)

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

(3)

(4)

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

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

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

(5)

其雅克比矩陣為

(6)

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

(7)

其雅克比矩陣為

(8)

r=f(r′)

(9)

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

(10)

雅克比矩陣為

(11)

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

(12)

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

圖4　標(biāo)準(zhǔn)FDTD和變換光學(xué)FDTD模型相對(duì)介電常數(shù)分布對(duì)比(a) 標(biāo)準(zhǔn)FDTD二維矩陣； (b) TO-FDTD二維矩陣； (c) 標(biāo)準(zhǔn)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模型測(cè)試

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

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

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

表1　不同有限差分方法計(jì)算時(shí)間對(duì)比

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

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

圖6　隨機(jī)介質(zhì)圓形體模型不同方法測(cè)試結(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　局部微小隨機(jī)介質(zhì)模型示意圖(虛框所示為局部微小目標(biāo)放大)Fig.7　The layer model with single local small random target (right figure is the zoom of the dotted box)

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

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

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

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

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

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

圖12　采用基于變換光學(xué)有限差分和常規(guī)有限差分單道信號(hào)對(duì)比.(a) 完整單道信號(hào)； (b) 局部目標(biāo)信號(hào)放大顯示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)格尺寸計(jì)算效率對(duì)比

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

6結(jié)論

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

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

(A1)

對(duì)于二維情況，可以寫成

(A2)

(A3)

可以得到

(A4)

式(A4)應(yīng)用于粗細(xì)網(wǎng)格邊界.對(duì)于亞網(wǎng)格，從當(dāng)前時(shí)刻t=nΔtc算起，需要沿時(shí)間軸往前推進(jìn)nf步.在這里設(shè)nf=4，則需要分別為n+Δtf，n+2Δtf，n+3Δtf，n+4Δtf時(shí)刻，其中每一步Ezf的計(jì)算都需要用到上式.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.

附中文參考文獻(xiàn)

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

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

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

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

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

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

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

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

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

(本文編輯胡素芳)

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

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

doi:10.6038/cjg20160630 中圖分類號(hào)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. 基于變換光學(xué)有限差分探地雷達(dá)數(shù)值模擬研究.地球物理學(xué)報(bào)，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.