郝靖瑋， 潘小敏， 盛新慶

(北京理工大學(xué) 信息與電子學(xué)院， 電磁仿真中心，北京 100081)

?

有限元中基于ID的快速掃角算法

郝靖瑋， 潘小敏， 盛新慶

(北京理工大學(xué) 信息與電子學(xué)院， 電磁仿真中心，北京 100081)

基于interpolative decomposition (ID)技術(shù)的有限元快速掃角算法，能有效地計(jì)算電目標(biāo)的單站雷達(dá)散射截面(RCS). 算法針對(duì)不同入射波對(duì)應(yīng)的右端項(xiàng)(RHSs)，構(gòu)成一個(gè)激勵(lì)矩陣. 將ID技術(shù)應(yīng)用到激勵(lì)矩陣用來(lái)選取對(duì)應(yīng)的Skeleton入射波的方向. 在采用快速算法得到Skeleton入射波對(duì)應(yīng)的解之后，所要求解的角響應(yīng)可以通過(guò)Skeleton對(duì)應(yīng)的解重構(gòu)出來(lái). 對(duì)電目標(biāo)的數(shù)值實(shí)驗(yàn)表明了該算法的效率和精度.

ID；Skeleton；角掃描；有限元

在許多電磁波散射的應(yīng)用中，需要計(jì)算大范圍角度的電磁場(chǎng)分布，本文稱這種計(jì)算為角響應(yīng). 每個(gè)入射角度對(duì)應(yīng)一個(gè)右端項(xiàng)，這就是說(shuō)每個(gè)入射角度都需要一次求解，因此對(duì)于計(jì)算大范圍角響應(yīng)的需求，計(jì)算量很大. 目前人們主要采用基于有理函數(shù)逼近的算法來(lái)實(shí)現(xiàn)大范圍角度響應(yīng)的計(jì)算. 有兩種方法來(lái)生成合理的近似函數(shù)：一種是插值模型，基于模型的參數(shù)估計(jì)(MBPE)[1]；另一種最先開(kāi)始應(yīng)用在高速電路分析，為漸進(jìn)波形估計(jì)(AWE)[2]. MBPE的關(guān)鍵在于建立優(yōu)良的擬合模型，但對(duì)于復(fù)雜目標(biāo)，很難找到一種通用的擬合模型. AWE固有的問(wèn)題在于其近似誤差不可控，并且對(duì)于大尺寸問(wèn)題，需要大量的插值點(diǎn)，而這些插值點(diǎn)的選取也是比較困難的. 為了解決這些問(wèn)題，本文基于Skeleton-based快速掃角理論，利用interpolative decomposi-tion (ID)[3]技術(shù)和Skeleton概念[4-7]來(lái)計(jì)算電目標(biāo)的RCS.

1 ID技術(shù)的簡(jiǎn)單介紹

首先對(duì)ID技術(shù)[3]做一個(gè)簡(jiǎn)單的介紹. 假定C是一個(gè)m×n的復(fù)矩陣，秩為q，q≤m，q≤n. 存在一個(gè)m×q的復(fù)矩陣B，B矩陣的列是C矩陣列的子集. 存在一個(gè)q×n型的矩陣P有以下性質(zhì)：

①P矩陣列的子集組成q×q的單位矩陣；

② P中元素的絕對(duì)值都不大于1；

④ P矩陣的最小奇異值至少為1；

基于上述性質(zhì)，近似的表達(dá)式可以寫成

(1)

當(dāng)Cm×n的秩大于q時(shí)，Cm×n的q+1階奇異值非常小. 在式(1)進(jìn)行分解之前，通常閾值εID用來(lái)控制近似誤差.εID在L2范數(shù)標(biāo)準(zhǔn)下，衡量C和B·P的誤差.

2 有限元方法

有限元方法(FEM)是一種非常有效的電磁計(jì)算方法，由于其使用靈活的四面體單元對(duì)計(jì)算區(qū)域進(jìn)行離散，因此該方法非常適合于求解復(fù)雜金屬和介質(zhì)目標(biāo)的電磁特性[8].

三維復(fù)雜目標(biāo)的散射屬于開(kāi)域問(wèn)題，使用有限元方法必須對(duì)無(wú)限區(qū)域進(jìn)行截?cái)? 本文采用一階吸收邊界條件(ABC)將無(wú)限區(qū)域截?cái)酁閹в羞吔鐥l件的有限計(jì)算區(qū)域Ω，在計(jì)算區(qū)域Ω內(nèi)散射電場(chǎng)滿足如下微分方程和邊界條件：

(3)

為使用有限元方法解決式(2)所定義的邊值問(wèn)題，整個(gè)計(jì)算區(qū)域需用四面體有限單元對(duì)其進(jìn)行離散，這樣計(jì)算區(qū)域內(nèi)的矢量電場(chǎng)可用邊緣元矢量基函數(shù)進(jìn)行展開(kāi). 使用矢量基函數(shù)對(duì)式(2)進(jìn)行迦略金匹配，可得到式(4)所示的有限元弱解方程為

(4)

有限元系數(shù)矩陣K和激勵(lì)向量f由下面公式計(jì)算得

