亞波長間距理想導體球陣列近區時間反演電磁場的快速求解?

龔志雙王秉中王任

(電子科技大學應用物理研究所,成都 610054)

(2017年11月22日收到;2017年12月27日收到修改稿)

1 引 言

時間反演技術能夠簡單地重構出輻射源的位置信息.在自由空間中,點源的反演場幅度在空間呈類似sinc函數的分布,其第一個過零點即對應半波長的位置,也就意味著極限聚焦分辨率為半波長[1?3].分辨率受限的本質原因是近場凋落波無法傳播至遠場,因此要從本質上突破半波長的極限,必須在近場區域對凋落波進行操作[4?7].在近場區域有散射體存在的情況下,時間反演技術能夠實現超分辨率(小于半波長)的聚焦[8?11].但是,前述研究主要通過仿真軟件進行反演場的計算,求解過程費時而且沒有嚴格的理論指導.因此,研究有散射體加載情況下反演場的解析解,就顯得尤為重要.金屬小球作為結構最簡單的散射體,也可以認為是復雜散射體的基本組成部分,是本文的主要研究對象.時間反演場的嚴格解析解可以通過時間反演腔理論進行求解,但前提是已知正向場對應的格林函數[12].

關于金屬球的散射格林函數早已有嚴格的解析解[13?17],但是,前述論文中給出的結果均是將不同區域的場分別展開為球矢量波函數疊加的形式,之后再利用邊界面上的場應該滿足的邊界條件對展開系數分別求解.如此求解出來的結果為無窮級數求和的形式,當空間中包含多個散射體時,上述方法的求解會非常復雜[18,19].用基于等效偶極子分析的思路可以較快速地求解復雜物體的散射格林函數[20,21],尤其在處理多散射體問題時優勢比較明顯.

本文在已有研究的基礎上,提出了球面波入射情況下基于等效偶極子分析的反演場快速求解方法.通過將金屬小球的散射場等效為電磁偶極子輻射場的疊加,極大地簡化格林函數的求解.在包含有多個散射體時,只需要將每一個金屬球的散射場均等效為相應的電、磁偶極子輻射場的疊加,在求解出等效偶極子強度之后,空間中總場對應的格林函數即可表示為所有等效電磁偶極子的貢獻之和.

本文第2部分詳細給出球面波入射情況下金屬小球的散射場求解及其相應的等效方法;第3部分給出總場對應并矢格林函數的表達式;第4部分給出幾種不同情況下相應時間反演場的解析求解結果以及商業軟件仿真結果,以驗證所提出的等效方法的正確性;最后,利用所提出的分析方法對加載散射體情況下源的聚焦分辨率做了初步研究.結果表明,利用所提等效模型可以對多金屬小球加載情況下的時間反演場進行快速求解,為實現時間反演超分辨特性的加載散射體結構設計提供理論指導.

2 偶極子源激勵情況散射場求解

2.1 水平電偶極子源激勵情況

求解模型示意圖見圖1.散射體為一位于坐標原點處、半徑為a的金屬球,激勵源為一處于r′=(0,0,b)位置處的水平電偶極子,并有b>a(即源在散射體外),相應的電偶極子矢量為pe=Ielex,其中Ie為電流強度,l為元的長度.場觀測點的球坐標位置為(r,θ,?).真空中介質的介電常數和磁導率分別為ε和μ.下面表達式中上標為“inc”的表示入射場,上標為“sca”的表示散射場,下標為“equ”的代表等效結果.k代表自由空間中傳播波數;ω為角頻率;η0為自由空間波阻抗;r為場觀測點到原點的距離;R為場觀測點到初始激勵源的距離.為第一類連帶勒讓德函數;分別為三種不同類型的球貝塞爾函數.

為了對凋落波進行操作使其能轉化為傳輸波,需要所加載的散射體尺寸與相應的空間頻率相對應,而凋落波對應較高的空間頻率分量,因此要求所加載的散射體尺寸相對自由空間波長較小.所以本文主要考慮加載散射體尺寸遠小于自由空間波長的情況(即ka<1).同時,因為凋落波隨傳輸距離呈指數衰減,需要在其衰落至太小之前對其進行操作,也即意味著需要在源的近場區進行加載.故本文只討論源與散射體之間間距小于半波長(即kb<π)的情形.

根據文獻[11]中的結果可以得到圖1所示場景散射場的表達式為

圖1 水平電偶極子激勵情況示意圖Fig.1.Diagram with an excitation source of electric dipole.

其中cen和den分別為

