龔志雙王秉中王任

(電子科技大學(xué)應(yīng)用物理研究所,成都 610054)

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

1 引 言

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

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

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

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

2 偶極子源激勵(lì)情況散射場(chǎng)求解

2.1 水平電偶極子源激勵(lì)情況

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

為了對(duì)凋落波進(jìn)行操作使其能轉(zhuǎn)化為傳輸波,需要所加載的散射體尺寸與相應(yīng)的空間頻率相對(duì)應(yīng),而凋落波對(duì)應(yīng)較高的空間頻率分量,因此要求所加載的散射體尺寸相對(duì)自由空間波長較小.所以本文主要考慮加載散射體尺寸遠(yuǎn)小于自由空間波長的情況(即ka<1).同時(shí),因?yàn)榈蚵洳S傳輸距離呈指數(shù)衰減,需要在其衰落至太小之前對(duì)其進(jìn)行操作,也即意味著需要在源的近場(chǎng)區(qū)進(jìn)行加載.故本文只討論源與散射體之間間距小于半波長(即kb<π)的情形.

根據(jù)文獻(xiàn)[11]中的結(jié)果可以得到圖1所示場(chǎng)景散射場(chǎng)的表達(dá)式為

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

其中cen和den分別為

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

圖2 水平電流源激勵(lì)情況下散射電場(chǎng)級(jí)數(shù)展開式 (a)前兩項(xiàng)絕對(duì)值最大值之比隨a和b的變化情況;(b)前三項(xiàng)絕對(duì)值沿x軸分布情況(其中a=0.03λ,b=0.3λ,觀測(cè)面為z=b平面.以下所有類似場(chǎng)幅度值分布圖所對(duì)應(yīng)的計(jì)算場(chǎng)景與此圖相同)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)式得到第一項(xiàng)的表達(dá)式為

對(duì)比沿x方向電偶極子及沿y方向磁偶極子的輻射場(chǎng)[20]可以發(fā)現(xiàn),此散射場(chǎng)等于一電偶極子和一磁偶極子輻射場(chǎng)的疊加.相應(yīng)的電磁偶極子矢量分別為

2.2 垂直電偶極子源激勵(lì)情況

對(duì)于圖1所示的模型,當(dāng)激勵(lì)源的極化方向?yàn)檠貁軸時(shí),電偶極子強(qiáng)度可表示為pe=Ielez.根據(jù)對(duì)稱性并且考慮小金屬球情況可得散射場(chǎng)表達(dá)式為

(6)式前三項(xiàng)絕對(duì)值分布的計(jì)算結(jié)果如圖3所示.仍以ρ>10為判斷標(biāo)準(zhǔn),當(dāng)b>1.0429a+0.1416λ

圖3 垂直電流源激勵(lì)情況下散射電場(chǎng)級(jí)數(shù)展開式 (a)前兩項(xiàng)絕對(duì)值最大值之比隨a和b的變化情況;(b)前三項(xiàng)絕對(duì)值沿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λ時(shí),可以只取(6)式的第一項(xiàng)代表整體散射場(chǎng),其表達(dá)式為

對(duì)比沿z方向極化電偶極子的輻射場(chǎng)可知,此時(shí)散射場(chǎng)等于一電偶極子的輻射場(chǎng).相應(yīng)的等效電偶極子矢量為

2.3 水平磁偶極子源激勵(lì)情況

當(dāng)激勵(lì)源為磁偶極子并且極化方向?yàn)檠豿軸時(shí),磁偶極子矢量可表示為pm=Imlex,其中Im為磁流強(qiáng)度.利用前面計(jì)算得到的電偶極子激勵(lì)情況的結(jié)果并根據(jù)電與磁的對(duì)偶性可得此時(shí)散射場(chǎng)的表達(dá)式為:

(10)式前三項(xiàng)絕對(duì)值分布計(jì)算結(jié)果如圖4所示.類似可得當(dāng)b>0.9938a+0.1248λ且a<0.08λ時(shí),其第一項(xiàng)的值遠(yuǎn)大于高階項(xiàng)的值.取(10)式的第一項(xiàng)可得散射電場(chǎng)的表達(dá)式為:

對(duì)比沿x方向極化磁偶極子及沿y方向極化電偶極子的輻射場(chǎng)可以發(fā)現(xiàn),此時(shí)散射場(chǎng)等于一磁偶極子與一電偶極子輻射場(chǎng)之和,相應(yīng)的偶極子矢量分別為:

圖4 水平磁流源激勵(lì)情況下散射電場(chǎng)級(jí)數(shù)展開式 (a)前兩項(xiàng)絕對(duì)值最大值之比隨a和b的變化情況;(b)前三項(xiàng)絕對(duì)值沿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 垂直磁偶極子源激勵(lì)情況