K=

(5)

f=

(6)

式中N為包含邊緣元矢量基函數(shù)的列向量. 由式(4)(5)可以看出，對(duì)于同一散射問(wèn)題，有限元系數(shù)矩陣是相同的，對(duì)于不同的入射角度，激勵(lì)向量有所不同，因此得到不同的散射電磁場(chǎng). 通過(guò)求解方程(4)便可得到整個(gè)計(jì)算區(qū)域的未知電場(chǎng)系數(shù)，從而獲取整個(gè)電場(chǎng)分布.

3 快速算法基本思想

近期，Skeleton的概念已經(jīng)成功運(yùn)用到開(kāi)發(fā)快速方法. 在這里，它被用來(lái)實(shí)現(xiàn)大規(guī)模電磁波散射角響應(yīng)的快速計(jì)算. 對(duì)于有限元方法，將不同入射方向?qū)?yīng)的右端項(xiàng)和解組合起來(lái)，可以得到一個(gè)線性方程組，這些方程組可以寫成

(7)

式中：K為N×N的有限元矩陣；f(λ)為N×Ninc已知平面波的激勵(lì)矩陣；E(λ)為N×Ninc的未知求解矩陣；N為未知數(shù)的數(shù)量；Ninc為激勵(lì)的個(gè)數(shù)；λ為入射波的方向.

(8)

把ID技術(shù)應(yīng)用f(λ)，可以得到

(9)

fs(λ)為Skeletonized以后得到的激勵(lì)矩陣，R(λ)為投影矩陣. 把式(9)代入到式(8)可得

(10)

所以

(11)

是對(duì)應(yīng)于fs(λ)的解矩陣. 從式(10)和式(11)可知，E(λ)變?yōu)?/p>

(12)

通過(guò)式(11)求解出Es(λ)之后，所有入射波所對(duì)應(yīng)的解可由式(12)求出. 特別的，求解全部λ∈[λmin，λmax](λ=θ，φ)的角響應(yīng)，只需求解對(duì)應(yīng)于skeleton的激勵(lì)項(xiàng)，從而實(shí)現(xiàn)角響應(yīng)計(jì)算的加速.

由于ID的特性，矩陣fs(λ)中的每一列都包含在f(λ)中. 因?yàn)閒(λ)的每一個(gè)列都對(duì)應(yīng)一個(gè)激勵(lì)向量，也就是有限元每個(gè)入射角度對(duì)應(yīng)的右端項(xiàng)，fs(λ)的列則是Skeleton對(duì)應(yīng)的激勵(lì)向量. 通過(guò)式(11)可知，Es(λ)中的列數(shù)為Skeleton入射角度的個(gè)數(shù). 這個(gè)特點(diǎn)把Skeleton-based快速掃角算法和基于SVD/QR[9]的方法以及其擴(kuò)展方法[10]區(qū)分開(kāi)來(lái). 利用這種特性，可以方便地恢復(fù)所有入射波對(duì)應(yīng)的解. 具體的，如果式(12)中的Es(λ)和R(λ)都已計(jì)算出來(lái)，那么所有入射波的解可以直接從式(12)中恢復(fù)出來(lái).

4 數(shù)值算例

為了證明ID算法應(yīng)用到有限元的高效性和準(zhǔn)確性，采用如圖1所示的一個(gè)邊長(zhǎng)為0.5λ的介質(zhì)立方塊為例，分別研究介質(zhì)特性為各向同性和各向異性時(shí)算法的精度與效率. 計(jì)算中采用了239 066 個(gè)未知數(shù). 各向同性時(shí)，介質(zhì)參數(shù)為ε=2，μ=1. 各向異性時(shí)介質(zhì)參數(shù)為μ=1，

(13)

分別采用直接法和本文方法計(jì)算了介質(zhì)立方體在θ=0°～180°，φ=0°平面波入射下的單站雷達(dá)散射截面. 計(jì)算中角度的步進(jìn)為1°. 為了考察算法精度，采用

(14)

式中：β為絕對(duì)誤差；R2為直接方法計(jì)算結(jié)果；R1本文算法計(jì)算結(jié)果.

表1給出計(jì)算的統(tǒng)計(jì)情況，圖1和圖2分別為各向同性時(shí)雷達(dá)散射截面和誤差，圖3和圖4分別為各向異性時(shí)雷達(dá)散射截面和誤差. 計(jì)算中V(λ)由181個(gè)激勵(lì)向量組成，應(yīng)用Skeleton-based掃角方法，無(wú)論是各向同性介質(zhì)還是各向異性介質(zhì)，都只需選取16個(gè)激勵(lì)向量來(lái)生成Vs(λ). 而且由于ID非常高效，尋找Skeleton入射方向的時(shí)間開(kāi)銷很小，從而使得本文算法計(jì)算的加速到達(dá)7倍多. 由于式(1)的近似計(jì)算，所以該算法與直接法存在誤差，誤差可由閾值εID控制.

表1 兩種算法對(duì)比圖

