采用分部外推BCGS-FFT方法快速求解電磁散射問題

韓曉冰 張瀟 王露潔 周遠(yuǎn)國 任強(qiáng)

（1.西安科技大學(xué)通信與信息工程學(xué)院，西安 710054；2.北京航空航天大學(xué)電子信息工程學(xué)院，北京 100191）

引 言

近年來，電磁散射在地球物理勘探[1]、遙感[2]、微帶天線[3]、微波成像[4]、無損檢測(cè)[5]等民用和軍用領(lǐng)域的應(yīng)用愈發(fā)廣泛.在這些領(lǐng)域的電磁散射研究過程中，經(jīng)常需要考慮一種層狀介質(zhì)的模型.由于異常體的非均勻性，傳統(tǒng)的表面積分方程[6]是不適用的，這種情況下應(yīng)用體積分方程更為合適.解決非均勻物體電場(chǎng)積分方程(electric field integral equation,EFIE)最直接的方法是矩量法(method of moments,MoM)[7].傳統(tǒng)求解體積分方程的一類方法是將迭代方法與快速傅里葉變換(fast Fourier transform,FFT)結(jié)合，比較典型的有共軛梯度快速傅里葉變換(conjugate gradient fast Fourier transform, CG-FFT)方法、雙共軛梯度快速傅里葉變換(biconjugate gradient fast Fourier transform, BCG-FFT)方法、穩(wěn)定雙共軛快速傅里葉變換(stabilized biconjugate gradient fast Fourier transform, BCGS-FFT)方法等.由文獻(xiàn)[8-9]中幾個(gè)均勻背景下的算例可知，BCG-FFT方法比CGFFT方法更高效.BCGS-FFT方法[10]在求解均勻背景下大規(guī)模散射問題比BCG-FFT方法和CG-FFT方法精確度更高，效率更快[11].

但是在求解EFIE的過程中，我們需要得到空域格林函數(shù)，這個(gè)過程一般是由索末菲積分完成的.在索末菲積分中，盡管層狀介質(zhì)譜域格林函數(shù)是解析的[12]，但是這個(gè)積分通常是震蕩的、奇異的、收斂很慢甚至是不收斂的，并且用常規(guī)的積分方法是費(fèi)時(shí)費(fèi)力的.因此，很有必要加速索末菲積分的收斂速度.目前計(jì)算索末菲積分的方法主要有復(fù)平面積分法[13]、復(fù)鏡像法[14]、近似法[15]、最陡下降法[16]，以及對(duì)索末菲尾部積分應(yīng)用分部外推法[17-19]等等.

本文我們?cè)趶?fù)平面中采用一種橢圓路徑來避開格林函數(shù)的奇異點(diǎn)并用分部外推法加速索末菲尾部積分的收斂速度.這種方法的優(yōu)勢(shì)是在尾部積分部分沿實(shí)軸積分時(shí)貝塞爾函數(shù)只有實(shí)參，并且不需要確定格林函數(shù)極點(diǎn)的精確位置.該方法比復(fù)鏡像[20]等近似方法更精確可靠，比直接積分法[21]速度更快，并且不用耗時(shí)去提取奇異點(diǎn)[22].一般來講這種積分方法是普遍適用的，但應(yīng)用該方法時(shí)盡量不要讓源點(diǎn)與場(chǎng)點(diǎn)之間的水平距離大于100個(gè)波長，否則會(huì)使貝塞爾函數(shù)在初始區(qū)間中有過多的震蕩，影響最終的積分結(jié)果[23].在各類分部外推法中，目前使用較多的有Euler變換[24]、Aitken迭代變換[25]、加權(quán)平均[26]、Shanks變換[27]，以及一般Levin變換[28]結(jié)合W算法[29]等，文獻(xiàn)[17]中的數(shù)值測(cè)試部分對(duì)比了這些方法，因?yàn)樗髂┓品e分的漸進(jìn)現(xiàn)象，一般Levin變換結(jié)合W算法與加權(quán)平均算法是加速索末菲尾部積分最通用有效的方法.本文的工作主要分為三個(gè)部分：①在復(fù)平面中使用一種簡(jiǎn)單實(shí)用的橢圓路徑來避開格林函數(shù)的所有奇異點(diǎn).②將分部外推法應(yīng)用于傳統(tǒng)計(jì)算索末菲尾部積分的高斯積分法[30]中，以此來計(jì)算多層介質(zhì)并矢格林函數(shù)，填充并矢格林函數(shù)矩陣.③將分部外推法與BCGS-FFT方法結(jié)合求解層狀介質(zhì)中的電磁散射問題，提高計(jì)算效率.

