快速計(jì)算一維分層粗糙面之間金屬目標(biāo)復(fù)合散射的互耦迭代算法

姬偉杰 童創(chuàng)明②

①(空軍工程大學(xué)導(dǎo)彈學(xué)院 三原 713800)

②(毫米波國(guó)家重點(diǎn)實(shí)驗(yàn)室 南京 210096)

1 引言

粗糙面和目標(biāo)的復(fù)合散射在雷達(dá)探測(cè)、目標(biāo)識(shí)別和遠(yuǎn)程遙感等領(lǐng)域有著重要的應(yīng)用。近年來(lái)，計(jì)算粗糙面和目標(biāo)的混合場(chǎng)的數(shù)值解法引起了諸多學(xué)者的興趣。例如：Johnson[1]曾用“四路徑”模型計(jì)算無(wú)限大介質(zhì)平板上介質(zhì)目標(biāo)的散射，文獻(xiàn)[2]利用GFBM/SAA研究低掠角入射時(shí)動(dòng)態(tài)分形粗糙面和2維艦船目標(biāo)的雙站散射模擬， 文獻(xiàn)[3,4]利用互耦迭代方法計(jì)算了金屬粗糙面與上方目標(biāo)的差場(chǎng)散射，文獻(xiàn)[5]等用矩量法計(jì)算了粗糙面與上方或下方目標(biāo)的復(fù)合散射，文獻(xiàn)[6]等基于互易定理和矩量法研究了平面上方2維介質(zhì)目標(biāo)對(duì)高斯波束的電磁散射，文獻(xiàn)[7]等研究了粗糙海面上方金屬目標(biāo)復(fù)合散射的解析-數(shù)值混合算法。

在實(shí)際自然界中，分層媒質(zhì)粗糙面更具有實(shí)際意義，比如覆蓋植被或雪的地面?？赡苡捎谟?jì)算量過于龐大，大多數(shù)學(xué)者都研究單層粗糙面和目標(biāo)的復(fù)合散射，只有少數(shù)學(xué)者研究了更具有實(shí)際意義的分層粗糙面和目標(biāo)復(fù)合散射。例如：文獻(xiàn)[8]研究了雙層隨機(jī)粗糙面和地下目標(biāo)的復(fù)合散射。文獻(xiàn)[9]研究了分界面是平面的分層媒質(zhì)上方目標(biāo)的散射。文獻(xiàn)[10]用FDTD研究了雙層粗糙面上方金屬目標(biāo)的散射，文獻(xiàn)[11]中用EBCM計(jì)算了分層粗糙面和下方目標(biāo)的復(fù)合散射，但是該算法只限于計(jì)算微粗糙面和小目標(biāo)的復(fù)合散射(例如半徑為0.16λ的圓柱)，不能計(jì)算中等程度粗糙面與較大目標(biāo)的復(fù)合散射。

由于不受粗糙面參數(shù)和目標(biāo)尺寸的限制，因此，數(shù)值算法的研究顯得更有意義。以前的研究都將目標(biāo)和分層粗糙面一起看成一個(gè)組合散射體，電磁積分方程離散后產(chǎn)生了巨大的未知量，需要龐大的計(jì)算量和存儲(chǔ)量，對(duì)計(jì)算機(jī)硬件和程序員編程技巧提出了很高的要求，這在很大程度上限制了數(shù)值算法的應(yīng)用。