設場觀測面為z=b平面,將(2),(3)式代入(1)式,計算觀測面上其前兩項絕對值的最大值之比ρ隨a和b的變化情況,如圖2(a)所示.由圖2(a)可知,ρ=10的等高線近似為一條直線,對其進行線性擬合的結果為b=1.2646a+0.1168λ.另外,由貝塞爾函數的特性可知,更高次項的絕對值會越來越小.以ρ值大于10為依據可知,當b>1.2646a+0.1168λ且a<0.08λ時可以近似采用(1)式的第一項代表整體散射場.當a=0.03λ,b=0.3λ時,(1)式前三項絕對值大小沿觀測面上x軸的分布如圖2(b)所示,可以看到此時第二、三項的值確實遠小于第一項.

圖2 水平電流源激勵情況下散射電場級數展開式 (a)前兩項絕對值最大值之比隨a和b的變化情況;(b)前三項絕對值沿x軸分布情況(其中a=0.03λ,b=0.3λ,觀測面為z=b平面.以下所有類似場幅度值分布圖所對應的計算場景與此圖相同)Fig.2.Series expansion of the scattered electric field with horizontal excitation electric dipole:(a)Relationship of the ratio between the maximum absolute value of the first two terms with a and b;(b)absolute value distribution of the first three terms along x-axis in the plane z=b(with a=0.03λ,b=0.3λ.The calculation scenery bellow is the same as this figure).

將(2)式和(3)式代入(1)式得到第一項的表達式為

對比沿x方向電偶極子及沿y方向磁偶極子的輻射場[20]可以發現,此散射場等于一電偶極子和一磁偶極子輻射場的疊加.相應的電磁偶極子矢量分別為

2.2 垂直電偶極子源激勵情況

對于圖1所示的模型,當激勵源的極化方向為沿z軸時,電偶極子強度可表示為pe=Ielez.根據對稱性并且考慮小金屬球情況可得散射場表達式為

(6)式前三項絕對值分布的計算結果如圖3所示.仍以ρ>10為判斷標準,當b>1.0429a+0.1416λ

圖3 垂直電流源激勵情況下散射電場級數展開式 (a)前兩項絕對值最大值之比隨a和b的變化情況;(b)前三項絕對值沿x軸分布情況Fig.3.Series expansion of scattered electric field with vertical excitation electric dipole:(a)Relationship of the ratio of the maximum absolute value of the first two terms with a and b;(b)absolute value distribution of the first three terms along x-axis in the plane z=b.

且a<0.08λ時,可以只取(6)式的第一項代表整體散射場,其表達式為

對比沿z方向極化電偶極子的輻射場可知,此時散射場等于一電偶極子的輻射場.相應的等效電偶極子矢量為

2.3 水平磁偶極子源激勵情況

當激勵源為磁偶極子并且極化方向為沿x軸時,磁偶極子矢量可表示為pm=Imlex,其中Im為磁流強度.利用前面計算得到的電偶極子激勵情況的結果并根據電與磁的對偶性可得此時散射場的表達式為:

(10)式前三項絕對值分布計算結果如圖4所示.類似可得當b>0.9938a+0.1248λ且a<0.08λ時,其第一項的值遠大于高階項的值.取(10)式的第一項可得散射電場的表達式為:

對比沿x方向極化磁偶極子及沿y方向極化電偶極子的輻射場可以發現,此時散射場等于一磁偶極子與一電偶極子輻射場之和,相應的偶極子矢量分別為:

圖4 水平磁流源激勵情況下散射電場級數展開式 (a)前兩項絕對值最大值之比隨a和b的變化情況;(b)前三項絕對值沿x軸分布情況Fig.4.Series expansion of scattered electric fi eld with horizontal excitation magnetic dipole:(a)Relationship of ratio between maximum absolute value of the fi rst two terms with a and b;(b)absolute value distribution of the fi rst three terms along x-axis in the plane z=b.

2.4 垂直磁偶極子源激勵情況

當激勵源為一沿z向極化的磁偶極子時,磁偶極子強度矢量為pm=Imlez.利用對稱性可以由磁Hertz矢量位求解得到散射電場的表達式為:

其中

相應的(14)式前三項的絕對值分布如圖5所示.

當b>4.1361a?0.0251λ且a<0.08λ時,其第一項的值遠大于高階項的值.取第一項代替整體散射場可得

對比沿z方向極化的磁偶極子的輻射場可以發現此時散射場就等于一磁偶極子的輻射場,相應的磁偶極子矢量為

圖5 垂直磁流源激勵情況下散射電場級數展開式 (a)前兩項絕對值最大值之比隨a和b的變化情況;(b)前三項絕對值沿x軸分布情況Fig.5.Series expansion of scattered electric fi eld with vertical excitation magnetic dipole:(a)Relationship of the ratio between maximum absolute value of the fi rst two terms with a and b;(b)absolute value distribution of the if rst three terms along x-axis in the plane z=b.