1 EFIE

如圖1所示的分層介質(zhì)的電磁散射模型中共有n層，在第m層中嵌入非均勻電介質(zhì)物體.因此，在第j層中總電場(chǎng)Ej(r)等于入射電場(chǎng)與散射電場(chǎng)的和，即：

圖1 平面分層介質(zhì)中電磁散射的非均勻物體的典型幾何模型Fig.1 A typical geometric model of an inhomogeneous object with electromagnetic scattering in a plane layered medium

式中：r是空間中的位置矢量；是在j層沒有異常體，僅僅有分層介質(zhì)情況下的入射電場(chǎng).關(guān)于散射場(chǎng)我們可以用磁矢勢(shì)Ajm(r)和 標(biāo)量勢(shì)?jm(r)來表示：

這里磁矢勢(shì)Ajm(r)是由嵌入第m層的非均勻物體產(chǎn)生的感應(yīng)電流源J引起的，表達(dá)式為

式中： μj是第j層的介磁常數(shù)；V是異常體體積；Gjm(r,r′) 是傳統(tǒng)形式的并矢格林函數(shù)[12]；J可以表示為J=iωχD， 其中是異常體的對(duì)比度函數(shù)，是復(fù)介電常數(shù)， εm與σm分 別是第m層介電常數(shù)、磁導(dǎo)率電通量通過洛倫茲條件，我們可以把磁矢勢(shì)Ajm(r) 和 標(biāo)量勢(shì)?jm(r)聯(lián)系起來，即：

將式(2)～(4)代入式(1)中，當(dāng)r∈V(也就是m=j)時(shí)，我們可以得到EFIE：

這里我們可以引入一個(gè) γ算子，即：

可以看出，只要電通量Dm(r)確定，我們就可以通過式(2)～(4)求出第j層任意位置的散射場(chǎng).我們的方案是使用Zwamborn和Berge提出的屋頂基函數(shù)弱離散EFIE方程的方法[31]，并用BCGS代替CG方法來求電通量Dm(r).

2 分層介質(zhì)格林函數(shù)的計(jì)算

在式(3)中，我們需要得到空域并矢格林函數(shù)Gjm(r,r′)，Gjm(r,r′)一般可被寫為[12]

譜域格林函數(shù)的封閉形式可以由傳輸線理論推導(dǎo)得出[2]，空域格林函數(shù)是由譜域格林函數(shù)經(jīng)過x和y方向上傅里葉逆變換得到的，這與索末菲積分或者漢克爾變換是等價(jià)的，即：

式中：n為貝塞爾函數(shù) Jn(x)的 階數(shù)；和ky為 波矢量k在x和y方向的分量.