從圖2和圖4的誤差看，算法的計(jì)算精度非常高. 雖然絕對(duì)誤差的最大值達(dá)到0.06 dBsm，但仔細(xì)觀察發(fā)現(xiàn)，對(duì)應(yīng)觀察角RCS的準(zhǔn)確值θRCS在-20 dBsm左右.

5 結(jié) 論

基于ID技術(shù)的有限元算法是一種計(jì)算電目標(biāo)單站雷達(dá)散射截面的有效算法. 為了獲得較高分辨率的雷達(dá)散射截面，往往需要對(duì)很多方向的入射波展開(kāi)計(jì)算. 不同入射波共用一個(gè)相同的有限元矩陣，把ID技術(shù)應(yīng)用于不同入射角度組成的激勵(lì)矩陣后，可以高效得到Skeleton入射方向. 一般情況下，只需較少的Skeleton入射方向，就能精確得到較好分辨率的雷達(dá)散射截面. 數(shù)值實(shí)驗(yàn)表明，該算法顯著提高了計(jì)算角效率，并且保證了計(jì)算的準(zhǔn)確性，該算法通過(guò)對(duì)不同入射角度的插值提取來(lái)加速計(jì)算，從而提高掃角的效率.

[1] Miller E K. Model-based parameter estimation in electromagnetics.iii. applications to EM integral equations[J]. IEEE Antennas Propag Mag， 1998,40(3):49-66.

[2] Pillage L T, Rohrer R A. Asymptotic waveform evaluation fortiming analysis[J]. IEEE Trans. Computer-Aided Design， 1990,9(4):352-366.

[3] Liberty E， WoolfeF， Martinsson P G， et al. Randomized algorithms for the low-rank approximation of matrices[J]. 2007,104:20167-20172.

[4] Ho K L, Greengard L. A fast direct solver for structured linear systems by recursive skeletonization[J]. SIAM J Sci Comput.， 2012,34(5):A2507-A2532.

[5] Pan X M, Sheng X Q. Hierarchical interpolative decomposition multilevel fast multipole algorithm for dynamic electromagnetic simulations[J]. Progr. in Electromagnetics Research， PIER， 2013,134:79-94.

[6] Pan X M, Sheng X Q. Preconditioning technique in the interpolative decomposition multilevel fast multipole algorithm[J]. IEEE Trans. Antennas Propag.， 2013,61(6):3373-3377.

[7] Pan X M， Wei J G， Peng Z， et al. A fast algorithm formultiscale electromagnetic problems using interpolative decompositionand multilevel fast multipole algorithm[J]. Radio Sci.， 2012,47:RS1011.

[8] Schroder A， Bruxns H D， Schuster C. A hybrid approach forrapid computation of two-dimensional monostatic radar cross section problems with the multilevel fast multipole algorithm[J]. IEEE Trans. Antennas Propag.， 2012,60(12):6058-6061.

[9] Zhen P， Stephanson M B， Lee J F. Fast computation of angular responses of large-scale three-dimensional electromagnetic wave scattering[J]. IEEE Trans. Antennas Propag.， 2010,58(9):3004-3012.

[10] 高紅偉，盛新慶.三維非均勻電大目標(biāo)散射的有限元撕接區(qū)域分解法計(jì)算[J].北京理工大學(xué)學(xué)報(bào)，2012，32(6):602-606.

Gao Hongwei， Sheng Xinqing. Computation of electromagnetic scattering by 3D large inhomogeneous targets utilizing finite element tearing and interconnecting method[J]. Transactions of Beijing Institute of Technology, 2012，32(6):602-606. (in Chinese)

(責(zé)任編輯：劉芳)

Fast Algorithm for Angular Sweeping in FEM Based on Interpolative Decomposition

HAO Jing-wei， PAN Xiao-min， SHENG Xin-qing

(Center for Electromagnetic Simulation， School of Information and Electronics，Beijing Institute of Technology， Beijing 100081， China)

A fast angular sweeping algorithm based on interpolative decomposition was proposed for the finite element method (FEM) to efficiently compute monostatic radar dispersion section (RCS). In the algorithm, an excitation matrix was constructed from multiple right-hand-sides (RHSs) associated with different incidents. Then, interpolative decomposition (ID) was applied to the excitation matrix to figure out the skeleton incidents. After the solutions corresponding to the skeleton incidents were obtained through FEM, the desired angular responses could be reconstructed from those of skeletons. The results of numerical experiments on isotropic and anisotropic dielectric targets show the efficiency and accuracy of the proposed algorithm.

ID; Skeleton; angular sweeping; finite element method(FEM)

2013-10-31

國(guó)家“九七三”計(jì)劃項(xiàng)目(2012CB720702)；國(guó)家自然科學(xué)基金資助項(xiàng)目(61371002)；國(guó)家教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃項(xiàng)目(NCET-12-0045)

郝靖瑋(1989—)，男，博士生，E-mail:hjw5299118@126.com；盛新慶(1968—)，男，教授，博士生導(dǎo)師，E-mail:xsheng@bit.edu.cn.

潘小敏(1978—)，男，副教授，博士生導(dǎo)師，E-mail:xmpan@bit.edu.cn.

O 441

A

1001-0645(2016)05-0498-04

10.15918/j.tbit1001-0645.2016.05.011