取上面4種情形滿足ρ值大于10的公共區間可得,當b>2.2875a?0.1043λ且a<0.08λ時,無論初始激勵源極化如何,均可采用散射場級數展開式的第一項代替整體散射場,即意味著可以將散射電場等效為電磁偶極子輻射場之和.

3 多金屬球散射體情況并矢格林函數分析

當有多個金屬小球散射體時,對每個金屬球而言,在求解總等效偶極子矢量時應考慮其余所有金屬球等效偶極子與初始激勵源的貢獻之和.將第j個散射體的等效電、磁偶極子矢量分別記為Pej和Pmj,其中j∈[0,1,···,N],N為散射金屬球的個數.初始激勵源視為0號散射體,其對應的電磁偶極子矢量大小即等于初始源的實際分量.分別將Pei和Pmi分解為相應的水平和垂直分量,即可得到每一個金屬小球的等效偶極子為

其中I為單位矩陣.因此可以得到此情況下的電并矢格林函數為

4 時間反演場數值求解結果

根據文獻[12]指出的電流源遠場時間反演場對應的格林函數與正向場對應格林函數的關系可知,初始電流源對應的時間反演格林函數為

其中imag[·]代表取虛部操作. 將(21)式代入(22)式即可得到時間反演格林函數的表達式.并由文獻[12]的結論可得相應的時間反演場表達式為

4.1 單散射體情況

首先考慮簡單情況,金屬球位于坐標原點,激勵源位于z軸上.假設金屬球半徑為a=0.03λ,所在的位置為(0,0,0),初始激勵電流源所在的位置為(0,0,h),其中h=0.1λ.相應的位置分布示意圖見圖6.如無特別說明,以下所有算例中初始激勵源的極化方向和所處位置均與此例一致.

圖6 初始源及散射體分布位置示意圖Fig.6.Distribution diagram of the excitation dipole and the scatterer.

根據(21)式和(22)式可得到此時的時間反演電并矢格林函數為

根據CST仿真和由(24)式計算得到的反演電場分布分別如圖7所示(以下所有場觀測面均默認為z=0.1λ平面).因為入射場的Ex和Ey分量遠大于散射場的相應分量,引入散射體并不會明顯改變其分布,因此以下只給出反演電場Ez分量的對比情況.

通過計算場值大小對應數值矩陣的相關性即可得到場解的相似度.矩陣的相關性Corr可以通過以下公式計算得到[22]:

其中xij和yij分別代表兩個矩陣中的不同元素;和為相應的平均值.

計算反演腔內反演電場的相似程度即可得知解析結果是否正確.因為反演場主要集中在±0.5λ的范圍之內,因此下面求解二者相似度時只考慮此范圍內的場.利用(25)式計算得到圖7(a)和圖7(b)中黑色框內反演電場的相似度為0.9817,可以發現二者幾乎完全一樣,這也說明了解析計算結果的正確性.同時圖7(c)給出了由CST和解析方法求解得到的反演場沿y=0直線的分布對比情況,可以看到,二者符合度非常高.

4.2 多散射體情況

利用所提方法計算反演場的最大優勢在于求解速度快,尤其是在散射體數量較多的時候,不需要進行復雜的邊值問題求解.待求解場景源與散射體的分布如圖8所示.金屬球所處的位置分別為r1=(d,0,0),r2=(0,0,0),r3=(?d,d,0),r4=(?d,?d,0),其中d=0.3λ.

圖7 反演電場Ez分布圖 (a)CST計算結果;(b)MATLAB解析計算結果;(c)兩者結果沿y=0直線分布對比Fig.7. Distribution of time reversal field Ez:(a)CST result;(b)theoretical result of MATLAB;(c)comparison of field distribution along y=0 line.

圖8 初始源及散射體分布位置示意圖Fig.8.Distribution diagram of the excitation dipole and the scatterers.

圖9 反演電場Ez分布 (a)CST計算結果;(b)MATLAB解析計算結果;(c)二者結果沿y=0直線分布對比Fig.9.Distribution of time reversal field Ez:(a)CST result;(b)theoretical result of MATLAB;(c)comparison of field distribution along y=0 line.

與單散射體情況類似,由CST計算得到的結果如圖9(a)所示,相應采用(21)—(23)式可以得到反演場分布如圖9(b)所示.同樣,根據(25)式可以計算得到此時解析解與仿真解的相似度為0.9348,二者相似度依舊非常高.同時圖9(c)給出的CST計算結果和解析結果沿y=0直線的分布對比圖也很好地表明了二者的相似程度.

加載金屬小球的最終目的是為了輔助實現源目標的超分辨率聚焦.綜合以上計算結果可以發現,解析方法計算結果與商業軟件CST仿真結果基本一致,這充分說明了解析推導的正確性.

圖10 初始源及散射體分布位置示意圖Fig.10.Distribution diagram of the excitation dipole and the scatterers.