在分層粗糙面與目標(biāo)的復(fù)合散射計(jì)算中，粗糙面的未知量占絕大多數(shù)，而目標(biāo)的未知量只占一小部分。因此如何有效計(jì)算粗糙面的電磁散射是一個(gè)很重要的問題。近年來(lái)發(fā)展的前后向迭代法(FBM)[12]對(duì)粗糙面這類散射問題有很好的收斂性，在文獻(xiàn)[13]中，Moss C D等人進(jìn)一步將其與譜加速(SAA)結(jié)合并應(yīng)用于計(jì)算分層粗糙面的電磁散射特性。本文結(jié)合FBM和Bi-CG，基于物理散射機(jī)理，提出一種計(jì)算雙層介質(zhì)粗糙面及下方金屬目標(biāo)復(fù)合電磁散射的快速互耦迭代數(shù)值算法(CCIA)。建立了目標(biāo)與分層粗糙面的耦合積分方程，將其分解為含有兩個(gè)方程的方程組，對(duì)分層粗糙面的散射方程用FBM求解，目標(biāo)的散射方程用雙共軛梯度法(Bi-CG)[14]求解，粗糙面與目標(biāo)的耦合作用通過更新兩方程的激勵(lì)項(xiàng)來(lái)實(shí)現(xiàn)，每次迭代過程中，粗糙面的激勵(lì)項(xiàng)中不僅有錐形入射波還有目標(biāo)對(duì)其的散射波；同樣，目標(biāo)的激勵(lì)項(xiàng)中也包括上層粗糙面的透射波和下層粗糙面的反射波。由于FBM和Bi-CG的計(jì)算量均為O(N2)，因此該算法的計(jì)算量為O(N2)，而傳統(tǒng)的矩量法(如高斯消元法)的計(jì)算量為O(N3)，因此計(jì)算效率大大提高。經(jīng)過計(jì)算，該算法具有良好的收斂性，當(dāng)目標(biāo)趨于無(wú)窮小時(shí)復(fù)合散射系數(shù)與只有分層粗糙面時(shí)的散射系數(shù)相吻合，驗(yàn)證了算法的有效性。最后，討論了目標(biāo)尺寸與深度變化對(duì)復(fù)合散射系數(shù)的影響。

2 基本理論

2.1 電場(chǎng)積分方程(EFIE)

設(shè)無(wú)限長(zhǎng)金屬圓柱目標(biāo)位于雙層粗糙面中間，如圖1所示。ψin為入射場(chǎng)，ψl為區(qū)域l中的場(chǎng)(l=0,1,2)。當(dāng)入射波為TE和TM波時(shí)，ψ分別表示電場(chǎng)和磁場(chǎng)(本文只考慮TE波入射)，在各區(qū)域中應(yīng)用格林定理：

n表示分界面表面法向向量。邊界條件為

2.2 耦合積分方程組

設(shè)粗糙面長(zhǎng)度為L(zhǎng)，離散密度Δx，離散總數(shù)為N，粗糙面Sl的高度用fl(xn)來(lái)表示。目標(biāo)表面用Δx0均勻離散，離散總數(shù)為M，目標(biāo)表面用z0(xn)表示。結(jié)合邊界條件，式(1)-式(3)整理為

整理式(4)便可得到分層粗糙面與下方金屬目標(biāo)耦合矩陣方程組：

其中式(5a)是求解分層粗糙面表面電流的方程，而式(5b)則是求解目標(biāo)表面電流的方程。該方程組表明粗糙面不僅受到初始入射波的影響，同時(shí)還要受到目標(biāo)對(duì)其散射的影響，目標(biāo)則只受到分層粗糙面對(duì)其散射的影響。=Cl?U，表示目標(biāo)對(duì)第l個(gè)0分界面的散射。ψSur=E(l)?Ul+F(l)?ψl，表示上層分界面的透射波與下層分界面反射波之和。

2.3 用CCIA求解耦合積分方程組

式(5)若直接用傳統(tǒng)矩量法求解則計(jì)算量為O(N3)，需要耗費(fèi)大量的時(shí)間。這里介紹一種新的迭代算法來(lái)快速求解該方程組，其中第i步的迭代方程如下：

其中Z0和Z1分別為分層粗糙面與目標(biāo)的阻抗矩陣，V0(i)表示第i次迭代過程中粗糙面的總?cè)肷洳?，包括錐形入射波和目標(biāo)對(duì)粗糙面的散射波。V1(i)表示第i次迭代過程中目標(biāo)總的入射波，包括上層粗糙面的透射波與下層粗糙面的反射波。