一般來說，傳統(tǒng)方法多采用在復(fù)平面區(qū)域沿著實(shí)軸進(jìn)行積分，這樣可以避開貝塞爾函數(shù)復(fù)雜的參數(shù)計(jì)算.但這類方法需要先精確定位平面波的極點(diǎn)，再計(jì)算剩余路徑的積分值.這種做法在只有一層結(jié)構(gòu)時(shí)比較容易做到，但是在多層結(jié)構(gòu)時(shí)會(huì)變得非常的耗時(shí).因此，本文采用一種橢圓路徑來避開所有極點(diǎn)[32]，當(dāng)然路徑的選擇并不是唯一的[33]，但是橢圓是簡(jiǎn)單并且靈活的路徑，這可以使我們不用去計(jì)算極點(diǎn)的位置，這種路徑可以最小化貝塞爾函數(shù)的振蕩和表面波極點(diǎn)的影響[34].格林函數(shù)的極點(diǎn)在的范圍內(nèi)，格林函數(shù)在這之間沿著實(shí)軸是無損的，沿著虛軸是指數(shù)變化的.因此，我們可以得到一條如圖2的橢圓積分路徑.橢圓積分路徑的中心在處，橢圓的短半軸平行于虛軸，與源點(diǎn)和場(chǎng)點(diǎn)之間的水平距離 ρ成反比關(guān)系.這種路徑避免了因貝塞爾函數(shù)自變量的虛部增長帶來的指數(shù)變化引起的數(shù)值問題.在數(shù)值積分時(shí)，可以先采用常規(guī)方法將橢圓參數(shù)化，為了保證積分的精度，積分點(diǎn)的數(shù)量應(yīng)該隨著場(chǎng)點(diǎn)與源點(diǎn)的水平距離 ρ的增大而增多，隨后采用高斯積分或其他積分方法.傳統(tǒng)的BCGS-FFT做法是在尾部積分部分繼續(xù)使用高斯積分，每得到一小段的積分就與前面所得到的積分和來對(duì)比，當(dāng)這段積分值小于積分和與截止誤差的積時(shí)停止積分.在本文中，我們?cè)谖膊糠e分部分使用高效計(jì)算的Levin變換[28]分部外推算法.

圖2 索末菲積分路徑Fig.2 Sommerfeld integral path

3 分部外推算法

關(guān)于各種具體的分部外推法在格林函數(shù)尾部積分段的表達(dá)形式，Michalski已經(jīng)做了詳細(xì)的說明[17].本文只介紹用到的Levin變換和Sidi的W算法.一般Levin變換被定義為一個(gè)模型序列，即：

式中： ωn為殘差估計(jì)系數(shù)；k為Levin變換的階數(shù)；ci為 未知系數(shù)；序列 {xn}是插值點(diǎn)的輔助序列，即圖2中的 ξ?1，ξ0,···.

Sn為前n項(xiàng)積分值的和，即：

式中，ui為第i段積分的值，u0即 為圖2中 ξ?1與 ξ0之間的積分值.令 ωn=un，通過式(9)我們可以看出，Levin變換的模型序列里的未知量有待求極限S和k個(gè)系數(shù)c?1,c0,···,ck?1.因此，只需要求出k+1個(gè)前n項(xiàng)積分和 {Sn,Sn+1,···,Sn+k}， 即可構(gòu)造一個(gè)擁有k+1個(gè)方程的方程組，得益于克萊姆法則，我們不需要求出未知系數(shù)ci， 就可以獲得所求的極限S：

式中，

可以看出，這種Levin變換的劣勢(shì)是非遞歸的，即每次設(shè)置k的值時(shí)都必須從新開始計(jì)算.為了解決這一問題，Sidi使用一種W算法使上面的算法變?yōu)檫f歸的[29].W算法通過引入一個(gè)除法差分算子δk：

這里的 δ0(un)≡u(píng)n，這樣式(9)可以改寫為

可以發(fā)現(xiàn)，式(14)等號(hào)右邊為k?1次 的關(guān)于的多項(xiàng)式，所以，在式(14)兩邊同時(shí)應(yīng)用δk算 子，可以得到

因此，利用 δk算子的線性性質(zhì)，我們可以把積分極限S表示為

4 數(shù)值算例

在數(shù)值算例部分，我們模擬了三種典型的平面分層結(jié)構(gòu)中嵌入異常體的電磁散射情況來驗(yàn)證上述方法.隨后，將本文提出的分部外推BCGS-FFT方法與傳統(tǒng)BCGS-FFT方法在數(shù)值精度與計(jì)算時(shí)間上做了對(duì)比.為了方便起見，將三種算例的源都設(shè)為一個(gè)單位電偶極子，極化方向?yàn)?1, 1, 1)，放置在(0,0, ?0.5)的位置.所有BCGS的迭代截止誤差和索末菲積分的截止誤差都設(shè)為10?5.所有算例的仿真都在Intel(R)Core(TM)i5-8250U的CPU以及16 G內(nèi)存的工作平臺(tái)上進(jìn)行.

4.1 半空間模型中嵌入一個(gè)圓柱散射體