圖11 反演電場Ez分布圖 (a)空間分布情況(MATLAB結果);(b)反演場幅值沿圖(a)虛線分布情況Fig.11.Distribution of time reversal field Ez:(a)Spatial distribution(MATLAB result);(b)field distribution along dotted line in Fig.(a).

下面通過一個具體算例來考察加載金屬小球之后源目標的聚焦分辨率情況.仿真模型的設置如圖10所示,三個小金屬球所處的位置分別為r1=(0,d,0),r2=(0,0,0),r3=(0,?d,0),其中d=0.3λ.相應反演場的求解結果如圖11所示.考慮到磁流源反演場與電流源反演場的差異性,其空間分布并不是呈現中間大周圍小的“sinc函數式”分布,而是呈現類似圖7所示的“8字形”分布,因此不能通過直接觀察整個空間場幅度大小分布來確定聚焦分辨率的大小.將一個“8字形”的分布看作一個整體,考察不同“8字形”整體最大場幅值的分布情況,即考察反演場幅度沿圖11(a)所示虛線的分布情況.相應結果如圖11(b)所示.可以看到,相距0.3λ位置處的反演場幅度值比最大幅值的1/2還小.以場幅值減小1/2作為可分辨的臨界標準,可得出結論,反演場的聚焦分辨率是小于λ/2的.

5 結 論

基于等效偶極子模型,提出了一種快速求解亞波長間距金屬小球陣列加載情況下時間反演場的快速求解方法.分析結果表明,在金屬球半徑及其與激勵源之間間距滿足一定關系時,近場區域的散射場能夠等效為電磁偶極子輻射場的疊加.在求解多散射體情況的散射場時,只要金屬小球之間間距與小球半徑大小之間也滿足前述關系,就可以采用該近似將所有金屬小球等效為相應的電磁偶極子.在利用該方法得到的并矢格林函數的基礎上,根據時間反演腔理論能夠快速地得到相應時間反演場的解析解.文章通過與由商業軟件CST求解得到的時間反演場分布結果進行的對比表明,兩種方法的求解結果數據符合度高達90%以上,這很好地驗證了本文所提出模型的正確性.同時結合所提出方法進行的計算結果表明,通過在源的近場加載金屬小球能夠實現0.3λ的源目標聚焦分辨率.該等效模型的建立為后續快速分析小金屬球散射體加載情況下的時間反演場分布提供了一種高效便捷的處理思路及理論指導.同時對后續高效研究如何進行合適的散射體加載以實現更高的聚焦分辨率提供了很好的分析工具.

[1]Parvulescu A,Clay C S 1965Radio Electron.Eng.29 223

[2]Fink M 1992IEEE Trans.Ultrason.Ferroelectr.Freq.Control39 555

[3]de Rosny J,Fink M 2007Phys.Rev.A76 065801

[4]Fang N,Liu Z,Yen T J,Zhang X 2003Opt.Express11 682

[5]Liu Z,Fang N,Yen T J,Zhang X 2003Appl.Phys.Lett.83 5184

[6]Grbic A,Eleftheriades G V 2004Phys.Rev.Lett.92 117403

[7]Pendry J B 2000Phys.Rev.Lett.85 3966

[8]Malyuskin O,Fusco V 2010IEEE Trans.Antenna.Propag.58 459

[9]Lemoult F,Lerosey G,de Rosny J,Fink M 2010Phys.Rev.Lett.104 203901

[10]Ourir A,Fink M 2014Phys.Rev.B89 115403

[11]Jouvaud C,Ourir A,de Rosny J 2014Appl.Phys.Lett.104 243507

[12]Carminati R,Pierrat R,de Rosny J,Fink M 2007Opt.Lett.32 3107

[13]Ioannidou M P,Skaropoulos N C,Chrissoulidis D P 1995J.Opt.Soc.Am.A12 1782

[14]Borghese F,Denti P,Toscano G,Sindoni O I 1979Appl.Opt.18 116

[15]Gouesbet G,Grehan G 1999J.Opt.A:Pure Appl.Opt.1 706

[16]Moneda A P,Chrissoulidis D P 2007J.Opt.Soc.Am.A24 3437

[17]Purcell E M,Pennypacker C R 1973Astrophys.J.186 705

[18]Moneda A P,Chrissoulidis D P 2007J.Opt.Soc.Am.A24 3437

[19]Fallahi A,Oswald B 2011IEEE Trans.Microwave Theory Tech.59 1433

[20]Harrington R F 2001Time-harmonic Electromagnetic Fields(New York:John Wiley&Sons)p293

[21]Wait J R 1960Geophysics25 649

[22]Rabiner L R,Gold B 1975Theory and Application of Digital Signal Processing(Englewood Clif f s,NJ:Prentice-Hall)p401