對(duì)式(6a)，采用FBM求解，將各子矩陣分解為上、下和對(duì)角矩陣，分別用U、D和L表示。例如，對(duì)上標(biāo)為(l, l, l)的子矩陣有

將未知向量分解為前后向分量，Ul=+，ψl=+，這里，和是前向分量，和是后向分量。于是，式(6a)中前向電流迭代公式寫為

類似地可得到后向電流的迭代公式。由迭代求解，第i次迭代的未知量為迭代算法以=0，=0，=0，=0為初始值，迭代至指定收斂精度。

用Bi-CG解式(6b)，F(xiàn)BM和Bi-CG求解均為循環(huán)迭代過程，稱為內(nèi)循環(huán)，設(shè)其各自的收斂精度為γ=10?4。而式(6a)和式(6b)的互耦迭代過程稱為外循環(huán)。由于FBM的計(jì)算量為O()(NSur為粗糙面的未知量個(gè)數(shù))，Bi-CG的計(jì)算量為O()(NTar為目標(biāo)的未知量個(gè)數(shù))，外循環(huán)迭代一次計(jì)算量為O()(N為總的未知量個(gè)數(shù))，因此該算法總Total的計(jì)算量為O()，較之傳統(tǒng)矩量法O()計(jì)算效率大為提高。

定義第n步外循環(huán)迭代誤差為

其中Ztotal和Itotal分別表示分層粗糙面與目標(biāo)總的阻抗矩陣和表面電流向量。從Itotal(n)=0開始迭代計(jì)算，此時(shí)粗糙面的入射波只有錐形波。計(jì)算得到粗糙面的表面電流分布，代入式(6b)右端激勵(lì)項(xiàng)計(jì)算目標(biāo)的總激勵(lì)，此時(shí)目標(biāo)的入射波包括上層粗糙面的透射波和下層粗糙面對(duì)反射波。計(jì)算得到目標(biāo)的表面電流再代入式(6a)右端更新粗糙面的激勵(lì)，如此反復(fù)迭代直至達(dá)到指定收斂精度。

2.4 錐形入射波和雙站散射系數(shù)