算例1：如圖3，考慮一個(gè)直徑d=1.6m 、高h(yuǎn)=1.6m的圓柱體嵌入到半空間的下層.半空間模型的分界面為z=0平面，上層空間為真空，下層空間的介質(zhì)為相對(duì)介電常數(shù)、磁導(dǎo)率、電導(dǎo)率分別為 εr2=4.0、μr2=1.0、 σ2=0.01S/m的類似粘性干土物質(zhì).圓柱散射體的電性參數(shù)設(shè)置為 εr=7.0 、 μr=1.0、 σ=0.001S/m的類似濕花崗巖物質(zhì).圓柱散射體的中心在(0, 0, 2.2)位置.

圖3 直徑d = 1.6 m，高h(yuǎn) = 1.6 m的圓柱散射體嵌入半空間的下層示意圖Fig.3 A cylindrical scatterer with a diameter of d = 1.6 m and a height of h = 1.6 m embedded in a semi-spatial background medium

本次仿真的工作頻率設(shè)置為200 MHz，異常體剖分如圖4所示，離散成 50×50×50的立方體單元，每個(gè)正方體邊長為0.04 m.因此每個(gè)波長采樣點(diǎn)數(shù)為37.5，設(shè)置61個(gè)接收點(diǎn).圖5和圖6展示了兩種方法計(jì)算出的接收點(diǎn)處散射電場(chǎng)與磁場(chǎng)的精度對(duì)比，可以發(fā)現(xiàn)新方法可以達(dá)到與傳統(tǒng)方法相同的計(jì)算精度.但是，在計(jì)算機(jī)內(nèi)存消耗相同的情況下，分部外推BCGS-FTT方法與傳統(tǒng)BCGS-FFT方法的計(jì)算耗時(shí)分別為11 238.3 s和17 819.3 s，加速后的BCGS-FFT可以節(jié)省36.7%的時(shí)間.

圖4 剖分圓柱體網(wǎng)格區(qū)域以及設(shè)置61個(gè)接收點(diǎn)的位置坐標(biāo)Fig.4 Sectional cylindrical grid area and the coordinates of 61 receiving points

圖5 算例1中接收點(diǎn)處的散射電場(chǎng)對(duì)比Fig.5 Comparison of the scattered electric field at the receiving points in example 1

圖6 算例1中接收點(diǎn)處的散射磁場(chǎng)對(duì)比Fig.6 Comparison of the scattered magnetic field at the receiving points in example 1

4.2 三層結(jié)構(gòu)中嵌入一個(gè)球體

算例2：如圖7所示，我們考慮一個(gè)直徑d=10m的球體嵌入到三層空間的中間層中.第一層介質(zhì)為真空，中間層為 εr2=4.0、 μr2=1.0、 σ2=0.01S/m類似粘性干土的物質(zhì)，第三層空間為 εr3=3.0 、μr3=1.0 、 σ3=0.01S/m類似干砂的物質(zhì).球體還是類似濕花崗巖的物質(zhì)，其電性參數(shù)為 εr=7.0、 μr=1.0、σ=0.001S/m.第一層、第二層的底部分界面分別在z1=0、z2=20.0處.球體的中心在(0, 0, 10)處.

圖7 直徑為d = 10 m的球體嵌入三層介質(zhì)空間中間層示意圖Fig.7 A spherical scatterer with a diameter of d = 10 m embedded in the middle layer of a three-layer medium space

本次仿真工作頻率設(shè)置為20 MHz，剖分結(jié)果如圖8所示，每個(gè)波長采樣點(diǎn)數(shù)為75，設(shè)置61個(gè)接收點(diǎn).圖9和圖10展示了兩種方法計(jì)算出的接收點(diǎn)處散射電場(chǎng)與磁場(chǎng)的精度對(duì)比，可以發(fā)現(xiàn)新方法達(dá)到了傳統(tǒng)方法的計(jì)算精度.但是，分部外推BCGS-FFT方法與傳統(tǒng)的BCGS-FFT方法的計(jì)算耗時(shí)分別為9 581.4 s和12 065.4 s，加速后的BCGS-FFT可以節(jié)省20.6%的計(jì)算耗時(shí).