當(dāng)激勵(lì)源為一沿z向極化的磁偶極子時(shí),磁偶極子強(qiáng)度矢量為pm=Imlez.利用對(duì)稱性可以由磁Hertz矢量位求解得到散射電場(chǎng)的表達(dá)式為:

其中

相應(yīng)的(14)式前三項(xiàng)的絕對(duì)值分布如圖5所示.

當(dāng)b>4.1361a?0.0251λ且a<0.08λ時(shí),其第一項(xiàng)的值遠(yuǎn)大于高階項(xiàng)的值.取第一項(xiàng)代替整體散射場(chǎng)可得

對(duì)比沿z方向極化的磁偶極子的輻射場(chǎng)可以發(fā)現(xiàn)此時(shí)散射場(chǎng)就等于一磁偶極子的輻射場(chǎng),相應(yīng)的磁偶極子矢量為

圖5 垂直磁流源激勵(lì)情況下散射電場(chǎng)級(jí)數(shù)展開式 (a)前兩項(xiàng)絕對(duì)值最大值之比隨a和b的變化情況;(b)前三項(xiàng)絕對(duì)值沿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的公共區(qū)間可得,當(dāng)b>2.2875a?0.1043λ且a<0.08λ時(shí),無論初始激勵(lì)源極化如何,均可采用散射場(chǎng)級(jí)數(shù)展開式的第一項(xiàng)代替整體散射場(chǎng),即意味著可以將散射電場(chǎng)等效為電磁偶極子輻射場(chǎng)之和.

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

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

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

4 時(shí)間反演場(chǎng)數(shù)值求解結(jié)果

根據(jù)文獻(xiàn)[12]指出的電流源遠(yuǎn)場(chǎng)時(shí)間反演場(chǎng)對(duì)應(yīng)的格林函數(shù)與正向場(chǎng)對(duì)應(yīng)格林函數(shù)的關(guān)系可知,初始電流源對(duì)應(yīng)的時(shí)間反演格林函數(shù)為

其中imag[·]代表取虛部操作. 將(21)式代入(22)式即可得到時(shí)間反演格林函數(shù)的表達(dá)式.并由文獻(xiàn)[12]的結(jié)論可得相應(yīng)的時(shí)間反演場(chǎng)表達(dá)式為

4.1 單散射體情況

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

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

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

根據(jù)CST仿真和由(24)式計(jì)算得到的反演電場(chǎng)分布分別如圖7所示(以下所有場(chǎng)觀測(cè)面均默認(rèn)為z=0.1λ平面).因?yàn)槿肷鋱?chǎng)的Ex和Ey分量遠(yuǎn)大于散射場(chǎng)的相應(yīng)分量,引入散射體并不會(huì)明顯改變其分布,因此以下只給出反演電場(chǎng)Ez分量的對(duì)比情況.

通過計(jì)算場(chǎng)值大小對(duì)應(yīng)數(shù)值矩陣的相關(guān)性即可得到場(chǎng)解的相似度.矩陣的相關(guān)性Corr可以通過以下公式計(jì)算得到[22]:

其中xij和yij分別代表兩個(gè)矩陣中的不同元素;和為相應(yīng)的平均值.

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

4.2 多散射體情況

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

圖7 反演電場(chǎng)Ez分布圖 (a)CST計(jì)算結(jié)果;(b)MATLAB解析計(jì)算結(jié)果;(c)兩者結(jié)果沿y=0直線分布對(duì)比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 反演電場(chǎng)Ez分布 (a)CST計(jì)算結(jié)果;(b)MATLAB解析計(jì)算結(jié)果;(c)二者結(jié)果沿y=0直線分布對(duì)比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計(jì)算得到的結(jié)果如圖9(a)所示,相應(yīng)采用(21)—(23)式可以得到反演場(chǎng)分布如圖9(b)所示.同樣,根據(jù)(25)式可以計(jì)算得到此時(shí)解析解與仿真解的相似度為0.9348,二者相似度依舊非常高.同時(shí)圖9(c)給出的CST計(jì)算結(jié)果和解析結(jié)果沿y=0直線的分布對(duì)比圖也很好地表明了二者的相似程度.

加載金屬小球的最終目的是為了輔助實(shí)現(xiàn)源目標(biāo)的超分辨率聚焦.綜合以上計(jì)算結(jié)果可以發(fā)現(xiàn),解析方法計(jì)算結(jié)果與商業(yè)軟件CST仿真結(jié)果基本一致,這充分說明了解析推導(dǎo)的正確性.

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

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

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

5 結(jié) 論

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

[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