為減小人為截?cái)嘁鸬恼`差，采用Thorsos錐形波，表達(dá)式為[15]

其中θi為入射角(相對(duì)于z軸逆時(shí)針方向)，k表示上半空間的波數(shù)，g是波束寬度。

定義雙站散射系數(shù)[5]為

其中θ為散射角，ψs(θo, θi)為散射場(chǎng)。

3 數(shù)值計(jì)算和討論

本文中所有分界面均采用高斯粗糙面，長(zhǎng)度為L(zhǎng)，相關(guān)長(zhǎng)度為li(i=1,2)，均方根高度為hi(i=1,2)，i=1,2分別表示上層和下層分界面，計(jì)算過程中應(yīng)避免上層分界面與下層分界面相重疊，金屬目標(biāo)位于兩層粗糙面之間。媒質(zhì)的相對(duì)介電常數(shù)為εr1，εr2，上層介質(zhì)厚度為d。求解時(shí)采用脈沖基函數(shù)離散粗糙面表面電流分布，剖分單元為Δx=λ/10，圓柱半徑為R，深度dp為目標(biāo)與上層粗糙面的距離，xp為目標(biāo)中心的水平位置，剖分單元總數(shù)為120。如無(wú)特殊說明，文中參數(shù)的設(shè)置如下：εr1=4.0+0.01i ，εr2=7.0，d=6.0λ，R=1.0λ，dp=2.0λ，xp=0.0λ，θi=30o，g=L/4 ，L=40λ，h1=h2=0.1λ，l1=l2=1.0λ，文中所有結(jié)果均是對(duì)50個(gè)樣本取均值得到的。所用計(jì)算機(jī)配置為：主頻：1.65 GHz ，內(nèi)存：1 GB。

算法有效性的驗(yàn)證。令目標(biāo)尺寸趨于零，比較此時(shí)的散射系數(shù)與由矩量法得到的相同介質(zhì)參數(shù)、粗糙度參數(shù)對(duì)應(yīng)的分層粗糙面(粗糙面之間無(wú)目標(biāo))電磁散射系數(shù)的差異即可驗(yàn)證該算法的有效性。圖2所示為圓柱目標(biāo)半徑退化10?4λ時(shí)的散射系數(shù)與只有粗糙面時(shí)的散射系數(shù)相比較。 由圖可知，目標(biāo)尺寸趨于零時(shí)的散射系數(shù)與無(wú)目標(biāo)時(shí)的結(jié)果吻合得很好，驗(yàn)證了算法的有效性。

算法的計(jì)算精度與收斂性。其它參數(shù)不變，令上下層分界面均方根高度分別取0.1 λ，0.3 λ，0.5 λ。圖3 給出了外循環(huán)迭代誤差τ隨迭代步數(shù)的變化關(guān)系。圖中還給出了當(dāng)目標(biāo)邊長(zhǎng)為a=1.0λ的無(wú)限長(zhǎng)立方柱，上下層分界面均方根高度為0.1 λ時(shí)的情況。 由圖可知，對(duì)不同粗糙度的粗糙面，經(jīng)過15步迭代，外循環(huán)收斂精度都能達(dá)到10?5以下，遠(yuǎn)遠(yuǎn)滿足計(jì)算精度要求。經(jīng)過驗(yàn)證，通常收斂精度為τ=10?2即可滿足精度要求，因此經(jīng)過4步迭代即可終止，并且不論是圓柱目標(biāo)還是立方柱目標(biāo)該算法都具有良好的收斂性，同時(shí)由圖可知，由于立方柱周長(zhǎng)小于圓柱周長(zhǎng)，對(duì)于相同參數(shù)的粗糙面，目標(biāo)為方柱時(shí)的收斂性較好。

圖4 給出了內(nèi)循環(huán)FBM的迭代次數(shù)隨外循環(huán)迭代步數(shù)的關(guān)系，由圖可知，經(jīng)過4步外循環(huán)之后，F(xiàn)BM只需迭代1次即可收斂，具有很好的收斂速度。由于計(jì)算分層粗糙面與目標(biāo)復(fù)合電磁散射過程中，計(jì)算分層粗糙面表面電流占用了大量時(shí)間，因此FBM迭代次數(shù)的減少大大節(jié)省了時(shí)間，提高了計(jì)算效率。值得指出的是，當(dāng)目標(biāo)為圓柱，h1=h2=0.1λ，收斂精度取τ=10?2時(shí)，用互耦迭代算法實(shí)現(xiàn)一次需要240 s，而用傳統(tǒng)高斯消元法需要485 s，節(jié)省了約一半時(shí)間，實(shí)現(xiàn)50次可節(jié)省3個(gè)多小時(shí)，較傳統(tǒng)矩量法計(jì)算速度大為提高。

圖5給出了分層粗糙面中間存在無(wú)限長(zhǎng)金屬圓柱和立方柱時(shí)的散射系數(shù)，其中立方柱邊長(zhǎng)為a=1.0λ，圓柱半徑為1.0 λ。圖中還給出了只有分層粗糙面時(shí)的散射系數(shù)，由圖可知，當(dāng)粗糙面中間存在目標(biāo)時(shí)，后向散射方向散射系數(shù)明顯增大，而前向散射方向散射系數(shù)有所減小，并且當(dāng)目標(biāo)為圓柱時(shí)變化更為明顯，這是因?yàn)閳A柱目標(biāo)的周長(zhǎng)比立方柱稍大(2πR＞4a)，所以圓柱與分層粗糙面的作用更加明顯。因此，當(dāng)分層粗糙面中間存在目標(biāo)時(shí)，由于目標(biāo)與粗糙面之間的相互耦合作用，散射系數(shù)明顯發(fā)生改變，當(dāng)計(jì)算目標(biāo)與環(huán)境的相關(guān)問題時(shí)，必須考慮兩者之間的相互作用。

圖2 算法驗(yàn)證

圖3 τ隨迭代步數(shù)的變化關(guān)系

圖4 內(nèi)循環(huán)迭代次數(shù)

討論圓柱目標(biāo)大小對(duì)散射系數(shù)的影響。令其他參數(shù)不變，目標(biāo)半徑取0.5 λ，1.0 λ，1.5 λ，收斂精度取τ=10?2。雙站散射系數(shù)見圖6，由圖可知，隨著圓柱目標(biāo)半徑的增大，目標(biāo)與分層粗糙面的相互作用增強(qiáng)，復(fù)合散射系數(shù)增大，尤其在鏡面反射方向附近以外的很大范圍內(nèi)這一結(jié)論表現(xiàn)的尤為明顯?？梢妶A柱目標(biāo)的半徑大小對(duì)復(fù)合雙站散射系數(shù)有明顯影響，同時(shí)也驗(yàn)證了該算法對(duì)不同尺寸目標(biāo)均具有良好的收斂性。

圖7給出了當(dāng)R=1.0λ時(shí)，由文獻(xiàn)[4]中方法計(jì)算的單層粗糙面下方金屬圓柱目標(biāo)復(fù)合散射系數(shù)與分層粗糙面之間金屬圓柱目標(biāo)復(fù)合散射系數(shù)的對(duì)比結(jié)果。由圖可知，與單層粗糙面相比較，由于存在下層粗糙面與上層粗糙面的相互作用，以及下層粗糙面與目標(biāo)的相互作用，使得復(fù)合散射系數(shù)增大，尤其在鏡面反射方向附近范圍內(nèi)變化更為明顯。

圖5 有無(wú)目標(biāo)時(shí)分層粗糙面的散射系數(shù)

圖6 不同尺寸目標(biāo)對(duì)應(yīng)的散射系數(shù)

圖7 與單層粗糙面結(jié)果對(duì)比

討論目標(biāo)位置對(duì)散射系數(shù)的影響。其余參數(shù)不變，目標(biāo)深度取1.5 λ，2.0 λ和2.5 λ收斂精度取τ=10?2。相應(yīng)的散射系數(shù)如圖8所示，圖中還給出了目標(biāo)深度取1.5λ，2.0λ時(shí)單層粗糙面下方金屬圓柱目標(biāo)的散射結(jié)果。由圖可知，當(dāng)目標(biāo)深度增大時(shí)，分層粗糙面與下方目標(biāo)的復(fù)合散射系數(shù)變化不明顯，而單層面與下方目標(biāo)的散射系數(shù)明顯減小。這是由于對(duì)于雙層粗糙面，當(dāng)目標(biāo)深度增加時(shí)，目標(biāo)與上層粗糙面的距離增大，相互作用減小，但與下層粗糙面的卻又距離減小，相互作用增大，因此目標(biāo)與分層粗糙面的相互作用變化不是特別明顯。

4 結(jié)論

圖8 不同目標(biāo)深度對(duì)應(yīng)的散射系數(shù)

用基于FBM和Bi-CG的互耦迭代算法(CCIA)快速獲取了1維分層媒質(zhì)粗糙面之間無(wú)限長(zhǎng)金屬圓柱目標(biāo)和無(wú)限長(zhǎng)立方柱目標(biāo)的復(fù)合電磁散射特性，并對(duì)計(jì)算結(jié)果進(jìn)行了驗(yàn)證，討論了算法的收斂性，當(dāng)目標(biāo)周長(zhǎng)較小時(shí)，算法具有更好的收斂速度。計(jì)算了不同目標(biāo)尺寸對(duì)應(yīng)的復(fù)合散射系數(shù)，結(jié)果表明，隨著目標(biāo)尺寸的增大，分層粗糙面與目標(biāo)的相互耦合作用增強(qiáng)，散射系數(shù)增大，反之則減小。將計(jì)算結(jié)果與單層粗糙面下方金屬目標(biāo)結(jié)果相比較。結(jié)果表明，由于下層粗糙面與上層粗糙面以及目標(biāo)的相互作用，散射系數(shù)明顯增大。同時(shí)，討論了目標(biāo)深度變化對(duì)復(fù)合散射系數(shù)的影響，結(jié)果表明，當(dāng)目標(biāo)深度變化時(shí)，分層粗糙面與下方目標(biāo)的復(fù)合散射系數(shù)變化不明顯。

[1] Johnson J T. A study of the four-path model for scattering form an object above a halfspace[J]. Microwave Optical Technology Letters, 2001, 30(2): 130-134.

[2] Jin Y Q and Li Z. Simulation of scattering form complex rough surface at low grazing angle GFBM/SAA method [J].IEEE Transactions on Fundamentals and Materials Society(A), 2001, 121(10): 917-921.

[3] Ye Hong-xia and Jin Ya-qiu. Fast iterative approach to difference scattering from the target above a rough surface[J].IEEE Transactions on Geoscience and Remote Sensing, 2006,44(1): 108-115.

[4] Ye Hong-xia and Jin Ya-qiu. A hybrid analytical-numerical algorithm for scattering froma3-D target above a randomly rough surface[J]. Acta Physica Sinica, 2008, 57(2): 839-846.

[5] Wang X, Wang C F, and Gan Y B. Electromagnetic scattering from a circular target above or below rough surface[J]. Progress in Electromagnetics Research, 2003, 40:207-227.

[6] Wang Yun-hua, Zhang Yan-min, and Guo Li-xin.Investigation of the scattered field form a two-dimensional dielectric target above the planar surface with a Guass beam incidence[J]. Acta Physica Sinica, 2008, 57(9): 5529-5535.

[7] Wang Rui, Guo Li-xin, Qin San-tuan, and Wu Zhen-sen.Hybrid method for investigation of electromagnetic scattering interaction between the conducting target and rough sea surface[J]. Acta Physica Sinica. 2008, 57(6):3473-3480.

[8] Magda El-Shenawee. Polarimetric scattering from twolayered two dimensional random rough surfaces with and without buried objects[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(1): 67-76.

[9] Chew W C, Jin J M, Michielssen E, and Song J M. Fast and Efficient Algorithms in Computational Electromagnetics[M].Norwood: Artech House Publishers, 2001: 40-65.

[10] Li J, Guo L X, and Zeng H. FDTD Investigation on Electromagnetic Scattering form the PEC Cylinder above Two_Layered Rough Surfaces[C]. ICMMT, Nanjing 2008:531-534.

[11] Kuo Chih-hao and Moghaddam M. Electromagnetic scattering from a buried cylinder in layered media with rough interfaces[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(8): 2392-2401.

[12] Holliday D, DeRaad L L, and StCyrGJ. Forward backward:a new method for computing low-grazing angle scattering[J] .IEEE Transactions on Antennas Propagation, 1996, 44(5):722-729.

[13] Moss C D, Grzegorczyk T M, Han H C, and Kong J A.Forward-backward method with spectral acceleration for scattering from layered rough surfaces[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(3): 1006-1016.

[14] Lin J H and Chew W C. BiCG-FFT T-Matrix method for solving for the scattering solution from inhomogeneous bodies[J]. IEEE Transactions on Microwave Theory and Techniques, 1996, 44(7): 1150-1155.

[15] Thorsos A. The validity of the Kirchhoff approximation for toygh surface scattering using a Guassian roughness spectrum[J]. Journal of the Acoustical Society of America,1988, 83(1): 78-92.