圖8 剖分球體網(wǎng)格區(qū)域以及設(shè)置61個(gè)接收點(diǎn)的位置坐標(biāo)Fig.8 Sectional spheriform grid area and the coordinates of 61 receiving points

圖9 算例2中接收點(diǎn)處的散射電場(chǎng)對(duì)比Fig.9 Comparison of the scattered electric field at the receiving points in example 2

圖10 算例2中接收點(diǎn)處的散射磁場(chǎng)對(duì)比Fig.10 Comparison of the scattered magnetic field at the receiving points in example 2

4.3 七層結(jié)構(gòu)中嵌入一個(gè)雙層結(jié)構(gòu)正方體

算例3：如圖11所示，我們考慮一個(gè)邊長l=60m的雙層正方體嵌入到七層空間的第六層中.第一層為真空，第二層到第七層的介質(zhì)分別為：

圖11 上層高h(yuǎn)1 = 30 m，下層h2 = 30 m的雙層正方體嵌入七層介質(zhì)空間第六層示意圖Fig.11 The upper level of the double cube with h1 = 30 m and the lower level with h2 = 30 m, and the double cube is embedded into the sixth layer of a seven-layer medium space

雙層正方體的上下層電性參數(shù)分別為：

第一層至第六層的底部分界面分別在z1=0,z2=4.0,z3=8.0,z4=12.0,z5=16.0,z6=86.0處.雙層正方體的中心在(0, 0, 51)的位置.

本次仿真工作頻率設(shè)置為2 MHz.分部外推BCGS-FFT算法與傳統(tǒng)BCGS-FFT算法都將雙層正方體離散成 6 0×60×60的立方體單元，每個(gè)正方體邊長為1 m，如圖12所示.因此每個(gè)波長采樣點(diǎn)數(shù)為150，設(shè)置81個(gè)接收點(diǎn).圖13和圖14展示了兩種方法計(jì)算出的接收點(diǎn)處散射電場(chǎng)與磁場(chǎng)的精度對(duì)比.可以發(fā)現(xiàn)新方法與傳統(tǒng)方法達(dá)到的計(jì)算精度相同.但是，分部外推BCGS-FTT方法與傳統(tǒng)的BCGS-FFT方法的計(jì)算耗時(shí)分別為22 651.9 s和36 128.9 s，加速后的BCGS-FFT算法可以節(jié)省37.3%的計(jì)算耗時(shí).

圖12 剖分雙層正方體網(wǎng)格區(qū)域以及設(shè)置81個(gè)接收點(diǎn)的位置坐標(biāo)Fig.12 Divide the spheriform grid area and set the coordinates of 81 receiving points

圖13 算例3中接收點(diǎn)處的散射電場(chǎng)對(duì)比Fig.13 Comparison of the scattered electric field at the receiving points in example 3

圖14 算例3中接收點(diǎn)處的散射磁場(chǎng)對(duì)比Fig.14 Comparison of the scattered magnetic field at the receiving points in example 3

5 結(jié) 論

本文采用了分部外推法計(jì)算空域格林函數(shù)的尾部積分，且與BCGS-FFT結(jié)合以便加速計(jì)算非均勻散射體嵌入多層空間的電磁散射問題，并與傳統(tǒng)BCGS-FFT方法進(jìn)行比較.通過三個(gè)不同類型異常體的算例驗(yàn)證了該算法的可靠性.將分部外推BCGSFFT算法與傳統(tǒng)BCGS-FFT的計(jì)算結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了算法精確度.結(jié)果表明，本文提出的分部外推BCGS-FFT算法適用于計(jì)算異常體在層狀結(jié)構(gòu)中的電磁散射問題.在精度相同的情況下，三個(gè)算例分別能節(jié)省36.7%、20.6%、37.3%的計(jì)算耗時(shí).

值得注意的是，本文公式推導(dǎo)僅適用于各向同性介質(zhì)，并且不考慮磁介質(zhì)材料.此外，文中提出的分部外推BCGS-FFT算法僅適用于散射體嵌入單一層中，暫不適用于跨層散射體.諸多待解決問題將成為我們今后的研究方向.